Richards' equation reconsidered

Semi-continuum modelling of unsaturated porous media ﬂow is based on representing the porous medium as a grid of non-inﬁnitesimal blocks that retain the character of a porous medium. Semi-continuum model is able to physically correctly describe diffusion-like ﬂow, ﬁnger-like ﬂow, 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 partial differential equation with a hysteresis operator of Prandl’s type. This limit differs from the standard Richards’ Equation (RE), which is not able to describe ﬁnger-like ﬂow. Since the physics behind both RE and the semi-continuum model is almost the same, we suggest a way to reformulate the RE so that it retains the ability to describe ﬁnger-like ﬂow. We conclude that RE should be reconsidered by means of appropriate modelling of the hysteresis and correct scaling of the retention curve.


Introduction
A porous medium is a solid material which contains pores. Pores are filled with one or more different fluids -liquids and gases. The skeletal structure of a porous medium (called the matrix) has usually a very complex geometry. In soil physics, the porous medium is called saturated if all the pores contain a liquid (usually water), and unsaturated if some pores are filled with liquid and some with gas (usually air). Saturation is defined as the ratio of the liquid volume and the pore volume. There are many important applications of porous media flow modeling, such as oil recovery 1 and soil physics 2,3 .
In contrast to "bulk" fluid flow (described e.g. by the Navier-Stokes equations), inertial forces in the porous medium can usually be neglected due to small flow velocities. However, the porous media flow may still be very complicated because the liquid movement in the porous matrix depends on the pressure in the wetting fluid, which is determined by the capillary forces 4,5 which, in turn, are directly affected by the geometry of the porous matrix 6 .
In an unsaturated porous medium, capillary pressure and saturation are related by a material characteristics -the so-called retention curve -which is known to exhibit substantial hysteresis 7 . The retention curve gives the matrix potential (or capillary pressure) as a function of saturation under equilibrium conditions 8 . It was shown, that the retention curve depends on the volume of the sample 9-12 . However, this important issue has not usually been considered in porous media flow modeling.
Three different approaches to porous medium flow modeling can be identified: (1) Continuum modelling, which is based on the idea of representing the porous medium as a continuum using the concept of the representative elementary volume (REV) 13 . In this concept, the key physical quantities (saturation and capillary pressure) are considered scalar fields smoothly varying in (continuous) time and space. (2) Percolation theory, in which the pore space is described as a network of nodes and bonds 14 . All physical quantities (saturation, pressure, and time) are usually considered discrete in these models. (3) Semi-continuum modelling, in which the porous medium is represented as a grid of non-infinitesimal blocks that retain the character of a porous medium 15,16 . Saturation and capillary pressure are assumed continuous but constant within each block, the time is either continuous or discrete.
In soil physics, three different descriptions of liquid transport through the porous material have been used. (A) At the continuum and semi-continuum level, the movement of a liquid through the porous body has usually been described by the Darcy-Buckingham law 17 . (B) At the level of capillary pores, various displacement rules have been proposed [18][19][20] . (C) At the level of large, non-capillary pores, the kinematic wave equation has been used 21 .
As a result, three different modelling approaches (1-3) may be combined with at least three flow rules (A-C). Thus, a plethora of flow models are used to describe fluid flow in various types of porous media 14,[22][23][24][25][26] . The oldest model is a combination (1A), which is expressed by the Richards' equation (RE) 27 . It has long been known that this model works well only for diffusion-like flow and fails to describe finger-like flow 28,29 . Therefore, especially in the field of soil physics and soil hydrology, many attempts were made to modify the RE so that it may also capture finger-like flow [30][31][32][33][34][35][36][37][38][39][40] . However, none of these modifications is able to reproduce all the observed features of finger-flow 2 . A different attempt to overcome these drawbacks is the combination (3A) proposed in 15,16,41,42 . It has been shown that some of these models describe diffusion-like flow, finger-like flow, and the transition between the two regimes 16,42 .
Models (1A) and (3A) are very similar -they both use a continuous method of modelling the discrete porous material and they both use an almost identical form of the Darcy-Buckingham law to describe liquid transfer. Therefore, the question arises: Why do they differ in their ability to model different types of flows? This article seeks answers to this question.

