Very generally speaking, the Method of Lines approach is to first spatially discretize the PDEs (in this case, the Boussinesq equations) and then use a time-integrator to advance the resulting system of ODEs. The choice of discretization method depends on the geometry and the Prandtl number of the Boussinesq equations. For example, people use finite differences and even Chebyshev pseudospectral methods for this purpose. However, the equations may have complicated boundary conditions, in which case people use things like the Influence Matrix method.

For solving Boussinesq equations with an infinite Prandtl number in a spherical shell, see Yoshida, M., and A. Kageyama (2004), Application of the Yin‐Yang grid to a thermal convection of a Boussinesq fluid

with infinite Prandtl number in a three-dimensional spherical shell. See also this paper by Wright, Flyer and Yuen which combines RBF-PS and Chebyshev-PS to accomplish this (http://scholarworks.boisestate.edu/cgi/viewcontent.cgi?article=1027&context=math_facpubs)

Very generally speaking, the Method of Lines approach is to first spatially discretize the PDEs (in this case, the Boussinesq equations) and then use a time-integrator to advance the resulting system of ODEs. The choice of discretization method depends on the geometry and the Prandtl number of the Boussinesq equations. For example, people use finite differences and even Chebyshev pseudospectral methods for this purpose. However, the equations may have complicated boundary conditions, in which case people use things like the Influence Matrix method.

For solving Boussinesq equations with an infinite Prandtl number in a spherical shell, see Yoshida, M., and A. Kageyama (2004), Application of the Yin‐Yang grid to a thermal convection of a Boussinesq fluid

with infinite Prandtl number in a three-dimensional spherical shell. See also this paper by Wright, Flyer and Yuen which combines RBF-PS and Chebyshev-PS to accomplish this (http://scholarworks.boisestate.edu/cgi/viewcontent.cgi?article=1027&context=math_facpubs)

There are many articles and theses which discuss numerical solution of Boussinesq equations. Please see the attachments.

http://www.hach.ulg.ac.be/cms/system/files/2004%20-%20J%20Hydraul%20Eng%20ASCE%20-%20Mohapatra%20Chaudhry%20-%20Numerical%20Solution%20of%20Boussinesq%20Equation%20to%20simulate%20Dam-break%20flow.pdf

http://www.cs.toronto.edu/pub/reports/na/heso05.pdf

http://cedb.asce.org/cgi/WWWdisplay.cgi?140062

http://etheses.whiterose.ac.uk/1291/1/walkley.pdf

Article Numerical Solution of Boussinesq Equations as a Model of Int...

Article Numerical studies on Boussinesq-type equations via a split-s...

The Boussines equations are implemented in the OpenFoam software packeage, which is open source. With some experience in CFD you will get very fast nice results. 

J.M. Sanz-Serna and his coauthors have several papers on numerically solving the Boussinesq equation (1988, 1990 and 1991).

Thank you very much.

An efficient numerical technique is used to obtain solutions of the Boussinesq equation for problems of groundwater recession and groundwater flow in response to changes in stream stage. Comparison of the numerical results with the analytical solution of Boussinesq [1904] attests to the accuracy of the former. The solution of a recession problem indicates that the form of the Werner and Sundquist [1951] model of groundwater discharge recession is appropriate. The solution of a problem involving changing boundary conditions on an aquifer because of a flood wave provides data relative to groundwater outflow. A comparison with the results of Cooper and Rorabaugh [1963] evaluates the applicability of a linear model of unconfined flow.

sincerey

Farzad

Thank you very much.