Solving the Boussinesq's (groundwater) equation by Cordano and Rigon

Boussinesq's equation describes the motion of a free aquifer under the Dupuit hypothesis, or, if you prefer, the motion of the surface of  the water table. This paper just submitted for reviewing treats its numerical solution using a new method derived from previous work of Brugnano and Casulli (2008) and Casulli (2009). Let's talk the abstract:


This work presents with a new conservative finite-volume numerical solution for the two-dimensional groundwater flow (Boussinesq) equation, which can be used for investigations of hillslope subsurface flow processes and simulations of catchment hydrology. The Boussinesq Equation is integrated for each grid element and can take account of the local variability of topography and soil properties within the grid elements. The numerical method allows for wetting and drying of the water-table, which has been successfully simulated.
The stability and convergence of the method is shown to be guaranteed apriori by the properties of the solver itself. 
The numerics are validated against some approximate analytical solutions, and  compared to another numerical solver of the Boussinesq equation. Finally, the solver capabilities are further explored with simulations of the Panola experimental hillslope where  the bedrock topography, which is accurately known, causes complex wetting and drying patterns;  in this situation the importance of a two-dimensional description of subsurface flows to obtain properly simulated discharges becomes clear.

What is relevant in our paper, among other thing, is:

• the cleanliness of the integration method (that makes the numerics fast)
• the ability to deal with wetting and drying zones
• the accurate  treatment of boundary conditions (which could, in theory, break the conditions of integrability a-priori)

The first appendix can be used as an introduction to the Newton's method for integration of non linear equations.

As usual, the code we provide is released with a GPL v3 license, and available for download, at the moment, in the GEOtop site. Emanuele Cordano wrote also some ancillary R code which plots the analytical solutions of the equation in some simplified setting described in the paper (under Packages - boussinesq in the CRAN site).

The model was used to create the reference conditions in the paper by Lanni et Al., 2011 [pdf]. The draft of the paper for curious can be downloaded from here. The final version can be retrieved through this post.