Continuum modelling of porous media flow -Richards equation
The Darcy-Buckingham law is the key constitutive relationship for modelling flow and transport in the saturated and unsaturated porous medium 13 . The flux q [m/s] modeled by the Darcy-Buckingham law takes the form: where κ [m 2 ] denotes the intrinsic permeability, ρ [kg/m 3 ] the fluid density, g [m/s 2 ] acceleration due to gravity, µ [Pas] the dynamic viscosity of fluid, and P [Pa] the capillary pressure in the unsaturated medium (caused by capillarity) and the hydrostatic pressure in the saturated medium (caused by gravity). There are two material characteristics in the unsaturated medium, the retention curve and the function k(S) [−], which is called the relative permeability. Both these material characteristics have to be measured. Many different models of the retention curve and the relative permeability function exist 13,43 . RE is then obtained by the combination of the Darcy-Buckingham law and the mass balance law equations 13 and is usually stated in the following form in 1D: where θ [-] denotes the porosity of the material. We use the notation From the mathematical point of view, RE is an elliptic partial differential equation of the second order for a saturated porous medium (P is the hydrostatic pressure in this case) and is a parabolic partial differential equation of the second order for an unsaturated porous medium 13 . In summary, RE consists of the law of mass conservation, together with a constitutive relationship for the liquid flow, and two material characteristics.

Semi-continuum modelling of porous media flow
The semi-continuum modelling of fluid transfer in porous media is based on representing the porous medium as a grid of non-infinitesimal blocks that retain the character of a porous medium. Many soil science researchers have tested this idea. As soon as 1989, Glass and Yarrington 15 proposed a cellular automaton under the title of "mechanistic modelling", or "Macro Modified Invasion Percolation". In Kmec et al. 16,42 , the authors introduced a semi-continuum model which will be described in this section. Let us consider a long narrow vertical test tube of cross-section A [m 2 ] and height L [m] filled with homogeneous and isotropic porous medium. We divide the tube into blocks of height dx. These blocks represent "pieces" of the original porous medium in the sense that each block is characterized by its retention curve, porosity, and permeability. The tube now consists of slices [idx, (i + 1)dx] with i = 0, 1, . . . , N. The key quantities that we want to track are: • Saturation S [-] in each block.
• Capillary pressure P [Pa] in each block.
• Fluxes q i, j [m/s] between the blocks i and j. Volumetric fluxes in [m 3 /s] can be recovered as Aq i, j . Naturally, only fluxes between pairs of neighboring blocks are nonzero in this setting.
At each instant, the pressure P = P(x,t) and saturation S = S(x,t) are considered constant within each block. We thus use the notations S i (t) = S(x,t) and P i (t) = P(x,t), x ∈ [idx, (i + 1)dx). Gravity is directed downward along the long axis of the tube, which is called the x-axis here. A constant (in time) influx q 0 [m/s] is set across the top boundary (at x = 0). Zero discharge q L = 0 (i.e. zero flux) is assumed at the bottom boundary (at x = L). Saturation and pressure in each block, and the fluxes across the block boundaries are considered continuous in time, however, to solve the model numerically, time is discretized with a step dt. At each time step, the saturation in each of the blocks is updated according to equation (4), which is an explicit discretized version of the semi-continuum model (equation (3)) A backward time discretization of equation (3) can be also used (see the Discussion in 42 ). Next, the capillary pressure in each block is updated according to the retention curve. The hysteresis is modelled by the simplest approach possible: If a block switches from draining to wetting, the capillary pressure starts moving from the draining branch toward the main wetting branch of the retention curve along a straight line with a very large gradient K PS 7,44,45 . Once the block (now in wetting mode) reaches the main wetting branch, it sticks to it and continues along it (see Fig. 1). This approach to hysteresis is motivated by the Prandtl-type hysteresis operator 46 . Modelling the scanning curves as lines with large slopes is a numerically feasible realization of the idea that the pressure in the block "jumps" from the draining branch to the wetting branch without any accompanying change in saturation. All it takes to produce a large increase in the liquid pressure (i.e., a decrease in matrix suction) is to relax (i.e., increase the radius of) the capillary menisci supporting the water body. This may be achieved with a negligible amount of liquid, thus keeping the saturation almost unchanged. This mechanism explains the hysteresis in the retention curve and offers the above-described modelling strategy. The last step of the modelling loop is the update of the fluxes. In the fingering regime, it is crucial how the model treats the conductance at the finger tip, where ∇P is large, and S changes abruptly from small values in front of the fingertip to large values inside the finger. The semi-continuum model uses the following discretization of the Darcy-Buckingham law: Thus, for the relative permeability across the boundary of two adjacent blocks, we simply take the geometric mean of the permeability of the respective blocks. The geometric mean √ ab has the desirable property of being small if one of the numbers a or b is small. It is also possible to use the harmonic mean with similar results. There is a theoretical justification for using this type of averaging 47 .
The semi-continuum model 16 works as follows: 1. The size of the blocks dx is chosen and an appropriate time step dt is set. Initial saturation is prescribed in each block, the corresponding capillary pressure is computed from the retention curve, and all fluxes are initially set to zero.

