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.

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 modelling, such as oil recovery[1] and soil physics[2,3]. The infiltration of rainwater into the soil and the recharge of groundwater reserved by soil water form the basis of the hydrological cycle of land. In this paper, we present a limit of the semi-continuum model[4] that describes the infiltration of water into a porous material such as soil.

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[5,6] which, in turn, are directly affected by the geometry of the porous matrix[7].

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[8]. The retention curve gives the matrix potential (or capillary pressure) as a function of saturation under equilibrium conditions[9,10]. It was shown that the retention curve and other physical properties depend on the volume of the sample[11-15]. However, this important issue has not usually been considered in porous media flow modelling.

Three different approaches to porous medium flow modelling 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)[16]. REV is the smallest volume for which the saturation and pressure fields (and other characteristics) can still be considered smooth, so that these characteristics provide a representative value of the whole. A representative value cannot be defined for smaller volumes than REV. 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[17]. All physical quantities (saturation, pressure, and time) are usually considered discrete in these models: The saturation is either one (full pore) or zero (empty pore), and the course of infiltration is followed in discrete time steps—n each step a single pore is filled or emptied. (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[4,18]. Saturation and capillary pressure are assumed continuous in time but constant within each block (i.e. piecewise constant in space), 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[19]. (B) At the level of capillary pores, various displacement rules have been proposed[20-22]. (C) At the level of large, non-capillary pores, the kinematic wave equation has been used[23].

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[17,24-30]. The oldest model is a combination (1A), which is expressed by the Richards’ Equation (RE)[31]. It has long been known that this model works well only for diffusion-like flow and fails to describe finger-like flow[32,33]. 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[34-44]. Another attempt is the combination (3A) proposed in[4,18,45,46]. It has been shown that some of these models describe diffusion-like flow, finger-like flow, and the transition between the two regimes[4,46].

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: How do they differ? This article seeks answers to this question. To make the argument as clear as possible, we present the one-dimensional form of the problem. The main characteristic of finger flow is the saturation overshoot[47,48]. Thus in 1D, we concentrate on the ability of the model to capture the saturation overshoot.

In this article, (1) we consider the limit of the semi-continuum model[4] as the block size goes to zero. Because the retention curve depends on the physical dimension of the block, (2) we introduce a physically relevant scaling of the retention curve so that the slope of the retention curve decreases with the sample size. (3) We introduce a novel type of porous media flow equation which is obtained as the limit of the semi-continuum model using appropriate scaling of the retention curve. (4) We analyse how this equation differs from the RE.

Methods

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[16]. The flux q (m/s) modeled by the Darcy–Buckingham law takes the form:where (m2) denotes the intrinsic permeability, (kg/m3) the fluid density, g (m/s2) 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). Capillary pressure in a liquid is created by the action of capillary menisci. The capillary meniscus is the curved surface forming the interface between the liquid phase (water) and gas (air) caused by surface tension. Pressure in the liquid phase is determined by the curvature of the menisci by means of the Young–Laplace equation[49]—tension in the liquid phase increases with the curvature of the menisci. The contact angle at the meeting points of all three phases (liquid, gas, and solid matrix) forms the boundary condition for the Young–Laplace equation.

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[16,50].

RE is then obtained by the combination of the Darcy–Buckingham law and the mass balance law equations[16] 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[16]. 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[18] proposed a cellular automaton under the title of “mechanistic modelling”, or “Macro Modified Invasion Percolation”.

In Kmec et al.[4,46], 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 (m2) and height L (m) filled with homogeneous and isotropic porous medium. We divide the tube into blocks of height . 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 [] with . The key quantities that we want to track are:

Saturation S (–) in each block.

Capillary pressure P (Pa) in each block.

Flux (m/s) between the blocks i and . Volumetric fluxes in [] can be recovered as . Naturally, only fluxes between pairs of neighboring blocks are nonzero in this setting.

At each instant, the pressure and saturation are considered constant within each block. We thus use the notations and , . Gravity is directed downward along the long axis of the tube, which is called the x-axis here. A constant (in time) influx (m/s) is set across the top boundary (at ). Zero discharge (i.e. zero flux) is assumed at the bottom boundary (at ). 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 . At each time step, the saturation in each of the blocks is updated according to Eq. (4), which is an explicit discretized version of the semi-continuum model [Eq. (3)] A backward time discretization of Eq. (3) can be also used (see the Discussion in[46]).

Next, the capillary pressure in each block is updated according to the retention curve. There are many approaches to modelling hysteresis[51-53], but here 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 [8,54,55]. 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[56] and is similar to the play-type hysteresis used e.g. in[57]. 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.

Figure 1

Model of the retention curve hysteresis. The coloured lines correspond to different parts of the retention curve.

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 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 i and , we simply take the geometric mean of the permeability of the respective blocks. The geometric mean 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[58].

