The difference between semi-continuum model and Richards’ equation for unsaturated porous media flow

Semi-continuum modelling of unsaturated porous media flow is based on representing the porous medium as a grid of non-infinitesimal blocks that retain the character of a porous medium. This approach is similar to the hybrid/multiscale modelling. Semi-continuum model is able to physically correctly describe diffusion-like flow, finger-like flow, and the transition between them. This article presents the limit of the semi-continuum model as the block size goes to zero. In the limiting process, the retention curve of each block scales with the block size and in the limit becomes a hysteresis operator of the Prandtl-type used in elasto-plasticity models. Mathematical analysis showed that the limit of the semi-continuum model is a hyperbolic-parabolic partial differential equation with a hysteresis operator of Prandl’s type. This limit differs from the standard Richards’ equation, which is a parabolic equation and is not able to describe finger-like flow.

The zero-flux boundary condition at the bottom boundary is not necessary and can be replaced by another one without a significant effect on the limit equation. Let us denote by x i+1 the point on the boundary between blocks i and i + 1. Points x 0 and x N+1 form the boundary of the one-dimensional sample. We put h = ∆x. The i−th block thus corresponds to the interval [x i , x i + h). Furthermore, we define the following notation. A functionf h (x) denotes a piecewise constant function that takes the value f i on each interval x ∈ [x i , x i + h). In this notation,S h (x,t) is the saturation, piecewise constant on each block. A piecewise linear (and thus continuous on [0, L]) approximation will be denoted without the tilde, i.e.
In the following paragraphs, we will use these piecewise constant and piecewise linear approximations to saturation, pressure, and flux. Equation (3) can be rewritten as a partial differential equation where we used This approach is inspired by Rothe's scheme which relates a partial differential equation to its discrete approximation 1 . Similarly, equation (2) can be rewritten as follows: where P h is the hysteresis operator generated by equations (4) and (5). Next, we employ the toolbox of modern partial differential equation theory and pass to the weak formulation. We take a smooth function ϕ with compact support, multiply both sides of equation (8) by it, and integrate over the interval [0, L] to get In this equation, the derivative is transferred to the test function ϕ which enables us to relax the assumptions on the limits performed below. The retention curve defined by equations (4) and (5) can be decomposed into two branches (wetting and draining) in the following symmetrical way: The scanning curves are modelled as line segments with the derivative K PS . Let us note that the constant K PS is fixed during the limiting process. Using the form of the two branches, we can see the similarity between equation (11) and the classical Prandtl model of elasto-plasticity (the stop operator), see equations (2.2) − (2.5) (p. 15 − 16) in Visintin for more details 2 . The hysteresis operator P(h) defined by equation (11) and the scanning curves can be thus expressed using the differential inequality To see the similarity with the classical Prandtl model of elasto-plasticity, we use the substitution The resulting operator form of the semi-continuum model takes the form of the following differential inequality We now check whether the differential inequality in equation (13) corresponds to the hysteresis operator defined by equation (11) with the scanning curves: where t 0 is an initial time such that P(h,t 0 ) ∈ [P d (h, S(t 0 )), P w (h, S(t 0 )]. Thus, we are located on the scanning curve for this case.

Let
Hence and from equation (12) it follows because the pressure corresponds to this branch only if ∂ t S ≥ 0 and moreover ∂ t f (h, S) = ∂ S f (h, S)∂ t S. Without the inequality ∂ t S ≥ 0 we are unable to connect the scanning curves with the wetting branch. From the inequality K PS ≥ ∂ S f (h, S) it follows that the scanning curves and the wetting branch are connected.
because the pressure corresponds to this branch only if ∂ t S ≤ 0 and moreover ∂ t f (h, S) = ∂ S f (h, S)∂ t S. Without the inequality ∂ t S ≤ 0 we are unable to connect the scanning curves with the draining branch. From the inequality K PS ≥ ∂ S f (h, S) it follows that the scanning curves and the draining branch are connected.
Finally, assume that as h → 0, the solutions of the model derived above and the corresponding operator converge in the following senseS The validity of this assumption is suggested by the numerical evidence described in the main article. Assume further we can apply Lebesgue's dominated convergence theorem to equation (10). Performing the limit in equations (10) and (11) yields where and Thus, the limit form of the semi-continuum model is a weak formulation (15) for a partial differential equation together with the classical Buckingham-Darcy law (16), containing a hysteresis operator of the Prandtl model of elasto-plasticity (17). Passing from the weak formulation to the classical one yields Equations (17)-(18) represent the classical form of the limit of the semi-continuum model. It is a partial differential equation with a hysteresis operator of Prandtl-type under the derivative. If we are located on the main wetting or draining branches, the limit will be a hyperbolic differential equation, because the pressure saturation relation is constant and thus independent on the space variable. This makes the limit switch between parabolic and hyperbolic type. Let us note that in the case of continuous saturation, it applies k(S) = k(S − ) k(S + ).