3/12
2. The top boundary condition is set: The flux into the topmost block is set and fixed to q 0 . The bottom boundary condition is set: The flux out of the bottom block is set and fixed to zero.
3. Using the current value of the fluxes q i, j , saturation S i in each block is updated according to equation (4).
4. Pressure P i in each block is updated according to retention curve and hysteresis model, keeping track whether the block is in the imbibition or draining mode.
5. Fluxes q i, j between neighboring blocks are updated according to equation (5), keeping the boundary fluxes fixed by step 2.
6. Time is updated to t + dt and the process goes back to step 3.

Scaling of the retention curve
The semi-continuum model uses only the physics of RE and in the limit dx → 0, it would reduce to RE. However, in this limiting process, the dependence of the retention curve on the block size 9, 11 would be missed. We argue that this is the key ingredient which is usually missing in the RE-based models and their modifications. Therefore, in the semi-continuum model, the dependence of the retention curve on the block size dx is implemented. This procedure is named scaling of the retention curve.
This procedure must respect a physically justified requirement that the nature of the flow does not change when the block size is changed. This means that the fluxes across the block boundaries must stay roughly the same when dx changes. The fluxes are given by equation (5) in which decreasing dx by half increases the flux by a factor of two. To compensate for this, the difference in pressure between the blocks must decrease. Decreasing the pressure difference without changing the saturation difference amounts to the flattening of the retention curve. Based on this idea, we propose the following simple scaling mechanism in which the main branches of the retention curve of a block take the form for the main wetting branch, and for the main draining branch, where C 1 [Pa] and C 2 [Pa] are constants. The parameter h [m] is the scaling parameter equals to the block size dx. Notice that in the scaling process, the distance between the two main branches does not change (although this feature is not a crucial requirement for further consideration). More standard models of the retention curve may also be used in the scaling process, e.g. the van Genuchten equation 43 . Let us explain, by means of an example of a porous material with a simple pore structure, how the shape of the retention curve changes during scaling. The retention curve given by equations (6) and (7) with C 1 = −700 Pa and C 2 = −1300 Pa roughly matches the main branches of retention curve of 20/30 sand in the experiments of DiCarlo 48 . Figure 2 illustrates the scaling of this retention curve as dx goes to zero. As dx decreases, the slope of both branches of the retention curve decreases, too. In the limit, both branches take the form of horizontal lines.
Let us show that the proposed scaling of the retention curve indeed does not change the nature of the flow. A simulation of liquid infiltration into a vertical column of porous material with a constant influx q 0 is shown in Fig. 3. The parameters used for the simulations are given in Table 1. Left panel of Fig. 3 shows the calculated moisture profiles for the initial saturation S i,in = 0.01 for a decreasing sequence of block size values. Right panel of Fig. 3 shows the numerical convergence for dx → 0 in case of a diffusion-like flow regime (the initial saturation S i,in = 0.14). We may observe that this retention curve scaling preserves the character of the flow across all levels of dx both in the finger regime and the diffusion-like regime. Therefore, the semi-continuum model allows for a physically reasonable scaling of the retention curve. Note that such considerations are impossible when deriving the RE. This numerical evidence suggests there should be a limit form of the semi-continuum model, i.e. a model should exist to which the semi-continuum model converges as dx → 0, if the retention curve is scaled in this appropriate way.