The semi-continuum model[4] works as follows:

The size of the blocks is chosen and an appropriate time step 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.

The top boundary condition is set: The flux into the topmost block is set and fixed to . The bottom boundary condition is set: The flux out of the bottom block is set and fixed to zero.

Using the current value of the fluxes and , saturation in each block is updated according to Eq. (4).

Pressure in each block is updated according to retention curve and hysteresis model, keeping track whether the block is in the imbibition or draining mode.

Fluxes between neighboring blocks are updated according to Eq. (5), keeping the boundary fluxes fixed by step 2.

Time is updated to and the process goes back to step 3.

The fundamental difference between a numerical scheme for the RE and the semi-continuum model lies in the scaling of the retention curve. In a numerical scheme for RE, the retention curve remains the same as the block size goes to zero. In the semi-continuum model, the slope of the retention curve vanished with vanishing block size. Thus, with , the semi-continuum model does not converge to the RE.

Results

In Kmec et al.[4], the authors showed that the model is able to reproduce well all the observed features in unsaturated flow in a narrow vertical test tube filled with sand in experiments reported by DiCarlo[47,48,59]. The model correctly predicts when a saturation overshoot effect will appear. Moreover, it captures well both the interesting aspects of the overshoot behavior: (1) the non-monotonic dependence of the overshoot magnitude on the influx, and (2) the transition from the overshoot regime to diffusion-like regime for increasing initial saturation. See[4] for all the details.

Fingering regime in a narrow tube does not allow the fingers to exhibit complicated spatial patterns. To observe these, one has to switch to two-dimensional experiments such as[18,32,60-65]. It is straightforward to extend the semi-continuum model to two spatial dimensions. However, to keep the setting as clear as possible, we stick to the simplest 1D formulation in this text. In[46], the authors used the two-dimensional version of the model to show that it is able to correctly reproduce the transition between a fingering regime with saturation overshoot (for small initial saturation) and a diffusion-like regime of a stable flat water front with a monotonic saturation profile (for large initial saturation)—see Fig. 6 in[46].

We can clearly see that the simulations of the semi-continuum model are in good agreement with the experimental results. However, the crucial question how to choose the block size is still missing. This is clearly a “parameter” of the semi-continuum model, and clearly a rather artificial one.

First, note that the time step is not a free parameter of the model—it is a discretization parameter. Figure 2 (left) shows the behavior of the 1D semi-continuum model for a range of values. The solution is stable in the limit . This is not the case for . If we let the block size go to zero ( has to go to zero, too), and kept the retention curve constant, the overshoot behavior would disappear (see Fig. 2 right).

Figure 2

Left panel: The dependence of the moisture profile at mins on time step , the profiles are almost identical for a range of between and . Right panel: The dependence of the moisture profile at minutes on . As , without proper scaling of the retention curve, the overshoot behavior would disappear. The parameters used for the simulations are given in Table 1.

Table 1

Parameters used for the simulations presented in Fig. 4.

Parameter	Symbol	Value
Porosity	θ	0.35
Density of water	ρ	1000 kg/m3
Dynamic viscosity of water	μ	9×10^-4Pas
Intrinsic permeability	κ	1×10^-10m^2
Relative permeability	k(S)	S^3
Acceleration due to gravity	g	9.81m/s^2
Large gradient	K_PS	10^5

In the limiting process mentioned above (Fig. 2), the dependence of the retention curve on the block size[11,13] is missing. 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 is implemented. This procedure is named scaling of the retention curve and will be discussed in the next section.

Scaling of the retention curve

The scaling of the retention curve 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 changes. The fluxes are given by Eq. (5) in which decreasing 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 formfor the main wetting branch, andfor the main draining branch, where (Pa) and (Pa) are constants. The parameter h (m) is the scaling parameter equals to the block size . 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[50].

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 Eqs. (6) and (7) with Pa and Pa roughly matches the main branches of retention curve of 20/30 sand in the experiments of DiCarlo[47]. Figure 3 illustrates the scaling of this retention curve as goes to zero. As decreases, the slope of both branches of the retention curve decreases, too. In the limit, both branches take the form of horizontal lines.

Figure 3

Scaling of the retention curve with the block size . Grid block sizes are colour-coded. The solid line denotes the main wetting branch and the dashed line denotes the main draining branch.

In practice, the retention curve will be measured for a sample of a known dimension . The scaling is then given by Eqs. (6) and (7) and there is no need to measure the retention curve repeatedly for blocks of different dimensions.

