The coupled equations for flow in unsaturated soil as proposed by Beliaev and Hassanizadeh  are described. These equations account for mechanism of dynamic capillary pressure via a first order relaxation function. A form of the relaxation function for the dynamic capillary pressure-saturation relation is proposed based on physical reasoning and a semi-analytical solution to the flow equations. A mass conservative and computational efficient numerical solution to the coupled equations in two space dimensions is derived and applied to the simulation of gravity-driven unstable flow. Simulated fingers have all the morphological features of fingers observed in laboratory experiments. The results demonstrate that the dynamic capillary pressure mechanism causes initial destabilization of the flow, while the mechanism of capillary hysteresis leads to finger persistence.
Bibliographical noteFunding Information:
There are two main requirements for modeling the phenomenon of gravity-driven fingering in unsaturated porous media: (i) the mathematical model must be capable of producing unstable perturbations, and (ii) the growing perturbations have to be persistent. In our previous work  a linear stability analysis was applied to the basic traveling wave solution, and we studied the conditions for the growth of perturbations. We analyzed three distinct models: (i) the conventional Richards equation (RE), (ii) a sharp front Richards equation (SFRE) , and (iii) an extended Richards equation with a nonequilibrium pressure-saturation function (RRE) . The analysis by Egorov et al.  shows that among the three models considered the RRE model is the only model capable of generating a structured field of fingers from a slightly perturbed uniform wetting flow. It was shown in [4,5] that persistence of fingers is dominated by hysteresis in the capillary pressure-saturation relation. Therefore, the hysteresis must be incorporated in the model. Beliaev and Hassanizadeh  recently extended the RRE model to what we refer to as the HRRE model, by accounting for hysteresis. In this paper we formulate the HRRE model, develop a numerical solution scheme, and then apply it to two-dimensional simulation of finger propagation. Two cases are considered, the case of a single finger *The authors wish to acknowledge support under NATO Collaborative Linkage Grant 978242. This work was supported in part by the Army High Performance Computing Research Center under the auspices of the Department of the Army, Army Research Laboratory cooperative agreement number DAAD19-01-2-0014, the content of which does not necessarily reflect the position or the policy of the government, and no official endorsement should be inferred.