Limit of the semi-continuum model
With this scaling of the retention curve, we can finally introduce the limit of the semi-continuum model. The limit of the model derived here is a formal one, i.e. we assume the solution of the semi-continuum model converges as dx → 0 to a function and show which equation this function should satisfy. This equation will be called called the limit of the semi-continuum model. First, let us complete the model (equations (4), (5)) by stating the initial conditions and boundary conditions q 0,1 (t) = q 0 > 0 and q n,n+1 (t) = 0.
Let us denote by x i+1 the point on the boundary between blocks i and i + 1. Points x 0 and x n+2 form the boundary of the one-dimensional sample. We put h = dx. 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 49 . Similarly, equation (5) can be rewritten as follows: 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 (10) 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 (6) and (7) can be decomposed into two branches (wetting and draining) in the following symmetrical way: and C = C 1 −C 2 2 . 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 (13) 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 46 . The hysteresis operator P(h) defined by equation (13) and the scanning curves can be thus express 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 (15) corresponds to the hysteresis operator defined by equation (13) with the scanning curves: Integration in time variable together with equation (14) results in 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.  (14) 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.
Hence and from equation (14) 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 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 operator form of the semi-continuum model (equation (15)) converge in the following sensẽ The validity of this assumption is suggested by the numerical evidence described above. Assume further we can apply Lebesgue's dominated convergence theorem to equation (12). Performing the limit in equations (12) and (13) yields where and Thus, the limit form of the semi-continuum model is a weak formulation (17) for a partial differential equation together with the classical Buckingham-Darcy law (18), containing a hysteresis operator of the Prandtl model of elasto-plasticity (19). Passing from the weak formulation to the classical one yields Equations (19)- (20) represent the classical form of the limit of the semi-continuum model. It is a partial differential equation containing a Prandtl-type hysteresis operator under the derivative. If we are located on the main wetting or draining branches, the limit will be a hyperbolic differential equation. In this case, the pressure saturation relation is constant and thus independent on the space variable. This makes the limit switch between parabolic and hyperbolic type. Similar notation as in equation (20) is also used for the classical RE (see equation (21)).
The two equations differ only in how the relationship between pressure and saturation is expressed. In the limit of the semicontinuum model (equation (20)), the pressure-saturation relation is defined by a Prandtl-type hysteresis operator P H , while for the RE (equation (21)), the pressure-saturation relation is described by a hysteretic operator P with wetting and draining branches represented by monotonically increasing functions. However, in a case of constant boundary influx q 0 , the hysteresis does not have any effect on the final solution of the RE, because the saturation is always a nondecreasing function of time 50 . Thus we are not able to reach the scanning curve of the retention curve. In light of the analysis shown above, it is clear that the RE arises from a physically unsound limit of the semi-continuum model.

