  Richards equation is often used to represent two-phase fluid flow in an
unsaturated porous medium when one phase is much heavier and more viscous than
the other. However, it cannot describe the fully saturated flow due to
degeneracy in the capillary pressure term. Mathematically, gravity-driven
variably saturated flows are interesting because their governing partial
differential equation switches from hyperbolic in the unsaturated region to
elliptic in the saturated region. Moreover, the presence of wetting fronts
introduces strong spatial gradients often leading to numerical instability. In
this work, we develop a robust, multidimensional mathematical and computational
model for such variably saturated flow in the limit of negligible capillary
forces. The elliptic problem for saturated regions is built-in efficiently into
our framework for a reduced system corresponding to the saturated cells, with
the boundary condition of the fixed head at the unsaturated cells. In summary,
this coupled hyperbolic-elliptic PDE framework provides an efficient,
physics-based extension of the hyperbolic Richards equation to simulate fully
saturated regions. Finally, we provide a suite of easy-to-implement yet
challenging benchmark test problems involving saturated flows in one and two
dimensions. These simple problems, accompanied by their corresponding
analytical solutions, can prove to be pivotal for the code verification, model
validation (V&V) and performance comparison of such simulators. Our numerical
solutions show an excellent comparison with the analytical results for the
proposed problems. The last test problem on two-dimensional infiltration in a
stratified, heterogeneous soil shows the formation and evolution of multiple
disconnected saturated regions.