Convergence of the moisture profiles

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 top boundary influx is shown in Fig. 4. The parameters used for the simulations are given in Table 1. Since we are not interested in the bottom part of the porous medium, this part is not depicted here. Left panel of Fig. 4 shows the calculated moisture profiles for the initial saturation for a decreasing sequence of block size values. Right panel of Fig. 4 shows the numerical convergence for in case of a diffusion-like flow regime (the initial saturation ). Moreover, the moisture profile convergence is not affected by the different top boundary flux. Figure 5 shows numerical convergence for in the case of two different top boundary fluxes varying in two orders of magnitude. A different scale of the y-axis is used in the right panel of Fig. 5 to show the details.

Figure 4

Left panel: Convergence of the moisture profile at minutes for for initial saturation , and constant top boundary flux m/s. The flux out of the bottom block is set to zero. Each colour corresponds to a particular grid block size. The moisture profile converges and retains the overshoot pattern. Right panel: Convergence of the moisture profile at minutes for for initial saturation , and constant top boundary flux m/s. The flux out of the bottom block is set to zero. Each colour corresponds to a particular grid block size. The moisture profile converges to a sharp water-front without saturation overshoot.

Figure 5

Left panel: Convergence of the moisture profile at minutes for for initial saturation , and constant top boundary flux m/s. The flux out of the bottom block is set to zero. Each colour corresponds to a particular grid block size. The moisture profile converges and retains the overshoot pattern. Right panel: Convergence of the moisture profile at minutes for for initial saturation , and constant top boundary flux m/s. The flux out of the bottom block is set to zero. Each colour corresponds to a particular grid block size. A smaller y-axis scale is used to show the details of the moisture profiles. The moisture profile converges and retains the overshoot pattern.

We may observe that this retention curve scaling preserves the character of the flow across all levels of 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.

The scaling of the retention curve described above can be also demonstrated in the 2D version of the semi-continuum model presented in[46]. Figure 6 shows the predicted moisture profile for 2D simulation for a decreasing sequence of block size values. Notice analogous behaviour to the 1D case: The saturation in the fingertip increases slightly with decreasing block size, but the oversaturated zone remains roughly the same. However, due to numerical errors that are not present in 1D, the fingers become more narrow and consequently slightly faster. The limiting process in 2D requires more attention and will be addressed in a subsequent paper. The numerical evidence of the 1D and 2D cases 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 , if the retention curve is scaled in this appropriate way.

Figure 6

Convergence of the moisture profile in 2D at time minutes for cm, 0.250 cm, 0.125 cm and 0.0625 cm from the left to the right for initial saturation . The moisture profile converges and retains the overshoot pattern. Saturation values are colour-coded according to the colour bar on the right.

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 to a function and show which equation this function should satisfy. This equation is called the limit of the semi-continuum model. A detailed derivation of this limit is given in the Supplementary Information. Equations (8) and (9) represent the classical form of the limit of the semi-continuum model. It is a partial differential equation containing a Prandtl-type hysteresis operator [defined by Eq. (8)] 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. Let us also note that if the saturation is continuous, then .

Similar notation as in equation (9) is also used for the classical RE [see Eq. (10)]. The two equations differ in how the relationship between pressure and saturation is expressed. In the limit of the semi-continuum model [Eq. (9)], the pressure-saturation relation is defined by a Prandtl-type hysteresis operator , while for the RE [Eq. (10)], the pressure-saturation relation is described by a hysteretic operator P with wetting and draining branches represented by monotonically increasing functions. Moreover, Eq. (9) preserves the geometric mean of the relative permeability.

In light of the analysis above, we argue that the RE arises from a physically unsound limit of the semi-continuum model which ignores the proper scaling of the retention curve. Thus, it is not surprising that the RE is unable to capture overshoot behavior. If the proper scaling of the retention curve is included, we obtain a different limit and the overshoot behavior is not lost.

Discussion

The limit of the semi-continuum model [Eqs. (3), (4), (5)] was found by means of mathematical analysis (see Supplementary Information) in the form of a partial differential Eq. (9). 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 Eq. (5) it follows that the shape of the retention curve must be scaled according to the block size [Eqs. (6), (7)].

The process of scaling the retention curve is not a common practice in flow modelling. 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[66]. The effects of sample dimension on capillary pressure have since then been pointed out e.g. in[11-13,67-69]. 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[70,71], 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[70,71]. 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[13]. Pražák et al.[70] 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[14].

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[49], the pressure drop across both menisci is where is the surface tension of the liquid. Setting the gas pressure to zero, the liquid drop is under tension (i.e., negative pressure) . Connecting such an empty pore to a liquid reservoir at a pressure lower than 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 , 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 . Let us continue this thought experiment and consider two pores of radii . At a certain pressure , the first pore will fill in completely, and at a higher pressure , the second one will fill in. Thus, the retention curve of this pair of two pores becomes a broken horizontal line (see Fig. 7). 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.

Figure 7

