AbstractThis paper investigates preconditioned iterative techniques for finite difference solutions of a high‐order Boussinesq method for modelling water waves in two horizontal dimensions. The Boussinesq method solves simultaneously for all three components of velocity at an arbitrary z‐level, removing any practical limitations based on the relative water depth. High‐order finite difference approximations are shown to be more efficient than low‐order approximations (for a given accuracy), despite the additional overhead. The resultant system of equations requires that a sparse, unsymmetric, and often ill‐conditioned matrix be solved at each stage evaluation within a simulation. Various preconditioning strategies are investigated, including full factorizations of the linearized matrix, ILU factorizations, a matrix‐free (Fourier space) method, and an approximate Schur complement approach. A detailed comparison of the methods is given for both rotational and irrotational formulations, and the strengths and limitations of each are discussed. Mesh‐independent convergence is demonstrated with many of the preconditioners for solutions of the irrotational formulation, and solutions using the Fourier space and approximate Schur complement preconditioners are shown to require an overall computational effort that scales linearly with problem size (for large problems). Calculations on a variable depth problem are also compared to experimental data, highlighting the accuracy of the model. Through combined physical and mathematical insight effective preconditioned iterative solutions are achieved for the full physical application range of the model. Copyright © 2004 John Wiley & Sons, Ltd.