Discussion
The limit of the semi-continuum model (equations (3), (4), (5)) was found by means of mathematical analysis (equations (8) to (19)) in the form of a partial differential equation (20). However, the limiting process is inspired by a numerical consideration that the flow between adjacent blocks should remain roughly the same when the block size decreases. From the Darcy-Buckingham equation (5) it follows that the shape of the retention curve must be scaled according to the block size (equations (6), (7)). The process of scaling the retention curve is not a common practice in flow modeling. Therefore, its physical justification should be addressed. The question is whether the shape of the retention curve depends on the volume of the sample. This issue has been examined in the literature for a long time. The discrepancy between retention curve models and the actual measurements was already mentioned more than 60 years ago by Fatt 51 . The effects of sample dimension on capillary pressure have since then been pointed out e.g. in [9][10][11][52][53][54] . In our view, the explanation is quite simple: The main draining branch is usually measured in a pressure plate extractor. In the case of the pressure plate test 55,56 , water is extruded by gas. During the test, gas first invades large pores and displaces water there. With increasing pressure, the gas enters smaller and smaller pores gradually 55,56 . However, the topology of the porous medium is usually ignored although it plays a crucial role. During drainage, certain pores are "candidates for draining" because the applied air pressure is greater than the capillary pressure holding the water inside the pore. However, some of these pores may not be accessible because they are not connected to the body of advancing air, or water is not able to leave them because they are not connected to the body of retreating water 11 . Pražák et al. 55 used percolation theory to model a simple capillary network. He showed that the retention curve is a constant for a homogeneous network (i.e. a network with zero variability in pore size). There is an analogy between the proposed scaling of the retention curve and such a homogeneous network. By decreasing the sample size, the variability of the pore sizes inside the sample also decreases (i.e, is more homogeneous) and thus the retention curve becomes flatter. This was experimentally confirmed in 12 .
Let us perform the following thought experiment: Imagine a single pore consisting of a cylinder of radius R. Assuming zero contact angle, a drop of liquid "sit" inside this capillary cylinder bounded by two hemispherical menisci of radius R. According to the Young-Laplace equation 57 , the pressure drop across both menisci is 2σ /R, where σ is the surface tension of the liquid. Setting the gas pressure to zero, the liquid drop is under tension (i.e. negative pressure) P = 2σ /R. Connecting such an empty pore to a liquid reservoir at a pressure lower than −P will yield zero saturation in the pore -the suction of the pore is not enough to draw any liquid inside. Once the pressure in the reservoir increases above −P, the pore will immediately fill with the liquid switching from zero to unit saturation. Thus, the dependence of saturation on pressure (i.e. the retention curve) is a horizontal line at −P. Let us continue this thought experiment and consider two pores of radii R 1 < R 2 . At a certain pressure −P 1 , the first pore will fill in completely, and at a higher pressure −P 2 , the second one will fill in. Thus, the retention curve of this pair of two pores becomes a broken horizontal line (see Fig. 4). A macroscopic sample of a porous medium contains many pores of various shapes. The resulting retention curve of the sample arises by assembling many horizontal lines at different levels of pressure. The main point of this excursion is to explain that as the sample size converges to zero, its retention curve has to converge to the retention curve of a single pore -i.e. to a horizontal line. It follows from the above considerations that the retention curve depends on the size of the sample. Therefore, the scaling we propose is supported by sound physical arguments. In continuum physics, material characteristics, such as the retention 8/12 curve, are related to the volume of REV. Thus, if blocks smaller than REV are used, the retention curve must be scaled. This causes substantial problems in the mathematical setting. The reconsidered Richards' equation switches between a hyperbolic equation (at points located on the main branches) and a parabolic equation (at points located on the scanning curves). It is the hysteresis operator in the limiting equation that enables to switch between two equations of different types. The situation is moderately similar to modelling the compressible fluid flow in subsonic and transsonic regions, which changes from hyperbolic to elliptic equation 58 . However, in the case of the reconsidered Richards' equation, the hysteresis operator appears inside the Laplace operator. Such equations seem to be unexplored in current mathematics [59][60][61] .
The use of the geometric mean of conductivity can be justified by the following consideration. When solving the fluid flow between adjacent blocks of porous material by the Darcy-Buckingham law (equation (1)), it is necessary to determine the conductivity between blocks at different saturation. The first idea could be to use the most common arithmetic mean. However, the arithmetic mean is applicable only for the extreme case of parallel pores 11 . In case of the other extreme -pores in a series -Hunt et al. 11 demonstrated that the harmonic mean should be used. In numerical experiments with random pore networks, Jang et al. 47 concluded that the geometric mean of both values is the most appropriate.

Conclusion
A semi-continuum model for the description of unsaturated homogeneous porous media flow is presented. The model is based on the idea of Macro Modified Invasion Percolation, in which the porous material is divided into blocks of non-infinitesimal size. Each block is represented by its retention curve and relative permeability. Saturation and pressure are considered continuous but constant within each block. Flow between adjacent blocks is described by the Darcy-Buckingham law. The limit of such a semi-continuum model is very similar to the standard Richards' equation. However, it differs in the way the pressure saturation relationship is captured. The retention curve has to be scaled appropriately to the size of the block. This results in a Prandtl-type hysteresis operator appearing under the derivative in the limiting equation. This limit differs from the standard RE, which is not able to describe finger-like flow. However, the physics behind both RE and the semi-continuum model is almost the same. Thus, the limit introduced above can be viewed as a reformulation of the RE in such a way that it does not lose the ability to describe finger-like flow. We conclude that the RE should be reconsidered by means of appropriate modelling of the hysteresis and proper scaling of the retention curve.
The limit of the semi-continuum model defined by equations (19) and (20) represents a rather interesting mathematical object. We are not aware of any research dealing with equations of such type. Since such equations seem to arise naturally by a limiting process of the semi-continuum model, we think they deserve more attention of the mathematical community.

Data Availability
No experimental data were generated or analysed during the current study.

Code Availability
The code that produced the simulations is available in MatLab upon request from the corresponding author.