The main draining branch of two pores of radii and the corresponding pressures and . Both pores are assumed to have the same volume, thus the horizontal line is broken at .

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 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 limit of the semi-continuum model 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[72]. However, in the case of the semi-continuum limit, the hysteresis operator appears inside the Laplace operator. Such equations seem to be unexplored in current mathematics[73-75].

The idea of taking REV into account in finite element discretization is completely unusual (at least in the mathematical community), surprisingly, it has long been used in porous media[76-79]. For instance, in[77], the authors proposed a method for the calculation of the effective permeability of sandstone with the lattice Boltzmann finite element method. The method consists of two models—pore-scale and macro-scale. First, the permeability tensor is calculated by the lattice Boltzmann method (pore-scale model). This permeability is then used for solving the diffusion equation (macro-scale model) on a finite element mesh. A crucial step in the model is determining an appropriate finite element discretization. The authors argue that a highly refined mesh cannot be used, because the size of the elements will then be so small that the continuum approximation is no longer appropriate. Thus, the size of the REV is estimated and used as a lower limit for the size of the finite elements. This means that the computational mesh takes into account the dependence of physical parameters on the size of mesh elements.

It has long been known that macro-scale models (i.e., models larger than REV) often fail to describe experimentally observed phenomena[80,81]. This usually happens because some physical properties (e.g. concentration) vary significantly within a REV. This may be an essential information for simulations[82]. Obviously, this cannot be captured by continuum-scale equations, which cease to be sufficiently accurate. A pore-scale simulations such as lattice Boltzmann methods are used in this case, however, there are two major drawbacks that make the model impractical as a predictive tool. (1) The simulations are computationally demanding and (2) lack detailed information about the pore geometry in porous media. Therefore, a logical solution is to use a pore-scale approach where the detailed structural information of the porous medium is necessary, and to use macro-scale approaches where such detailed information are not needed. This approach to modelling is represented, e.g., by multiscale[83-86] and hybrid algorithms[80-82,87,88]. Multiscale algorithms use pore-scale models to evaluate physical properties such as permeability, density, and/or velocity. These are then upscaled into a macro-scale model. The model presented in the previous paragraph[77] is an example of a multiscale algorithm. On the other hand, hybrid algorithms use a slightly different approach: In regions where continuum methods fail, they are replaced by the pore-scale methods. Thus, the porous medium is divided into two disjoint domains, in which different algorithms are used. Sometimes it is even needed to use several different methods across several scales[86]. For both multiscale and hybrid algorithms, the methods are strongly dependent on the size of the porous medium. There are similarities between our semi-continuum approach and hybrid/multiscale modelling in porous media. We use a retention curve that scales with the block size; thus we use a different formula dependent on the size of the porous medium. It is possible to implement an adaptive block size; in the region where the finger-like flow occurs, small blocks can be used, and in the region where the diffusion occurs, larger blocks can be used. This is a formal analogy to hybrid modelling.

At the pore scale, the fluid motion is characterized by capillary displacement and the Darcy–Buckingham law is not longer appropriate[16,20]. The transition between two different models of water displacement is captured by the semi-continuum approach. In the limit, the main wetting and draining branches degenerate into parallel lines which results in the Darcy’s law taking the character of capillary displacement. Indeed, water fills and drains the pores almost instantly, similarly as in percolation theory. Therefore, we actually use “different” physics for Darcy and sub-Darcy scales. This is achieved by varying the slope of the retention curve, not by explicit changes of the governing equations.

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 [Eq. (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[13]. In case of other extreme—pores in a series—Hunt et al.[13] demonstrated that the harmonic mean should be used. In numerical experiments with random pore networks, Jang et al.[58] concluded that the geometric mean of both values is the most appropriate.

A comparison of the RE (10) and the limit of the semi-continuum model (9) shows that there are two major differences between them. In the semi-continuum model (1) a proper scaling of the retention curve is used, and (2) the geometric mean of the conductivity between adjacent blocks is maintained. Note that the hysteresis has to be included in the retention curve because it is known[64] that the saturation overshoot is a consequence of the retention curve hysteresis. With (1) and (2), any reasonable numerical solution of the RE will converge to the limit of the semi-continuum model. In this case, it will correctly simulate the diffusion-like flow, the finger-like flow, and the transitions between them, just as the semi-continuum model can.

Conclusion

A semi-continuum model for the description of unsaturated homogeneous porous media flow is presented. One can see the similarity between the proposed model and the multiscale/hybrid modelling. 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 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. Moreover, the geometric mean of adjacent blocks is maintained. 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 Eqs. (8) and (9) represents a rather interesting mathematical object. From a mathematical point of view, the Richards’ equation is a parabolic Eq. (10), but the limit of the semi-continuum model switches between parabolic and hyperbolic type. 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 soil science and mathematical community.

Supplementary Information.