Combining mesh-less finite difference method and least square approximation, a new numerical model is developed for highly dispersive and fully nonlinear Boussinesq equations in two horizontal dimensions. The 3 rd order truncated series solution of the Laplace equation is employed to approximate the velocity distribution in the vertical plane. The linear properties of the wave model are discussed with Fourier analysis. It is shown that the model is suitable to predict the propagation of water waves at the range of 0 ≤ kh ≤ 10 for both the linear dispersion characteristic and shoaling gradient. Preliminary verifications of the numerical model are given for nonlinear wave shoaling problems, wave run-up on conical island. The numerical results agree well with the experimental data available in the literature.