Transient Two-Phase Flow in Porous Media: A Literature Review and Engineering Application in Geotechnics

Abstract

This work reviews the transient two-phase flow in porous media with engineering applications in Geotechnics. It initially overviews constitutive relationships, conventional theories, and experiments. Then, corresponding limitations are discussed according to conflicting observations and multiphase interfacial dynamics. Based on those findings, the dynamic nonequilibrium effects were so defined, which could also be abbreviated as dynamic/transient effects. Four advanced theories have already been developed to resolve these effects. This review collects them and discusses their pros and cons. In addition, this work further reviews the state-of-art in terms of experimental methods, influential factors in dynamic/transient effects, and modelling performance, as well as micromodel and numerical methods at pore-scale. Last, the corresponding geotechnical applications are reviewed, discussing their applicability in effective stress, shear strength, and deformation. Finally, the entire review is briefed to identify research gaps in Geotechnics.

1. Introduction

Multiphase flow in porous media is a complex engineering problem. It covers various disciplines, including agriculture, hydrogeology, geotechnical engineering, petroleum engineering, biological engineering, and medical engineering. All of these disciplines require a sound understanding of multiphase motion in porous media, such as air-water moving into the soil, gas-oil-brine moving in the fractured rock, blood flowing into tissues, solution flowing into porous synthetics, and ink being dipped on papers. Every phenomenon of fluid anti-gravitational motion in porous media is dominated by the capillary effect on the interface between two immiscible fluids. As for the geotechnical engineering perspective, most concerns are mainly given in two cases. The first is to study air-water motion in unsaturated soil. In this case, the air is the nonwetting phase, and water is the wetting phase due to the hydrophilic character of the natural ground. Thus, the air-water interface governs the two-phase flow in the soil matrix and further impacts soil strength, stiffness, and deformability [1,2]. The second instance is a gas, oil, and brine system in the oil/gas reservoir. A large amount of crude and natural gas are stored in natural geological formations under overpressure. With the ongoing oil production, overpressure decrease leads to a reduction of oil production. For enhancing oil production, petroleum engineers often utilise gas or water injection [3]. The theory of multiphase flow in porous media offers geotechnical engineers solutions to determine the hydro-mechanical behaviour of unsaturated soil under various environmental and hydraulic conditions [4]. It also can predict oil recovery efficiency during the gas/water injection process in deep reservoirs [3].

The first pioneer studying multiphase flow in porous media is the famous soil physicist Edgar Buckingham, who invented the Buckingham π-theory in dimensional analysis. Narasimhan [5] reviewed their work on laying the foundation of the soil moisture retention behaviour and capillary flow conductivity formularization. Although their work did not involve scientific measurement of soil suction due to the constraints of the measuring method of that period, the findings from their experimental results and scientific intuition have significantly influenced the orientation of the multiphase flow study. In general, during the period between the end of the nineteenth century and the beginning of the twentieth century, inspired by studies in heat transfer (Fourier’s Law) and electrical conduction, soil physicists commenced constructing seepage theory for saturated soil and partially saturated soil based on the diffusion theory [5,6]. This theory was also applied to study solute transportation in surface water and groundwater seepage [7]. Therefore, the state variables and constant parameters selected for groundwater modelling are somehow similar with the diffusion theory’s state variables and constants in both heat and electricity transfers, such as the hydraulic potential to the electrical potential, the hydraulic conductivity to the electrical conductivity (reverse of electrical resistivity), and the water storage capacity (specific storage) to the electrical capacity, etc. [5].

Since the single-phase and multiphase flow theories are based on diffusion, two challenges were left to engineers and applied mathematicians. The first one is the diffusion coefficient determination, which was later decomposed by the chain rule of calculus into different constants, including permeability, specific storage, and capillary specific storage [7]. As a result, diligent engineering experiments are required to reveal the nature of the assumed constants. Under the single-phase flow seepage theory, this diffusion coefficient (hydraulic diffusivity), which is hydraulic conductivity divided by specific storage, can be treated as a constant with incompressible fluid and rigid soil matrix assumptions [7]. However, as for the multiphase flow seepage in porous media, this coefficient is a non-constant parameter, depending on phase potential or saturation [1,2,6,7,8]. The same circumstance also happens for the single-phase flow if the porous media deformation or fluid compressibility were counted [2,7,9]. This formulation is usually defined as the hydro-mechanical coupling model for fully saturated porous media. The second challenge is to derive the solution of diffusion partial differential equation (PDE) under the circumstance that the diffusion coefficient is a constant or a parameter continually changing with the state of phase. With the increasing complexity of the constitutive relationship to determine the diffusion coefficient, the original diffusion PDE becomes nonlinear [6,7,10]. The analytical solution can only be derived by simplifying fluid phase potential-saturation-permeability constitutive relationships into the exponential form [11]. Provided complex and flexible constitutive relationships, the nonlinear PDE can only be solved using numerical methods [2,12]. In summary, nonlinear PDE has been successfully solved by computational mechanical engineers and applied mathematicians [4,9,11,12]. Moreover, different solutions have been given and are reasonable under the state-of-art of various constitutive relationship formulations [13,14,15,16,17,18,19,20].

With the development of the sensor technique, Topp, et al. [21] experimentally discovered the dynamic effects in the constitutive relationship of the nonlinear diffusion theory under transient flow conditions. Such observations cannot be simulated by solving the conventional theory both analytically and numerically. Hassanizadeh, et al. [22] reviewed some early validations of theory to find that the diffusivity depends on the speed of the wetting front. In addition, some of the literature reviews from the earlier research concluded that the diffusion equation could not be verified for moisture transport, and the application of Darcy’s Law is questionable [23,24]. Moreover, the relationship between diffusivity and moisture content loses uniqueness for transient flow conditions [22]. Since the first observation of dynamic effect by Topp, Klute and Peters [21], many soil column tests were carried out to validate the conventional theory [23,24,25,26,27]. Meanwhile, more attention has been paid to investigating the dynamic effects in soil suction or capillary pressure for each saturation or effective permeability [27,28,29,30,31,32], and its dependence on testing scales [33,34], intrinsic properties of porous media [35,36,37], fluid properties [38,39,40], and other impact factors [41,42], etc. In order to simulate the transient two-phase flow in porous media under nonequilibrium conditions, several works developed a series of novel theories from various physical bases [22,43,44,45]. By far, the investigations of dynamic effects have been still continually conducted using experimental [46], theoretical [47], and numerical methods [48] in the continuum scale. Diamantopoulos and Durner [49] comprehensively reviewed the so-defined dynamic nonequilibrium effects from a soil science and hydrology perspective, merely focusing on the air-water system in natural soil. Another most recent review can also be sourced from Li, et al. [50] in petroleum engineering, more focusing on the oil-water/gas system in deep tight reservoirs. Beyond the investigations on the continuum scale, their work also briefly reviewed microscale physical pore network models [51], numerical pore network models [38] and Computational Fluid Dynamic (CFD) simulations [52,53,54,55,56]. In summary, all reviews draw similar conclusions on the non-negligibility of dynamic nonequilibrium effects under transient flow conditions, the dependence of dynamic effects on fluid and porous media intrinsic properties, the importance of further experimental development, theoretical expansion, determination of constitutive parameters using inversion analysis, microscale physical and numerical modelling, etc.

Previous reviews of dynamic nonequilibrium effects in soil water flow from Hassanizadeh, Celia and Dahle [22] and Diamantopoulos and Durner [49] were sound and adequate in soil science and hydrology. However, those two works were published almost ten and twenty years ago and, therefore, could not involve recent contributions to this research objective. Moreover, the latest review from Li, Luo, Li, Liu, Tan, Chen and Cai [50] was given briefly without expanding details and merely focused on the dynamic effects in the tight reservoir with ultra-low permeability. Thus, there is a lack of a comprehensive overview of the state-of-art by far, especially from Geotechnics. Several decades ago, Geotechnics expanded saturated soil mechanics to include unsaturated soil in terms of unsaturated soil water seepage and retention behaviour by introducing the unsaturated soil studies in Agriculture [2]. Later, the unsaturated soil hydrology enriched the geomechanics in terms of unsaturated soil effective stress, shear strength [57,58,59,60], and deformation [61,62,63,64,65]. Although the foundations of unsaturated soil mechanics were solidly laid above the unsaturated soil hydrology, new findings from soil hydrology, such as the dynamic effects for transient seepage, were rarely introduced into Geotechnics. Geotechnics embraces the hydro-mechanical behaviour of a variety of geological materials, e.g., unsaturated soil and deep reservoir. Therefore, advancing the transient seepage flow from instantaneous equilibrium to nonequilibrium conditions is essential by introducing advancements from soil hydrology. Moreover, more dynamic effects were recently detected in natural slopes by Bordoni, et al. [66] and low permeable reservoirs from Tian, et al. [67]. Hence, it is convincing that the dynamic nonequilibrium effects, usually neglected in Geotechnics, will facilitate a better understanding of the unsaturated soil suction and shear strength reduction. It will also contribute to a more accurate temporal prediction of the geohazards (e.g., dam, shallow slope failure, etc.) triggered by extremely intensive rainfall and flooding. In addition, the more precise estimation of deep reservoir production by considering dynamic effects can be applied to conduct appropriate reservoir development plans in reservoir Geotechnics [50].

Based on the considerations mentioned above, this work aims to deliver a comprehensive and state-of-art review of transient two-phase flow in porous media under both instantaneous equilibrium and nonequilibrium conditions for any potential applications in Geotechnics. This review starts with an illustration of conventional two-phase seepage from fundamental multiphase physics to traditional theories, appending with standard testing methods. Then, the constraints and criticism of those theories and experiments are summarised as a lead-in of the dynamic nonequilibrium effects. Subsequentially, this work overviews the state-of-art of nonequilibrium transient two-phase seepage. It mainly covers advanced theories for modelling with novel experiments in the continuum scale and microscale numerical modelling. Lastly, a few studies for applications in Geotechnics are briefly reviewed. This section could orientate the future development of geomechanics with the transient two-phase flow in porous media. Finally, a summary and conclusion enclose this review as referencing principles for future experimental, numerical, and theoretical explorations and applications in Geotechnics.

2. Soil Suction and Soil Water Retention Behaviour

2.1. Multiphase Physics and Definition of Soil Suction

The total suction is defined as the free energy state of soil water. A thermodynamic-based expression of total suction could be given by Kelvin’s Equation, which provides a relation between this free energy of the soil water and the partial pressure of the pore water vapour [2]. According to Fredlund and Rahardjo [2], the total suction consists of matric and osmotic suction. The matric suction can be further decomposed into suction induced by the capillary effect and short-range adsorption effects, including electrical and van der Waals force fields [1]. Lu and Likos [1] also stated that short-range adsorption effects are most significant for fine soil with a large surface area and most related to low water content conditions. As for the non-cohesive granular soil, the capillary effect, shown in Figure 1a, will be taken as the only one that contributes to matric suction. The Equation for soil suction can be given as

$${\psi }_t=\frac{RT}{v_w}\ln (RH)={\psi }_m+{\psi }_o=(u_a-u_w)+{\psi }_o=\frac{2{\sigma }_s\cos \mathsf{\theta }}{r}+CRT(1+B_2{C}^2+B_3{C}^3+\dots )$$

(1)

where ψ_t = total soil suction (kPa), R = universal constant [8.31432 J/(mol K)], T = absolute temperature (K), v_w = partial molar volume of liquid water (1.8 cm³/mol at 20 °C), RH = relative humidity [RH = u_v/u_v0, u_v = partial pressure of vapour (kPa), u_v0 = pressure of water vapour over a flat surface of pure water at the same temperature (kPa)], ψ_m = matric suction (kPa), ψ_o = osmotic suction (kPa), u_a = air pressure (kPa), u_w = water pressure (kPa), σ_s = surface tension of soil water (N/m), θ = contact angle of water-air interface, r = pore radius (m), C = molar concentration of solute in the pore solution (mol/m³), and B_i = viral coefficients [68].

2.2. Soil Water Retention Curve

In unsaturated soil mechanics, the pore matric conceptual model is usually considered a bundle of capillary tubes of tube radius, statistically distributed in a soil Representative Elementary Volume (REV). The soil suction estimated from this REV scale is an upscale matric suction by spatial averaging. The Soil Water Retention Curve (SWRC) concept is the relationship between this soil matric suction and water saturation/moisture content for a single soil REV. An S-shape curve can be fitted into these data points after measuring the moisture content at different equilibrium soil suction. Figure 2 shows a comprehensive sketch of SWRC, including drainage, imbibition, and hysteresis scanning loops.

In the domain of static SWRC (Figure 2), there are mainly two boundary curves: primary drainage and imbibition curves, which constrain all possible SWRC under the varying history of hydraulic loading. Due to the ink-bottle effect inside the pore matrix (Figure 1b), according to Equation (1), capillary water prefers to pass through smaller pores rather than large pores, which requires more water pressure to penetrate. Every SWRC, except primary drainage curves, own this nature, so the hysteresis scanning loops exist between two boundaries curves. There is a small amount of irreducible film water near the high suction range because of the short-range adsorption effect. In this condition, soil water exists in molecule film strongly attracted by mineral particle surface charge and van der Waal forces. Gravity or high gas pressure cannot drain them out unless the high-temperature oven is applied. The water content (or saturation) in this suction range is residual moisture content (or residual saturation). Air entry value is a suction data point where full saturation starts to decrease. This concept is usually adopted when the capillary flow is modelled using drainage SWRC. Based on the conceptual model of the pore matrix, this point can be identified as a suction threshold where the meniscus initiates in the largest pores. In fact, a single drainage curve cannot fully recover SWRC of natural soil due to the previously mentioned ink-bottle effect and gas trapping mechanism. However, many numerical modellers still prefer a single drainage curve rather than including hysteresis.

2.3. Soil Water Retention Function

Soil Water Retention Function (SWRF) is a continuous function fitted into experimental data. This function allows continuously extracting the relation between water content and soil suction for numerical modelling. In

Table 1. Soil Water Retention Function (SWRF).

SWRF Authors	Fitting Functions
Gardner [71]	$S_e=\frac{1}{1+{\alpha }_G{\psi }^{n_G}}$	(2)
Brooks [13]	$\begin{array}{l}{S}_e=\{\begin{array}{c}\begin{array}{cccc}1& & & \psi <{\psi }_{AEV}\end{array}\\ \begin{array}{cccc}{\left(\frac{{\alpha }_{BC}}{\psi }\right)}^{n_{BC}}& & \psi >{\psi }_{AEV}& \end{array}\end{array}\\ {\alpha }_{BC}={\psi }_{AEV}\end{array}$	(3)
Van Genuchten [14]	$S_e={\left[\frac{1}{1+{\alpha }_{VG}{\psi }^{n_{VG}}}\right]}^{m_{VG}}$	(4)
Fredlund and Xing [15]	$S_e=\left[1-\frac{\ln (1+\frac{\psi }{{\psi }_r})}{\ln (1+\frac{{10}^6}{{\psi }_r})}\right]\frac{1}{{\left\{\ln \left[e+{(\frac{\psi }{{\alpha }_{FX}})}^{n_{FX}}\right]\right\}}^{m_{FX}}}$	(5)

where ψ = the soil suction, $S_e=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}$ is the effective saturation; ${\theta }_s,{\theta }_r$ are the saturated and residual water contents; $\theta ,\theta \in [{\theta }_s,{\theta }_r]$ is water content; ${\psi }_{AEV}$ is air entry value; ${\alpha }_i,n_i,m_i(i=G,BC,VG,FX)$ are fitting parameters for each model; ${\psi }_r$ is the soil suction for irreducible water content.

, there are mainly four SWRF applied for nondeformable soil in geotechnical engineering and soil hydrology [1,7,15]. The fitting parameters determine the shape and air entry value of SWRF.

Soil water retention function initially appeared as an empirical function fitted into data. One of the physical interpretations can be sourced from Fredlund and Xing [15]. In this work, Fredlund and Xing [15] derived each previous SWRF by assuming different pore size distribution functions (PSD). This connection is based on both the Young–Laplace Equation and the statistical distribution of the pore matrix. The derivation of SWRF from PSD is given by

$$\theta (R)={\displaystyle {\int }_{R_{\min }}^{R}f(r})dr=={\displaystyle {\int }_{\psi }^{{\psi }_{\max }}f(\frac{2{\sigma }_s\cos \mathsf{\theta }}{\delta }})\frac{2{\sigma }_s\cos \mathsf{\theta }}{{\delta }^2}d\delta =\theta (\psi )$$

(6)

where θ = volumetric water content; f(r) = the function of pore size distribution; R = the pore radius (m); R_min = the smallest pore radius; ψ = the soil suction (kPa); ψ_max = the maximum soil suction (kPa); σ_s = surface tension of soil water (N/m); θ = contact angle of water-air interface; δ = dummy variable of soil matric suction [15]. This derivation opened a gate for studying SWRC of deformable soil, where the variation of PSD can be measured using Mercury Intrusive Porosimetry (MIP) [65].

For nondeformable soil, SWRC can be fitted by these four models. Leong and Rahardjo [72] reviewed these four equations against many soil water characteristic data from previous publications and eventually suggested that the Fredlund and Xing [15] model provides the best fit among all. Still, the van Genuchten [14] model and the Fredlund and Xing [15] model without a suction correction factor have better-fitting performance for sandy soil. Later, many parameter studies agreed with this recommendation: both models provide similar fitting curves except for the high suction range. Moreover, they did not recommend using the correction factor in the Fredlund and Xing [15] model because it lacks physical support for absolute zero water content at the suction of 10⁶ kPa.

2.4. Problems in Soil Water Retention Behaviour

In fact, several factors influence the uniqueness of SWRC. Zhou and Yu [73] reviewed three effects impacting the uniqueness of SWRC: the initial void ratio, initial water content, stress state, and high suction values. Malaya and Sreedeep [74] reviewed another impact in the suction measurement methodology, stress history, additives, ageing, and measuring range of suction. Regarding these two review studies, the effects impacting the uniqueness of SWRC are mainly governed by two characteristics. The first one is geometry characteristic of soil (e.g., soil structure, pore structure), influenced by sample fabrication, stress history, and wetting-drying history. The other one is chemical characteristic, which is affected by soil mineralogy, solute concentration, and wetting phase temperature variation [75,76]. Jotisankasa [69] concluded that SWRC is firstly dominated by soil matric at the low suction range, and the other characteristics have more influence on SWRC at the high range. Thus, to study one factor impacting the non-uniqueness of SWRC, other conditions have to be maintained through the entire experimental operation.

Through following a review from Malaya and Sreedeep [74], several impact factors are worth to be further studied:

• The SWRC variation of deformable soil because of hydro-mechanical loading;
• Hysteresis of SWRC additionally involving stiffness variation of deformable soil caused by hysteresis and densification of collapsing soil induced by hysteresis;
• The time-dependent change of SWRC because of wetting phase reconfiguration at transient state.

The first research question draws more attention from geotechnical researchers because it highly concerns the determination of stress state variables, stress-strain behaviour, and the hardening behaviour of unsaturated plastic soil. Many studies were given to study SWRC variation under changing the void ratio and soil plasticity [65,77,78,79,80,81,82]; all of them have one thing in common, which is to construct a three-dimensional space of suction-saturation-void ratio (or specific volume). In this way, the usual SWRC measured using conventional pressure cells is just a projection of the state surface of the suction-saturation-void ratio in the two-dimensional suction-saturation plane. With the addition of the void ratio axis, the variation of SWRC caused by both mechanical and hydraulic loading can be tracked from the Soil Water Retention Surface (SWRS), rather than a single curve having a blurry specification of soil deformation history. Gallipoli, Wheeler and Karstunen [80] embedded specific volume change into the air entry value of the van Genuchten [14] model to capture this behaviour. The other method is to derive SWRS according to PSD variation [65]. Comparing both of them, it seems that the SWRS derived from PSD is more physically reasonable because PSD can be directly measured from the deformed specimens. Simply using varying air entry values to account for soil deformation is not able to reveal PSD variation. However, corresponding experimental effort simultaneously increases with the theory expansion (many MIP tests instead of standard tests) and, therefore, limits the application of new research findings.

The second research question is a concern for both geotechnical researchers and hydrogeologists because hysteresis behaviour changes the stress state variables of unsaturated soil and the specific capillary storage. Fredlund, et al. [83] studied the shift between a drying curve and a wetting curve and finally indicated that a median curve should be fitted between two curves for engineering application. They confirmed the shortage of this method but also stated the variability of in situ SWRC, so they encouraged further research on this topic. As for the hysteresis SWRC model of nondeformable soil, Pham, et al. [84] reviewed many physical-based and empirical models. Through comparing these models with 34 soil datasets, they concluded that Mualem [17] appears to be the most accurate and straightforward model for scanning loops prediction in engineering practice [84]. A novel model for SWRC hysteresis was recently developed by Pedroso and Williams [19] at the University of Queensland. This model can be calibrated using genetic algorithms [20] and has been successfully applied to a numerical simulation of unsaturated soil stress-strain with SWRC hysteresis [9]. Nevertheless, these hysteresis models neglected the volume change effect, which might be more significant than the hysteresis effect. Moreover, these hysteresis models were usually developed solely to model multiphase seepage, assuming that rigid soil matric does not affect capillary storage. For solving the unsaturated soil stress-strain problem in the geotechnical engineering domain, soil deformation and hysteresis simultaneously happen, which is considered less in the numerical simulation because the 3D hysteresis SWRS and corresponding experimental procedure are still in the stage of development. Three available hysteresis models accounting for void ratio variation are given in Gallipoli [81], Hu, Chen, Liu and Zhou [65], and Tsiampousi, et al. [85]. In their works, hysteresis curves are constructed between primary drying and wetting surfaces. Comparing the three studies mentioned above, the prediction of SWRS from PSD variation owns more physical meaning because it captures the air-entry value changing with specific volume and the SWRC gradient changing with the void ratio.

The last research question is one of the most interesting, which is more explored by soil scientists, petroleum engineers, and hydrogeologists but less concerned with geotechnical engineering. Time dependence of SWRC can also be defined as SWRC for transient multiphase flow or SWRC for dynamic (unsteady state) multiphase flow. SWRC is defined as a constitutive relationship that intrinsically exists in the unsaturated soil. If there is a non-negligible discrepancy between dynamic and static SWRC, the entire unsaturated soil mechanics theory might still miss this critical piece. Moreover, there should be no doubt that this effect could be coupled with soil shear strength, deformation, and hysteresis effects, if possible.

Despite these three popular research questions, there are still other studies working on simplifying the procedure of SWRC determination. Due to the nature of SWRC being strongly related to void structure, SWRC can be estimated from Grain Size Distribution (GSD) with appending parameters, such as soil density, fine content, organic content, and Atterberg limits. According to the literature review from Fredlund, et al. [86], there are two types of Pedotransfer Functions (PTF) to transfer GSD to SWRC. The first one is an entirely empirical approach, seeking the correlation between SWRF fitting parameters and grain size index or grain size fraction, combined with a density or plastic index [87,88,89,90,91]. The other one is the physical-empirical-model approach, for which estimation of SWRC is given based on the geometric assumption of granular particle, pore or constriction size distribution, and Young–Laplace equation [86,92,93,94]. The SWRC estimated by these methods cannot be perfectly fitted into experimental data, but the physical-empirical models show a better prediction for poorly graded uniform sand [86].

SWRF was initially designed for multiphase flow seepage simulation, in which the Darcy seepage equation is often introduced to express momentum conservation. In the next section, the governing equation derivation of multiphase flow in porous media will be briefly illustrated with the Hydraulic Conductivity Function (HCF) of the unsaturated soil and the association between SWRF and HCF.

3. Steady-State and Transient Two-Phase Flow Seepage Theories

3.1. Richards Model

The theory of transient two-phase flow in porous media was initiated by Richards [6]. Richards [6] constructed this model in three main steps. First, an explanation of negative water pressure or suction idea was invented to determine the state variable for the mathematical model. This state variable is

$$h=h_c+z=\frac{{\psi }_m}{{\rho }_{w}g}+z=\frac{u_a-u_w}{{\rho }_{w}g}+z$$

(7)

where h = the total water head (m); h_c = the soil suction head or negative water pressure head (m); z = the gravitational potential in elevation (m); ψ_m = matric suction (kPa); ρ_w = soil water density (kg/m³); g = gravitational accelerator (m/s²); u_a and u_w = soil gas and water pressure (kPa); and, here, u_a is assumed to be zero (u_a = 1 atm).

Second, Richards [6] convinced readers that capillary flow is a capillary and viscous predominated flow (laminar flow). Hence, the Hagen–Poiseuille law can be applied to the single capillary channel, and the empirical Darcy law can be used to describe capillary flow motion in a soil REV. The two-phase Darcy equation is

$$q=-k_s{k}_r\nabla h$$

(8)

where q = the volume of water passing through a unit area in unit time (m/s), k_s = the hydraulic conductivity of saturated soil REV (m/s); k_r = k_unsat/k_s the relative hydraulic conductivity; k_unsat = the hydraulic conductivity of unsaturated soil REV; h = the total water head (m).

Last, the two-phase Darcican form seepage equation, usually named the Darcy–Buckingham law or the two-phase Darcy seepage equation, is inserted into the continuity equation (mass conservation) to represent the flux. Finally, the Richards [6] model was constructed to become a nonlinear Partial Differential Equation (PDE) in diffusion form:

$$\nabla (k_s{k}_r\nabla h)=\frac{\partial \theta }{\partial t}=\frac{\partial nS}{\partial t}$$

(9)

where θ = soil volumetric water content of soil REV; n = porosity of the soil REV; S = saturation of soil REV; t = time (s); other parameters are the same as the notation above; ρ_w is neglected because of incompressible fluid assumption. This Equation is the entire formularization of the groundwater flow equation for saturated (k_r = 1 and S = 1) and unsaturated soil (0 ≤ k_r ≤ 1 and 0 ≤ S ≤ 1). By ignoring different terms in Equation (9), different types of groundwater PDE can be derived. The details could be sourced from Bear [7]. For steady-state and unsteady-state multiphase flow in rigid porous media, Equation (9) can be written in the following forms:

$$\frac{\partial }{\partial x}\left[k_s{k}_r(\theta )\frac{\partial h_c}{\partial x}\right]+\frac{\partial }{\partial y}\left[k_s{k}_r(\theta )\frac{\partial h_c}{\partial y}\right]+\frac{\partial }{\partial z}\left[k_s{k}_r(\theta )\frac{\partial (h_c+z)}{\partial z}\right]=0$$

(10)

$$\frac{\partial }{\partial x}\left[k_s{k}_r(\theta )\frac{\partial h_c}{\partial x}\right]+\frac{\partial }{\partial y}\left[k_s{k}_r(\theta )\frac{\partial h_c}{\partial y}\right]+\frac{\partial }{\partial z}\left[k_s{k}_r(\theta )\frac{\partial (h_c+z)}{\partial z}\right]=C_m(h_c)\frac{\partial h_c}{\partial t}$$

(11)

where C_m(h_c) is the specific capillary moisture capacity (m⁻¹), other notations are the same as above. Equation (10) is for steady-state flow conditions in the unsaturated soil, while Equation (11) is for the transient flow conditions.

Richards [6] model can be rewritten in three forms: the soil suction head (h_c), volumetric water content (θ), and two-state-variable mixing-based forms [1]. All of them can be transformed into the nonlinear diffusion form. They are individually listed in

Table 2. Richards model forms on different state variable basis [12].

Forms	PDE
hc base	$\frac{\partial }{\partial x}\left[k_s{k}_r(h_c)\frac{\partial h_c}{\partial x}\right]+\frac{\partial }{\partial y}\left[k_s{k}_r(h_c)\frac{\partial h_c}{\partial y}\right]+\frac{\partial }{\partial z}\left[k_s{k}_r(h_c)\frac{\partial (h_c+z)}{\partial z}\right]=C(h_c)\frac{\partial h_c}{\partial t}$	(12)
θ base	$\frac{\partial }{\partial x}\left[D(\theta )\frac{\partial \theta }{\partial x}\right]+\frac{\partial }{\partial y}\left[D(\theta )\frac{\partial \theta }{\partial y}\right]+\frac{\partial }{\partial z}\left[D(\theta )\frac{\partial \theta }{\partial z}\right]+\frac{\partial k(\theta )}{\partial z}=\frac{\partial \theta }{\partial t}$ where diffusivity $D(\theta )=\frac{k(\theta )}{C_m(\theta )};k(\theta )=k_s{k}_r(\theta );C_m(\theta )=\frac{\partial \theta }{\partial h_c}$	(13)
Mixing	$\frac{\partial }{\partial x}\left[k_s{k}_r(h_c)\frac{\partial h_c}{\partial x}\right]+\frac{\partial }{\partial y}\left[k_s{k}_r(h_c)\frac{\partial h_c}{\partial y}\right]+\frac{\partial }{\partial z}\left[k_s{k}_r(h_c)\frac{\partial (h_c+z)}{\partial z}\right]=\frac{\partial \theta }{\partial t}$	(14)

.

3.2. Hydraulic Conductivity Function

In order to solve the Richards equation, two functions have to be offered to account for the variation of hydraulic conductivity by volumetric water content, k(θ) and specific capillary moisture capacity varying with soil suction head, C_m(h_c). Until recent years, those two demands promoted the development of SWRF and Hydraulic Conductivity Function (HCF). However, since estimating HCF from SWRF, those two demands could be merged into SWRF only. Hence, a summary of HCF frequently applied in numerical simulation is listed in

Table 3. Hydraulic Conductivity Function (HCF).

Model Authors	Model Equations		Notations
Childs and Collis-George [95]	$k_r=\frac{{\displaystyle {\int }_{{\theta }_r}^{\theta }\frac{\theta -x}{\psi {(x)}^2}dx}}{{\displaystyle {\int }_{{\theta }_r}^{{\theta }_s}\frac{{\theta }_s-x}{\psi {(x)}^2}dx}}$	(15)	Se based statistical model; Fredlund, Xing and Huang [16] rewrote it into continuum form.
Gardner [96]	$k_r=\exp ({\alpha }_G{h}_c)$	(16)	Simplifying analytical solution derivation but having poor-fitting performance; hc based empirical model.
Brooks [13]	$k_r=\{\begin{array}{c}\begin{array}{cccc}1& & & h_c<h_c{}_{AEV}\end{array}\\ \begin{array}{cccc}{\left(\frac{{\alpha }_{BC}}{h_c}\right)}^{2+3n_{BC}}& & h_c>h_c{}_{AEV}& \end{array}\end{array}$	(17)	hc based empirical model; hcAEV = ψAEV/ρwg = air entry value in water head.
Brooks [13]	$k_r=S_e{}^{3+2/n_{bc}}$	(18)	Se based empirical model; BC SWRF inserted into Equation (18).
Mualem [17]	$k_r=S_e{}^{0.5}{\left[\frac{\left({\displaystyle {\int }_0^{{\theta }_w}\frac{d{\theta }_w}{h_c}}\right)}{\left({\displaystyle {\int }_0^{{\theta }_s}\frac{d{\theta }_w}{h_c}}\right)}\right]}^2$	(19)	Statistical model; Requiring well-developed SWRF inserted into Equation (19).
Van Genuchten [14]	$k_r=S_e{}^{0.5}[1-{(1-S_e{}^{1/m_{VG}})}^{m_{VG}}]$	(20)	mVG = 1 − 1/nVG; Se based model; VG SWRF inserted into Eqaution (20).
Fredlund, Xing and Huang [16]	$k_r=\frac{{\displaystyle {\int }_{\psi }^{{\psi }_r}\frac{\theta (y)-\theta (\psi )}{y^2}{\theta }^{\prime }(y)dy}}{{\displaystyle {\int }_{{\psi }_{aev}}^{{\psi }_r}\frac{\theta (y)-{\theta }_s}{y^2}{\theta }^{\prime }(y)dy}}$	(21)	Insert Fredlund and Xing [15] SWRF into Equation (21); Suction (ψ) based.

.

According to Table 3, available HCF can be divided into HCF based on soil suction/suction head and HCF based on effective saturation/effective volumetric water content. Equation (16) was usually used to derive the analytical solution of the Richards [6] model in one dimension [97] or two and three dimensions [11]. Nevertheless, this selection is just for simplifying the analytical derivation. Due to the hysteresis nature of SWRC, suction-based HCFs also have hysteresis behaviour [2]. The unique relation between relative hydraulic conductivity and effective saturation is assumed because the hydraulic paths are determined by fluid fractions filled into pores. Fredlund and Rahardjo [2] also presented a series of experimental data that wetting and drying HCFs collapse into a unique HCF. Therefore, S_e-based HCF is often considered for numerical solving Richard’s Equation. Leij, et al. [98] investigated the performance of a large number of HCFs against 346 S_e(h_c)-K(h_c) and 557 S_e(h_c)-K(S_e) datasets and recommended using Equation (19) for HCF and Equation (5) for SWRF. However, the hysteresis of S_e-based HCF was questioned by Childs [99] and Lu and Likos [1]: different hydraulic loading paths might not guarantee the same hydraulic conduits in one saturation. Moreover, those HCFs seem only satisfied with sandy soil with negligible soil deformation during the drying-wetting process. As Equation (19) is determined from PSD [17], soil deformation induced by purely hydraulic loading might result in the non-uniqueness of HCF. Thus, HCF encounters the same problem with SWRF for deformable soil. Moreover, Li, Luo, Li, Liu, Tan, Chen and Cai [50] pointed out that the time-dependence of SWRC and HCF coexist. In summary, the three issues previously mentioned for SWRC, including deformation, hysteresis, and time-dependence, could also be encountered for HCF.

3.3. Green-Ampt Model

The second semianalytical method to solve the transient process of water invading into the unsaturated soil is the one-dimensional Green-Ampt model. Green and Ampt [100] proposed a transient infiltration model by assuming that there is a sharp wetting front clearly separating the zone of saturation and dry zone. In Figure 3, this sharp wetting front replaces the water distribution along the vertical axis. The water content of the saturated zone is assumed to be effective porosity (θ_e). An initial water content (θ_i) is assigned to the dry zone. The wetting front infiltration rate, therefore, can be derived using Darcy law as

$$q=\frac{dQ}{dt}=\frac{({\theta }_e-{\theta }_i)dz}{dt}=k(1+\frac{h_0+h_s}{z})$$

(22)

where q = infiltration rate (m/s); t = time (s); k = effective hydraulic conductivity (m/s); h₀ = water head of ponding water above top soil surface (m); h_s = capillary suction head at wetting front (m); and z = depth of wetting front (m). Solving this ordinary differential equation gives cumulative infiltration displacement as

$$Q-(h_0-h_s)({\theta }_e-{\theta }_i)\ln \left(1+\frac{Q}{(h_0-h_s)({\theta }_e-{\theta }_i)}\right)=kt$$

(23)

where Q = cumulative infiltration displacement (m); other notations are the same as above. Equation (23) can be solved if ponding depth (h₀), effective hydraulic conductivity (k), and initial water content of dry zone (θ_i) are known.

Philip [102] stated that the Green-Ampt model is an exact solution of the Richards model if the diffusivity D(θ) is a Dirac-Delta function with non-zero water content in the saturated zone. Hence, the Green-Ampt model is just a simplified tool for approximating wetting front advancing depth. Kale and Sahoo [101] reviewed the Green-Ampt model with modified versions and concluded that the prediction of this model is sensitive to the effective hydraulic conductivity. The main advantage of this model is less demand of input parameters, compared to the Richards model in need of SWRF and HCF. For the Richards model, SWRF and HCF have to be experimentally measured in the laboratory for a long time due to equilibrium condition achievement. In contrast, the Green-Ampt model only needs effective hydraulic conductivity, ponding head, porosity, and dry zone initial water content. Each of these can be easily determined using standard soil mechanic tests in a shorter period. However, ignoring water content distribution in unsaturated soil constrains the accuracy of this model. Furthermore, due to the inhomogeneous nature of the ground, this averaged effective hydraulic conductivity can hardly be determined [101]. Therefore, the Richards model is better for predicting the unsaturated soil water flow, whereas the Green-Ampt model is more suitable for engineering applications under certain conditions.

3.4. Buckley and Leverett Model

The third method for transient two-phase flow seepage simulation was originated from petroleum engineering as the Buckley and Leverett [8] model. This method is almost identical to the Richards model and is also based on the Darcy–Buckingham seepage motion equation. The difference is that, instead of only considering the wetting phase fluid-soil water, both wetting and nonwetting phase fluid are considered in the continuity equations (mass conservation). So, the governing equations turn to be two sets of Equation (24) in

Table 4. The full governing equations for transient two-phase flow in porous media under equilibrium.

Equations	Forms
Mass balance of phase i	$\frac{\partial {\rho }_{i}n{S}_i}{\partial t}+\nabla ({\rho }_i{q}_i)=0$	(24)
Momentum balance of phase i	$q_i=-\frac{K_{r,i}(S_i)K}{{\mu }_i}(\nabla P_i-{\rho }_{i}g)$	(25)
Capillary pressure and Leverett J function	$P_c(S_w)=(P_n-P_w)={\sigma }_s\sqrt{\frac{n}{K}}J(S_w)$	(26)

where i = w, n (wetting and nonwetting fluid phase); ρ_i = density of fluid phase i; n = porosity of a designated REV of porous media; S_i = saturation of fluid phase i; t = time; q_i = volumetric flux for fluid phase i in unit area (Darcy flux); K_r,i(S_i) = relative permeability for fluid phase i; K = intrinsic permeability of a designated REV of porous media; µ_i = dynamic viscosity of fluid phase i; g = gravational accelerator; P_i = pressure of fluid phase i; P_c(S_w) = capillary pressure in the function of saturation S_w; σ_s = surface tension at two-phase interface; J(S_w) = dimensionless P_c, also defined as the Leverett J function by Buckley and Leverett [8].

. As this method is just another reformulation of the Richards model considering nonwetting phase fluids, the model performance still depends on SWRF and HCF. Last but not least, it should be noted that SWRF and HCF are redefined as capillary pressure-saturation P_c(S_w) and relative permeability-saturation K_r,i(S_i) constitutive relationships in petroleum engineering because of the nonwetting phase fluid considered.

According to the conventional theories reviewed above, it is not sophisticated to find that SWRF and HCF determine the modelling solution on each REV in the simulating domain. Therefore, the experimental determinations of SWRF and HCF play critical roles in the continuum-scale theory of transient two-phase flow seepage.

4. Review of Conventional Experiments in Continuum Scale

4.1. The Conventional Experiments of Soil Water Retention Curve

Two variables have to be measured to determine a Soil Water Retention Curve (SWRC). The first is water content or degree of saturation, while the second is soil suction.

Traditionally, the gravimetric water content can be easily measured using the oven drying method [103]. To further transfer gravimetric water content into volumetric content or degree saturation, the soil packing condition, including specific unit weight [104] and density [105], needs to be measured. However, as for deformable soil, dry density (also the porosity) will not keep constant, so shrinkage behaviour must be measured to determine the volumetric water content [106]. However, this is only necessary for the deformable soil such as clay, silt, or a mixture of clay and silt, but is not needed for coarse and sandy soil. Other methods also involve directly measuring a specimen using a scale (Method C of ASTM D6836-02 [107]), measuring water expelled from a specimen by a burette (Method A and B of [107]), and using water content sensors such as a Time Domain Reflectometer (TDR) [108].

The laboratory-scale suction measuring methods can be divided into two streams: suction control and suction measurement. The mainstream one is a suction control test based on Axis Translation Technique (ATT). The water content is measured for each static soil suction or capillary pressure under equilibrium during the suction control process. According to ASTM D6836-02 [107], there are two standard ATTs: the Buchner funnel, also named the hanging column method [109], and the pressure chamber method [6]. The only difference is due to their suction control approach. As for the hanging column method, the soil specimen is located on a saturated High Air Entry (HAE) ceramic plate embedded into a hanging column. The top of a soil specimen directly contacts the ambient atmosphere. A measuring burette is attached to the bottom output. Ensuring good moisture contact between the specimen and HAE ceramic disk, soil suction can be applied by changing the burette elevation in the water table. The soil suction head can be vertically measured via adjusting the burette to different water elevations. The water expelled from soil specimen is measured from the incremental volume in this burette. A sketch of the hanging column method is shown in Figure 4a.

As for the other ATT, suction can be controlled by increasing air pressure. The soil specimen is placed on the ceramic disk prelocated at the bottom of a gas-liquid proof chamber. Moreover, the bottom disk is attached to a water burette for measuring the volume of expelled water. Then, the high air pressure can be injected from the top of the cell. After each air pressure increase, a static soil suction or capillary pressure will be achieved until no further water is expelled from the specimen. The water content decrements can be either measured from the burette or measured the entire weight of the pressure chamber by a bench scale. An example of the pressure chamber method is shown in Figure 4b.

Although the Axis Translation Technique has been confirmed as a standard test for determining SWRC of unsaturated soil, it still has some limitations. The most critical issues for the Buchner funnel are the evaporation from the top of the specimen and a short range of measurable suction. Evaporation can be reduced by decreasing the ambient temperature. Nonetheless, the surface tension variation induced by temperature change has to be calibrated back to normal surface tension under laboratory temperature. In addition, it is questionable that the specimen might have different water distribution under various surface tensions induced by the temperature differences. Moreover, due to the limitation of vertical space in the laboratory, the standard Buchner funnel can only measure suction from 0–20 or 30 kPa (0–2 or 3 m in the suction head) [110]. Hence, this method is usually for sandy soil having a maximum suction of less than 30 kPa. This method is recently further developed in need of less vertical space using a vacuum control system [107,111]. This new method increases the upper limit to about 40 kPa [110]. However, this range is still insufficient for fine soil, so the pressure chamber is an irreplaceable apparatus in the laboratory’s SWRC determination.

According to ASTM D6836-02 [107], the pressure chamber method covers the suction range from 0 to 1500 kPa, and the recommended suction loading sequence is 10, 50, 100, 300, 500, 1000, and 1500 kPa. Nevertheless, according to the previous users’ experience, pore water will be quickly drained out by adjusting the low resolution of air pressure on the air pressure control panel. So, it is hard to collect enough data points to plot an SWRC of coarse soil using the pressure chamber. This operation will subsequently result in no coarse soil Air Entry Value (AEV) identification. Thus, the pressure chamber is more suitable for fine soil within a suction range of 10–1500 kPa. Vanapalli, Nicotera and Sharma [110] reviewed the ATT limitations, especially the pressure chamber method. The main problem leading to a discrepancy between the suction controlled by the pressure chamber and measured in the field is that the in situ pore water pressure is negative but positive in the pressure chamber. The high air pressure causes air diffusion into the soil water. Besides, the soil specimen in the pressure chamber will be compressed because of low air permeability in the initial step of air pressure loading. The theory for pressure chamber is only valid for a soil specimen under three assumptions. These assumptions include incompressible soil particles, interconnected pore-air voids, and continued air-water interphase [110]. Moreover, more compressed air bubbles will expand during and after passing through the HAE ceramic disk, with more air diffused into soil water. This phenomenon is because the low permeability of the HAE ceramic disk leads to a significant pressure drop between the pore water pressure in the chamber bottom and the positive hydrostatic pressure in the outlet burette. Therefore, air bubble flushing is essential under the bottom of the HAE ceramic disk [110]. The vital shortage of ATT is that the ceramic disk’s low permeability further induces a low speed of equilibrium achievement. After a one-step suction increase, the pressure difference between the air and back pressure (hydrostatic pressure in burette) is not the soil matric suction under transient capillary water flow conditions. The back pressure will be much smaller than the real pore pressure above the ceramic disk. Hence, even though ATT is accepted as the standard tests, it is still unable to determine the dynamic soil suction, negative water pressure, or capillary pressure under transient flow conditions.

Except for the mainstream test, the Axis Translation Technique, some other tests are recorded into Lu, et al. [112] and ASTM D6836-02 [107], such as the chilled-mirror hygrometer, contact and non-contact filter paper methods, centrifuge method, and tensiometer. The chilled-mirror hygrometer is a total suction control device utilizing humidity control techniques. By constantly controlling the ambient temperature to dewpoint temperature by chilled-mirror sensing technology, the total soil suction can be determined using Kelvin’s Equation (Equation (1)). This method is usually used for measuring a high suction range of fine soil, e.g., silt and clay (1000–450,000 kPa) [60]. Gubiani, et al. [113] compared this method with other methods to conclude that the lowest limit of a dewpoint meter should be 7000 kPa instead of 1000 kPa.

The contact and non-contact filter paper methods are indirect suction measuring methods. They measure the unsaturated soil suction from water content transfer from the sample to itself in an enclosed space. They cannot directly measure the soil suction, so a calibration between soil suction and moisture absorbed into filter papers should be accomplished before experiment implementation. Contact filter paper has intimate contact with ambient soil to associate absorbed soil moisture and matric suction. The non-contact filter paper is placed into the ambient environment of a soil specimen, so it can indirectly measure the total suction. Details of using this method can be sourced from Lu and Likos [1]. The advantage of the contact filter paper method is that it covers the entire suction range from 0 to 10⁶ kPa. However, this method is quite time-consuming because it usually takes ~7–10 days to reach the following equilibrium condition after the current one [1].

The centrifuge method is also a suction control technique in ASTM D6836-02 [107]. The soil matric suction can be controlled in the range of 0~120 kPa by varying angular velocities of the centrifuge. The advantage of this method is less temporal consumption, but the relatively high cost constrains its popularity.

A tensiometer is a powerful suction measuring technique initially developed by Gardner, et al. [114]. It is a water-filled shaft with a High Air Entry (HAE) ceramic tip at one end, and the other is connected to a pressure sensor. A micro-sized tensiometer (UMS T5 Tensiometer^® Umweltanalytische Mess-Systeme GmbH) and corresponding diagram are separately shown in Figure 5a,b. After inserting a tensiometer into a soil specimen, the water in the glass shaft will be absorbed into the soil by soil matric suction. As a result, this soil matric suction can be transferred to a negative pressure vacuuming the water in the glass shaft. This vacuum pressure, finally, is measured by the pressure sensor at the other end. This method only measures soil matric suction because the HAE ceramic tip is not solute impermeable. The accuracy of tensiometer measurement depends on the apparatus itself and users’ prudent installation. The key to measuring suction precisely is to set good contact between unsaturated soil and the HAE ceramic tip [60]. The measuring range of a conventional tensiometer is usually 0 to 100 kPa. Recently, the high-suction tensiometer measuring range has been increased up to 1500 kPa. Toll, et al. [115] reviewed the characteristics of newly developed high suction tensiometers, finding that the pressure transducer range (highest is 15 MPa) is much higher than the HAE ceramic tip (highest is 1.5 MPa). The suction range is constrained by the highest Air Entry Value (AEV) of this tip. Compared with all other methods, tensiometers have the following advantages:

• Fast and ease of installation (ensure tip contact condition)
• Applicability in both laboratory and field (in situ) condition
• Short responding time (less than 1 min for low suction capacity micro-tensiometer)
• Long-term measurement (ensure no drain out of water in the shaft)

According to the authors’ usage experience, a tensiometer should be carefully prepared before taking a measurement. Water in the shaft should be bubble-free above the highest suction encountered in the selected sample. Suppose a bubble clogs the glass shaft under high suction. In that case, the hydraulic path in the glass shaft will be significantly reduced, subsequently leading to a prolonged responding speed or even an incorrect reading. Due to the fast responding speed of the T5 Tensiometer, it is beneficial for studying dynamic soil suction or capillary pressure under the transient two-phase flow conditions.

4.2. The Conventional Experiments of Hydraulic Conductivity Function

Unsaturated soil hydraulic conductivity (k_unsat) is the hydraulic conductivity of unsaturated soil varying with suction or water content. It can be calculated by saturated hydraulic conductivity (k_s) and timed by relative hydraulic conductivity (k_r). In Section 3, “Steady-state and transient two-phase seepage theory”, the prediction of k_r from Soil Water Retention Function (SWRF) was already introduced. However, this approach is just an approximation of k_unsat to reduce experimental complexity and time consumption. Therefore, the most reliable result is still based on direct experimental measurement. Here, the experiments for measuring k_unsat will be summarised and discussed.

There are two types of experiments based on the two-phase flow seepage theory: steady-state and transient flow measurement. The conventional Axis Translation Technique (ATT) is often introduced in the steady-state experiment with some shortcomings. The centrifuge method is another one included in steady-state methods. The determination of k_unsat using ATT is relatively straightforward. Through applying suction by varying air chamber pressure and measuring outflow from the pressure chamber, the k_unsat can be calculated by the Darcy equation and then plotted against the corresponding average suction of soil specimen. ASTM D7664-10 [108] offers two subcategories of ATT: rigid wall (odometer) type ATT in Figure 6b and flexible wall (triaxial cell) type ATT in Figure 6a. The former is more suitable for coarse-grain soil, and the latter is more applicable for deformable soil containing fines. However, as the discussion of ATT shortages in the previous content in terms of air diffusion, air compressing specimen, and impedance of HAE ceramic disk with bubble expanding after pressure drop, ATT can only be used for an approximation of the Hydraulic Conductivity Function (HCF) of unsaturated soil under steady-state flow conditions [108]. According to ASTM D7664-10 [108], the air-water seepage is controlled by applying air pressure on the specimen top and back pressure (external water pressure) at the bottom underneath the HAE ceramic disk. It is recommended to insert tensiometers at different specimen depths, but soil moisture sensors with previously tested SWRC could replace them. The volumetric water content variation in a specific time interval can be integrated to calculate the water flux in this period. The hydraulic gradient can be provided between each pair of depths in the soil specimen where soil suction and water content are measured. With water flux and the chamber’s cross-section area and hydraulic gradient, the k_unsat can be calculated by the Darcy equation. More details can be sourced from ASTM D7664-10 [108].

In Lu and Likos [1] and Masrouri, et al. [117], other steady-state methods are introduced based on ATT. One of them is the constant head method, additionally setting a constant water head on the top of the soil specimen, as shown in Figure 6b [6]. In this method, suction is maintained during a two-phase flow seepage. The k_unsat is calculated by the Darcy equation under various suctions cautiously controlled by ATT. Then, the water content can be measured by the destructive method (the oven drying method) or the nondestructive method (e.g., outflow, soil moisture sensors). For the destructive method, more identical specimens have to be prefabricated. In this way, the HCF can be finally plotted with either a suction or water content basis. The other is the constant flow rate method developed by Olsen, et al. [118]. The core of the experimental setup is similar to the constant head method (suction controlled by ATT). Still, a highly complex steady-state flow control system is required to generate a constant flux in a triaxial cell [118]. After measuring suction, water content, and controlled flux, the k_unsat is also calculated by the Darcy equation.

The other steady-state method, except for ATT, is the centrifuge method. It uses a centrifugal force as a driving force instead of gravity [108]. The water is allowed to infiltrate into unsaturated soil specimens and finally reach the outflow reservoir. Higher horizontal field force can be produced by increasing angular velocity in centrifugal operation. With such a high horizontal field force, vertical gravity is neglected. Instead, the horizontal force is the only force-generating water pressure in the system. The hydraulic gradient can be calculated by angular velocity and seepage radius. Meanwhile, the outflow divided by a corresponding time interval determines the infiltration rate. The k_unsat is eventually calculated by the Darcy equation. Nevertheless, due to the high cost of the geotechnical centrifuge, it is not often available in most geotechnical laboratories. Masrouri, Bicalho and Kawai [117] recommended that the centrifuge method is only suitable for nondeformable soil specimens with a pore structure insensitive to the state of stress because of high net stress applied by a centrifugal operation.

For transient flow experiments shown in Figure 7, there are mainly three methods: the multistep in-outflow method [119], infiltration method [120], and instantaneous profile method (IPM) [10]. The first transient flow experiment still relies on ATT shown in Figure 6. By changing the suction within a small interval that is large enough for flux measurement, the volume change of pore water expelled for applied suction increment can be recorded for later calculating unsaturated soil water diffusivity using an analytical solution of the 1D Richards model in h_c-based diffusion form. Finally, with previously determined specific capillary moisture capacity (C_m), the k_unsat can be calculated by the notation in Equation (13). However, this experiment also has the same limitations for SWRC and steady-state HCF using ATT. Lu and Likos [1] summarised the six assumptions of the multistep outflow methods:

• Each suction increase interval must be small enough, so k_unsat can be assumed as a constant in this interval (which requires meticulous suction control);
• The relation between soil suction and water content is linear (but, in fact, it is not only nonlinear but also hysteresis);
• HAE ceramic disk does not cause any hydraulic resistivity (but it is a significant impedance, especially for high permeable sandy soil);
• Flow is just one dimension;
• The gravity effect can be ignored;
• The testing specimen is homogeneous and nondeformable (which is only available for sandy soil).

Therefore, it is not hard to find that transient flow experiments based on ATT are just methods for approximating k_unsat. Masrouri, Bicalho and Kawai [117] commented on this method because it owns simplicity and is good at mass control, but there are few reliable results compared with other methods. ASTM D7664-10 [108] also records this method as one transient flow method in the standard but mentions its limitations.

The second method [120] indirectly measures k_unsat by measuring the unsaturated soil diffusivity D(θ). The soil water content along advancing profile, distance, and the corresponding time interval can be recorded during soil water invading through the unsaturated soil. The Boltzmann variable (λ = xt^1/2), which is a function of both invading displacement (x) and duration (t), can be plotted against corresponding water content as a soil water content function of Boltzmann variable. The diffusivity D(θ) can be calculated from the Boltzmann transformation of the 1D Richards model (θ based diffusion equation). It needs integration of D(θ) function against the Boltzmann variable. With the addition of the previous measured SWRC, the k_unsat can be finally calculated by the notation in Equation (13).

The last method is the Instantaneous Profile Method (IMP), one of the transient methods recommended in ASTM D7664-10 [108]. This experiment is to replicate a natural vertical soil column in the laboratory. The soil suction and water content can be continuously measured under different wetting and drying paths by spatially inserting suction sensors and moisture sensors along the soil column. This method allows the flexible application of boundary conditions. ASTM D7664-10 [108] offers four types of hydraulic loading paths: infiltration, evaporation, drainage, and imbibition. This experiment can produce data for a dynamic nonequilibrium soil suction and water content under transient flow conditions. The water flux can be easily calculated by integrating soil water content variation by the corresponding duration. The suction head differences to corresponding vertical depths also determine the suction gradient. The Darcy equation is finally used to calculate k_unsat. Moreover, it is a good setup for investigating the dynamic effects in SWRC. However, Masrouri, Bicalho and Kawai [117] pointed out that this method lacks stress state and volumetric measurements. Therefore, sandy soil is often selected to conduct soil column tests of IPM to avoid deformation-induced stress variations.

5. Limitations of Conventional Theories and Experimental Methods

This section will initiate the dynamic nonequilibrium effects in soil water retention behaviour from the genesis of experimental findings. Then, the paradoxes embedded in the conventional two-phase flow seepage theory will be summarised according to prior studies. In the third part, microscale dynamic capillary physics will be reviewed from interfacial physics. Other possible reasons proposed for dynamic nonequilibrium effects will also be discussed.

5.1. The Experimental Exploration of Dynamic Nonequilibrium Effects

The unsaturated soil water flow was developed with a continuum-scale-constitutive relationship Soil Water Retention Curve (SWRC). However, many interfacial physics were omitted without rigorously experimental techniques on the microscale. The narrow view of investigating scales somehow limited our understanding of two-phase flow through a unit cubic-Representative Elementary Volume (REV). To calculate the two-phase flow in the soil domain made up of many soil REVs, a constitutive relationship, for instance, an intrinsic mechanical-geometry mathematical relationship such as a stress-strain elastic system, has to be assumed in this REV. SWRC is just another constitutive relationship to account for the inherent characteristic of water retention in porous media. By far, whether SWRC is an intrinsic characteristic or not has still been being continually debated in recent experimental studies [22,25,27,28,30,31,39,122].

Since traditional experimental techniques constrained early studies, ATT became the standard testing method for SWRC. With the development of soil sensors, Topp, Klute and Peters [21] studied the SWRC under unsteady state flow conditions using a tensiometer for suction measurement and gamma-ray system for moisture measurement and compared it with static SWRC. As tensiometers could provide an instant response of soil suction, the SWRC measured under transient flow conditions showed a significant discrepancy with static SWRC. An instance from Topp, Klute and Peters [21] is shown in Figure 8a. The dynamic suction is higher than static suction for water content, around 25% of static suction in a sand column [21]. However, the study of dynamic effect merely focused on drainage in early studies. This experimental observation strongly questioned the validity of the Richards model for simulating transient flow using SWRC without dynamic nonequilibrium effects. This issue was discovered by experimental studies on SWRC and the experimental validations of the Richards model. Hassanizadeh, Celia and Dahle [22] reviewed some early studies of Richards equation validation to find that the diffusivity depends on the speed of the wetting front but is not a material intrinsic. Some of the literature reviewed from the early study concluded that the diffusion equation could not be verified for soil moisture transport, and the application of Darcy’s law is questionable [23]. Rawlins and Gardner [123] experimentally found that the relationship between diffusivity and moisture content lost uniqueness for transient flow conditions (see Figure 8b). The speed dependence of diffusivity raised the concern that SWRC and HCF have dynamic nonequilibrium effects (see Figure 8a–d).

The dynamic nonequilibrium effects in drainage SWRC founded by Topp, Klute and Peters [21] were also experimentally measured from later studies [37,125,126,127,128,129]. These works used a similar setup, although, except Wana-Etyem [128], most of them solely focused on the drainage curve. Stauffer [37], Weller, et al. [130], and Weller and Vogel [131] also checked the relative permeability or HCF and found that it also had dynamic nonequilibrium effects. Smiles, Vachaud and Vauclin [126] confirmed that dynamic drainage curves have higher suction than static drainage curves, but they concluded that the imbibition curve had no observable dynamic effects. Nevertheless, this was later challenged by Wana-Etyem [128] that dynamic effects happen for both drainage and imbibition SWRC, as shown in Figure 8a,b. The nonequilibrium soil water retention behaviour was explored not only for hydraulic loading paths but also for unsteady-state evaporation by Bohne and Salzmann [132]. Hassanizadeh, Celia and Dahle [22] summarised their findings into the following points:

• The dynamic effects are not significant in fine-textured soil;
• The higher rate of water content variation, the more significant the dynamic effects;
• The dynamic effects are more significant in coarse-textured sand;
• The dynamic effects in primary drainage curves are more significant than the dynamic effects in main drainage curves.

5.2. The Theoretical Paradoxes of Transient Two-Phase Flow Seepage

In the original findings of dynamic nonequilibrium effects in SWRC, experimental results were the evidence for demonstration. However, there was no theoretical mathematical formulation to capture such behaviour due to a lack of understanding of two-phase thermodynamics. Moreover, most of the works were purely experimental exploration. With an aim to derive a mathematical equation counting dynamic nonequilibrium effects, the conventional theory must be critically reviewed to identify the shortages. Thus, the Richards model can be further improved for simulating transient soil water flow in unsaturated soil under nonequilibrium conditions.

The conventional theory was constructed in three steps: state variables, momentum conservation, and simplified momentum balance into mass balance. However, this nonlinear diffusion theory failing for the dynamic suction or water redistribution effect during a transient process is not due to the utilization of conservation laws but two other reasons: the poor expression of state variables and extending the groundwater diffusion of saturated soil to unsaturated soil water diffusion by simply adding parameters [133]. The state variables (e.g., positive water pressure or total hydraulic head) and simplified version of momentum balance (Darcy’s law) in the groundwater diffusion theory were solidly validated by many experimental studies [7]. However, the unsaturated soil water nonlinear diffusion theory does not have sufficient experimental studies to verify the reasonability of applying static capillary pressure as soil matric suction and two-phase Darcy seepage equation with addition of relative permeability concept. Hence, the paradoxes just happen in these two points.

As for the selection of state variables, the matric suction is always assumed to be equal to macroscale averaged capillary pressure as given in Equation (1). Then, this so-called negative water pressure (matric suction) is empirically correlated with wetting phase saturation by the concept of SWRC based on some early experimental findings. Gray and Hassanizadeh [133] criticised that the negative water pressure (the suction under zero reference atmosphere pressure) should not be simply determined by wetting phase saturation, while it should be rigorously derived from Equation of State (EOS) with the addition of phase saturation. If this is separated, there is a paradox to define if a negative pressure is a function of saturation or a function of wetting phase density and temperature.

The second paradox point Gray and Hassanizadeh [133] argued is the range of negative water pressure. With zero reference pressure of the atmosphere, soil matric suction should be bounded inside 0–1 atm. However, in the actual measurement, matric suction is usually above 100 kPa, which cannot be physically conceptualised. So, they stated that this high suction is energy that soil extract from a water molecule, and simply assuming the Equation between high matric suction and meniscus geometry cannot be proved. Bolt and Miller [134] reported that the pressure of most soil water in unsaturated soil remains positive. The significant values of soil water tension (high suction) are an artificial variable to account for the energy generated by the attraction of soil particles. Hence, Gray and Hassanizadeh [133] did not recommend that pore pressure and adsorption effects be lumped into one negative pore pressure.

The third paradox point exists in the relationship between matric suction pressure and suction head. Gray and Hassanizadeh [133] argued that the suction head could not be determined from matric suction divided by wetting phase density and gravitational accelerator. As the argument mentioned in the last paradox, this energy generated by the water-grain adsorption effect still depends on soil mineralogy and water molecule inter-attraction. Moreover, wetting phase density variation was not experimentally investigated. Gray and Hassanizadeh [133] suspected applying constant density into hydrostatic form to account for matric suction in elevation potential due to binding was causing significant density variation. Hence, simply using the hydrostatic form to represent the elevation potential of soil suction has never been physically proved.

The last paradox is on interface dynamics. The Young–Laplace Equation only provides a relationship between meniscus geometry and surface tension that is a microscale artificial tension variable accounting for the combination of solid surface adsorption and molecule attraction effects for the single capillary tube. It can only be experimentally validated without external actions. When capillary water starts to expel other nonwetting phase fluids or be expelled by such fluids, this Equation has been experimentally verified to be a failure [135,136,137,138,139,140,141,142]. These physical studies indicated that dynamic capillary pressure was influenced by advancing velocity, which is usually lumped into the dimensionless Capillary number (Ca) as an additional term to original static capillary pressure. Gray and Hassanizadeh [133] also proposed this paradox in their study. However, instead of looking into meniscuses in capillary conduits, they expected to add the specific interfacial area (air-water interface area per bulk volume) to quantify the interface variation because they expect a unique relationship between interfacial area and saturation. This involvement is also in their improved multiphase flow theory that will be introduced in Section 6.

5.3. The Physical Causes for Dynamic Nonequilibrium Effects

Before introducing the current advanced continuum-scale theories of two-phase flow in porous media, the recent research on dynamic capillary flow will be briefly reviewed to identify dynamic capillary theory’s state of the art in multiphase physics. The microscale mentioned here is not the molecular scale in chemistry, but rather a pore-scale among soil particles that can be simply conceptualised as a capillary tube. Although natural porous media has complex and pore structures as well as statistically distributed paths, each flow conduit can be simplified as a capillary tube to look into pore-size physics. Whether a theory is rigid should be verified by physics and experimental observation. In addition to the dynamic capillary effect in single flow conduit, other hypothesised reasons for dynamic nonequilibrium effects in porous media will be presented and discussed.

The Richards model using the two-phase Darcy equation to represent immiscible phases displacement into a soil body can be simplified as a two-phase movement in a series of capillary tubes. In this case, the assumptions, that viscosity force purely dominates fluid flow condition and conductivity of each phase is due to phase fraction in tubes, overlook the real flow phenomenon in even one single capillary tube. Compared with single-phase flow in either porous media or single tube, two-phase flow in a pore channel has more driving forces dominated flow condition in the transient state. Unfortunately, some early studies simply extend Hagen–Poiseuille’s law for studying two-phase displacement by adding capillary pressure in the form of the Young–Laplace Equation. The representative theory derivations in this way can be sourced from Bell and Cameron [143], Lucas [144], and Washburn [145]. Their eventually derived equation is

$$l=\sqrt{\left(\frac{{\sigma }_{s}r\cos \mathsf{\theta }}{2\mu }\right)t}=\sqrt{Dt}$$

(27)

where l = advancing distance (m); σ_s = surface tension (N/m); r = capillary tube radius (m); θ = contact angle; μ = dynamic viscosity (kg/ms or Pa·s); D = diffusivity (m²/s); t = advancing period (s). Equation (27) is the famous Lucas–Washburn (LW) equation for the mathematical formulation of dynamic capillary flow at the beginning of the twentieth century. There is an obvious overestimation of physics in this equation, which is the content angle is assumed to be constant during a moisture diffusion process. The Young–Laplace equation was assumed to be held under transient capillary flow but is not actually. Hoffman [136] stated that four forces control the capillary flow: viscous, inertial, liquid-gas interfacial, and liquid-gas-solid juncture. This set of combining forces dominates capillary flow depending on the testing system and flow rate. Thus, an additional term was added upon static capillary pressure to form the total capillary pressure [146]. This term was experimentally determined as K(Ca)^x in the function of capillary number Ca, and K and x are empirical coefficients. In order to continually use Hagen–Poiseuille’s law as a basis for two-phase flow in a capillary tube, some studies focus on verification of dynamic capillary pressure or dynamic capillary contact angles against Ca, so that the physically validated dynamic term can be added into Hagen–Poiseuille’s law [136,139,140,142,147]. However, the coefficients of such an empirical relationship are more determined by experimental results for the single capillary tube, which varies from study to study. It somehow constrains upscaling to a continuum-scale seepage equation for engineering application in natural porous media.

Zhmud, Tiberg and Hallstensson [138] reviewed the unphysical assumptions used for deriving the LW equation. Starting from the Newton dynamics equation assisted with capillary, gravity, viscous, and turbulent drag terms, Zhmud, Tiberg and Hallstensson [138] gave an experimental, analytical, and numerical analysis of dynamic capillary flow in a single capillary tube. It proved that the LW equation fails for prediction under short timescales, small viscosity limit and turbulent drag induced meniscus damping oscillations. As for the damping effect, Zhmud, Tiberg and Hallstensson [138] concluded that it is required for long capillary conduits while it can be neglected for short capillary tubes in a few times of tube radius. Zhmud, Tiberg and Hallstensson [138] explained the dynamic capillary flow process: at the beginning of the capillary flow, the wetting phase absorbed into a capillary conduit is dominated by capillary force, which violates the WS equation (x~t^1/2) and instead gives an x~t² relation; after viscous drag balancing capillary effect, two-phase flow reach to quasi-steady state obeying the LW equation; finally, flow is eased by gravity. Kim and Kim [137] reviewed recent physic studies on the dynamic capillary rise to find that the power number of advancing time gives different values for different packing beads. They suspected that this is because pores are not fully filled with wetting phase fluids. Moreover, not only for a short instant time-scale, but the LW equation fails to match experimental results for a considerable long enough time [137].

Despite the dynamic advancing interface theory and corresponding empirical function of Ca with endless fitting into various experimental results, the energy method is another method to account for dynamic capillary pressure. Yang, Krasowska, Priest, Popescu and Ralston [135] developed a modified LW equation by considering capillary driving force generated from an unleashed free energy. In this way, the capillary driving force can be divided into a static capillarity given by the Young–Laplace equation and a dynamic free energy term, referring to the thickness of the meniscus. Inserting dynamic capillary pressure terms into the Hagen–Poiseuille law, Yang, Krasowska, Priest, Popescu and Ralston [135]’s model shows good agreement with their experimental results. Moreover, the dynamic energy term could be used for continuum-scale capillary dynamics in porous media. It also offers a chance for upscaling from the thermodynamic basis [148]. Hence, the continuum-scale dynamic two-phase flow in porous media commenced with theories built upon physics rather than experimental empirical findings.

The aforementioned multiphase physics in a single capillary conduit merely focus on the microscale. Other reasons could cause the dynamic nonequilibrium effect for a soil specimen at the continuum scale. Those were proposed without experimental validations on the microscale. Diamantopoulos and Durner [49] summarised those reasons: immiscible two-phase interface reconfiguration, dynamic contact angle, air entry, air and water entrapment and blockage, hydrophilicity variation with time, and micro and macroscale inhomogeneities. The interface configuration and dynamic contact angle belong to the capillary dynamics previously reviewed. Besides, organic substances in natural soil inducing the time-dependent hydrophilicity are beyond the scope of Geotechnics and refer more to Agriculture. It sometimes also happens in the deep reservoir because the wettability shifts in the rock mass after surfactant flooding.

The reason for air entry and entrapment challenged the instantaneous atmosphere pressure for soil gas assumed in the Richards model. However, it can be easily tackled by the Buckley and Leverett [8] model, considering nonwetting phase fluid with experiments of two-phase displacement in oil/gas reservoirs. Nevertheless, the discontinuity of gas or nonwetting phase cannot be simulated due to scale constraints.

The water entrapment and blockage are due to a combination of microscale local heterogeneity (e.g., ink-bottle effect) and film water flow after the main capillary flow drained out (Marangoni effect). This effect can be frequently found as tears of wine on the wine glass and is a capillary convention process. After shaking wine in the glass, most wine drops while other very thin wine films are left on the glass and finally form many small liquid spots on the wine glass. The same phenomenon also occurs in the soil matrix during a nonequilibrium drainage process.

As for the macroscale heterogeneity, preferential water or wetting phase flow in cracking ground or fissures, fractures and faults in deep reservoirs can occur, leading to finger flow for gas or nonwetting phase fluids. This issue can be categorised into large-scale inhomogeneity, which can be resolved by assigning different intrinsic hydraulic properties in the simulating domain. It can also be tackled by applying the dual-porosity and dual permeability model. In this model, two porosities (specific storages) and two permeability are set for a single REV of soil to separately count seepage through homogeneous porous media and nonequilibrium seepage in preferential flow paths [149].

Similar fingering flow also happens in a single REV of unsaturated soil when microscale local heterogeneity is considered. Mirzaei and Das [35] investigated the influence of microscale heterogeneity on dynamic effects using numerical simulation. They concluded that the intensity of microscale heterogeneity magnifies the dynamic effects. Joekar-Niasar and Hassanizadeh [150] generated dynamic effects using pore-network numerical modelling without considering dynamic capillarity in capillary conduits. Therefore, they concluded that local heterogeneity, such as the ink-bottle effect, also significantly contributes to the dynamic effects.

In summary, there are various physical reasons for dynamic nonequilibrium effects. Each of them can be separately investigated using physical experiments or numerical simulations. However, each cause contributing to the dynamic effects cannot be quantified and compared, especially for porous media. Furthermore, when the microscale state variables are upscaled to the continuum scale, all reasons will be combined together and finally lumped into the REV-scale state variables. Therefore, theoretically, it is impossible to derive the dynamic effects from a single microscale physical cause. Instead, multiple reasons in multiscale should be involved in theoretical development.

6. Advanced Theories of Transient Two-Phase Flow Seepage

To the authors’ best knowledge, four theories have been developed to simulate nonequilibrium transient two-phase flow in porous media since the 1960s. They are, separately, the theories of dynamic fluids redistribution [151,152], dual-fraction with dynamic fluids redistribution [153], dual-porosity and dual permeability [154,155,156], and dynamic nonequilibrium capillary pressure [43,148,157]. In this section, those four theories will be summarised with discussions regarding their physical basis.

6.1. The Theory of Dynamic Fluids Redistribution

In petroleum engineering, Barenblatt [151] was one of the pioneers finding that an aiming saturation could not be achieved immediately with a fast two-phase fluids displacement in porous media; instead, this process takes a finite relaxation period (τ_B). Based on the experimental phenomenon and this saturation relaxation assumption, the first dynamic theory is proposed as the difference between the equilibrium saturation (S_equ) and dynamic nonequilibrium saturation (S_dyn) is equal to a saturation relaxation term as

$$S_{dyn}-S_{equ}={\tau }_B\frac{\partial S_{equ}}{\partial t}$$

(28)

where τ_B is a fluid redistribution time (s) and t is time (s). So, the capillary pressure-saturation relationship (P_c-S) is a capillary pressure function of S_dyn as

$$P_n-P_w=P_c(S_{dyn})=={\sigma }_s\sqrt{\frac{n}{K}}J(S_{dyn})$$

(29)

where P_c is dynamic capillary pressure (kPa), other notations are the same as Table 4. Since relative permeability K_r(S) is a function of saturation, therefore, the dynamic relative permeability can be given as

$$K_{r,i}(S_{dyn})=K_{r,i}\left(S_{equ}+{\tau }_B\frac{\partial S_{equ}}{\partial t}\right)$$

(30)

where K_r,i is the relative permeability of fluid phase i = w, n (wetting and nonwetting phase). By adding Equation (28) into the mass balances and two-phase Darcy equations in Table 4, the theory of nonequilibrium two-phase flow seepage was constructed in petroleum engineering.

This model has been solved by Barenblatt, et al. [158], which shows that the width of stabilised displacement front is not a linear relationship with the inverse of advancing velocity. Details of a modified version of their original model and modelling result can be sourced from Barenblatt, Patzek and Silin [158], and numerical solvers are available in the reference list of this work. Although the modelling results agree with experiments in their studies, the phase distribution term does not have a physical basis. Still, it is only a phenomenological term to account for experimental findings. Moreover, in their studies, the two-phase Darcy equation is still assumed to be valid against the physics of various flow regimes for dynamic capillary flow in previous review.

In soil hydrology, Ross and Smettem [152] also developed a theory of nonequilibrium water flow in unsaturated soil beyond the Richards model. It shares a similar physical concept (Equation (28)) with Barenblatt [151] as

$$\frac{\partial \theta }{\partial t}=f(\theta ,{\theta }_{equ})=\frac{{\theta }_{equ}-\theta }{{\tau }_R}$$

(31)

where θ = actual soil water content, t = time, θ_equ = equilibrium soil water content and τ_R = an equilibrium time constant. The differences from previous petroleum engineering theory include neglecting the nonwetting fluid phase-soil gas and treating soil water redistribution as an independent process from equilibrium transient unsaturated soil water flow described by the Richards model. As Equation (31) does not need to be in Soil Water Retention Function (SWRF) and Hydraulic Conductivity Function (HCF), this theory is much simplified, and a semi-analytical approximation can, therefore, be derived as an asymptotic solution:

$${\theta }^{t+1}={\theta }^t+({\theta }_{equ}^{t+1}-{\theta }^t)[1-\exp (-\frac{\Delta t}{{\tau }_R})]$$

(32)

where t and t + 1 indicate the discrete-time steps and next step. Ross and Smettem [152] successfully modelled infiltration in a field-aggregated soil with the significant preferential flow using this nonequilibrium model and also proposed a dural-fraction model:

$$\theta ={\theta }_{equ}+{\theta }_{dyn}$$

(33)

where θ_dyn = dynamic nonequilibrium soil water content. In this model, Equation (32) approximates θ_dyn, and the Richards model can solve θ_equ. However, they did not develop a systematic method with experiments to split the porosity or total soil water content.

6.2. The Theory of Dual-Fraction with Dynamic Fluids Redistribution

Diamantopoulos, Iden and Durner [153] continued the theoretical work left from Ross and Smettem [152] to arrive at a dual-fraction governing equation:

$$\nabla (k_s{k}_r({\theta }_{equ})\nabla h)=f_{equ}\frac{\partial {\theta }_{equ}}{\partial t}+f_{dyn}\frac{\partial {\theta }_{dyn}}{\partial t}=(1-f_{dyn})\frac{\partial {\theta }_{equ}}{\partial t}+f_{dyn}\frac{{\theta }_{equ}-{\theta }_{dyn}}{{\tau }_D}$$

(34)

where k_s = hydraulic conductivity of saturated soil, k_r = relative hydraulic conductivity, h = total water head, θ_equ = equilibrium soil water content, θ_dyn = dynamic nonequilibrium soil water content, f_equ = θ_equ/(θ_equ + θ_dyn) fraction of equilibrium soil water content, f_dyn = θ_dyn/(θ_equ + θ_dyn) fraction of dynamic nonequilibrium soil water content, and τ_D = soil water equilibration time.

In addition to infiltration simulation by Ross and Smettem [152], Diamantopoulos, Iden and Durner [153] succeeded in the simulation of multistep outflow by numerical solutions of Equation (34). Using Hydrus-1D software code from Šimůnek, Jarvis, et al. [149], Diamantopoulos, Iden and Durner [153] also developed a systematic method to determine τ_D and f_dyn using an inversion analysis of the dual fraction model against experimental data.

6.3. The Theory of Dual-Porosity and Dual-Permeability

The dual-porosity theory was initially created by Philip [154] to simulate soil water flow in a structured soil containing two domains: the inter- aggregate and intra-aggregate pore matrices. The Richards model can describe the water flow in the inter-aggregate pore matrix, while the water in the intra-aggregate pore matrix is immobile. The water exchange between inter-aggregate and intra-aggregate pore matrices is described by introducing an additional term of water transfer rate. However, the dual-porosity model was not developed for nonequilibrium transient unsaturated soil water seepage due to one domain described by the Richards model and stagnant flow in the other domain.

Based on this two-domain concept, the theory of dual-porosity and dual-permeability was later invented by Gerke and van Genuchten [155] to model the nonequilibrium unsaturated soil water flow. There are two improvements added to the dual-porosity theory. First, the soil water flow in both domains can still be simulated using the Richards model. Second, the water exchanging term between two domains is divided by two fractions. The governing equations include two sets of the Richards equation:

$$\begin{array}{l}\nabla (k_{\mathrm{m}}(h_m)\nabla h_m)-\frac{{\Gamma }_w}{w_m}=\frac{\partial {\theta }_m}{\partial t}\\ \nabla (k_{im}(h_{im})\nabla h_{im})+\frac{{\Gamma }_w}{w_{im}}=\frac{\partial {\theta }_{im}}{\partial t}\end{array}$$

(35)

where k_m = hydraulic conductivity for inter-aggregate pore matrix, k_im = hydraulic conductivity for intra-aggregate pore matrix, h_m = water head in inter-aggregate pore matrix, h_im = water head in intra-aggregate pore matrix, θ_m = soil water content in inter-aggregate pore matrix, θ_im = soil water content in intra-aggregate pore matrix, Γ_w = soil water exchange between two domains, w_m = θ_m/(θ_m + θ_im) fraction of soil water in inter-aggregate pore matrix, and w_im = 1 − w_m = θ_im/(θ_m + θ_im) fraction of soil water in intra-aggregate pore matrix. In later work, Gerke and van Genuchten [155] also evaluated the first-order water transfer term between two domains as Γ_w = α_w(h_m − h_im), α_w = a first-order mass transfer coefficient. Šimçnek, et al. [159] successfully applied this theory to simulate an upward imbibition experiment. The newly proposed parameters were also determined using inverse modelling. A comprehensive review of the two-domain model simulating the nonequilibrium unsaturated soil water flow can be sourced from Šimůnek, Jarvis, Genuchten and Gärdenäs [149].

There is common incapability amongst the aforementioned three theories. Any of them cannot simulate the nonequilibrium soil suction and capillary pressure. Therefore, they are only applicable for agricultural irrigation but less applicable in Geotechnics because the hydro-mechanical behaviour of unsaturated soil cannot be predicted in terms of wetting-induced soil shear strength reduction and soil collapse. Besides, all newly developed physical concepts with corresponding coefficients and parameters cannot be determined straightforwardly using a specifically designed experiment but only achieved using an inverse modelling technique.

6.4. The Theory of Dynamic Nonequilibrium Capillary Pressure

The last theory was developed based on thermodynamics. Kalaydjian [157] and Hassanizadeh and Gray [148] all started from microscale multiphase physics in single capillary conduit and finally ended up by upscaling to dynamic nonequilibrium capillary pressure equation for continuum-scale porous media REV as

$$P_c{}^{dyn}-P_c{}^{stat}=-\tau \frac{\partial S_w}{\partial t}$$

(36)

where P_c^dyn = dynamic capillary pressure (Pa), P_c^stat = static capillary pressure (Pa), τ = dynamic coefficient (Pa·s) as the same as the unit of viscosity, S_w = wetting phase saturation, and t is the advancing duration (s). Those two works have commons for deriving the macroscale dynamic capillary pressure. It was achieved by introducing the Gibbs energy for each fluid phase and later concluded that the Helmholtz free energy could not be neglected under the nonequilibrium capillary flow condition [148]. The difference between these two theoretical works is their upscaling methods. Hassanizadeh and Gray [148] used the capillary tube as a unit microscale system and upscale the units by volume averaging method, which is based on the concept of REV [160]. Instead, Kalaydjian [44] used the weight function method. The key findings from those theoretical works were that the specific interfacial area (total interfacial by the volume of porous media REV) and saturation variation play critical roles in Helmholtz energy, which causes the difference between equilibrium and nonequilibrium capillary pressures. Instead of merely considering interface dynamics in single capillary conduit in microscale, the saturation as a variable on REV-scale should also be involved in the REV-scale Helmholtz free energy terms because of local heterogeneity issue in single REV (e.g., ink-bottle effect and Haines jumps [109]).

Both studies derived the mass balance, momentum balance, and energy balance for each fluid phase and the interface. By constructing the energy transfer between the fluid phase and interface, the system is finally enclosed to provide a similar framework of conventional two-phase flow theory with some additional parameters accounting for so-called Helmholtz free energy released in a dynamic capillary flow phenomenon [43]. The entire theory is very physical-based and much beyond the scope of engineering practice. The advanced multiphase flow theory system simplified from Hassanizadeh and Gray [43] is summarised in

Table 5. A simplified system of advanced theory for dynamic two-phase flow in porous media [43].

Equations	Forms
Mass balance phase	$\frac{\partial {\rho }_{i}n{S}_i}{\partial t}+\nabla ({\rho }_i{q}_i)=0$	(37)
Momentum balance	$q_i=-\frac{S_i{}^{2}K}{{\mu }_i}\left[(\nabla P_i-{\rho }_{i}g)+\frac{{\lambda }_{ii}}{a_{wn}}\nabla a_{wn}+\frac{{\Omega }_i}{S_i}\nabla S_i\right]$	(38)
Dynamic capillary pressure	$P_c{}^{dyn}-P_c{}^{stat}=(P_n-P_w)-P_c{}^{stat}=-\tau \frac{\partial S_w}{\partial t}$	(39)
Equation of state	$P_c{}^{stat}=P_c{}^{stat}(S_w,a_{wn},T)$	(40)

where i = w, n (wetting and nonwetting phase), ρ_i = the density of fluid phase i, n = the porosity, t = time, S_i = the saturation of fluid phase i, q_i = the volumetric flux in unit area (discharge velocity), K_i = the intrinsic permeability of porous media, µ_i = dynamic viscosity of fluid phase i, g = gravitational accelerator, P_i = the pressure of fluid phase i, a_wn = the specific interfacial area (interfacial area per REV), λ_ii, Ω_i= the material coefficients (S_w·λ_wn=S_n·λ_nw), P_c^dyn = dynamic capillary pressure, P_c^stat = static capillary pressure, τ = dynamic coefficient, T = the absolute temperature.

to promote the application.

However, the new parameters such as specific interfacial area and two new material coefficients can hardly be measured straightforwardly using currently available experiments but can only be processed with inverse modelling. Thus, even though this theory has a rigid physical basis, its application requires improvement of measurement techniques and numerical techniques to determine coefficients of this theory and test the performance of prediction. For further simplifying this theory for unsaturated soil water transient seepage, Hassanizadeh, Celia and Dahle [22] adopted the zero-reference of soil gas pressure and, therefore, arrived at a modified Richards model with a third-order term:

$$\begin{array}{l}\nabla (k_s{k}_r(\theta )\nabla h_{dyn})=\nabla (k_s{k}_r(\theta )\nabla h_{equ})+\nabla (k_s{k}_r(\theta )\nabla (\frac{\tau }{n\rho g}\frac{\partial \theta }{\partial t}))=\frac{\partial \theta }{\partial t}\\ h_{dyn}=h_{c,dyn}+z=h_{c,equ}+\frac{\tau }{n{\rho }_{w}g}\frac{\partial \theta }{\partial t}+z=h_{equ}+\frac{\tau }{n{\rho }_{w}g}\frac{\partial \theta }{\partial t}\end{array}$$

(41)

where k_s = hydraulic conductivity of saturated soil, k_r = relative hydraulic conductivity, θ = soil water volumetric content, h_dyn = total dynamic nonequilibrium water head, h_c,dyn = total dynamic nonequilibrium soil suction head (Gibbes energy in head), z = elevation potential head, h_c,equ = equilibrium soil suction head (static capillary pressure in head), h_equ = total equilibrium water head (total hydrostatic water head), τ/nρ_wg∙∂θ/∂t = dynamic nonequilibrium overpressure head (Helmholtz free energy in head), τ = dynamic coefficient, n = porosity, ρ_w = soil water density, g = gravitational accelerator, and t = time. In this theory, the material coefficients for Helmholtz free energy were reduced to the dynamic coefficient only. All other previously defined SWRF and HCF can still be applied to the simplified theory. It also alleviates the experimental effort so that the dynamic coefficient can be directly determined using a soil column test or pressure cell apparatus without a ceramic disk underneath.

In comparison to the previous three theories, the thermodynamic-based theory owns more rigorous physical origins. It provides a continuum-scale vision, looking into the dynamic nonequilibrium capillary effect. It also gives a Darcian form flow motion equation with clear identification of the contribution of driving force released from interfacial energy transformation. In fluid redistribution theory, two-phase Darcy law is still applicable with the requirement of relative permeability for both two phases in Buckley and Leverett [8] model. However, the contribution of dynamic capillary pressure lumped into HCF or relative permeability overlooks flows regimes that occurred in the transient seepage process. Moreover, the dynamic SWRC and HCF cannot be readily determined using currently available experimental techniques. Thus, there is always a conflict of interest between rigorous physical-based theories and available experimental methods.

7. Novel Experimental and Numerical Contributions in Multiscale

In the latest decade, the mainstream of experimental or numerical studies is validating the thermodynamic-based theory and investigating the testing methods with inverse modelling of previously reviewed theories. They can be separated into mainly three groups:

• Macroscale 1D soil column experiments supported by pressure and moisture sensors with inversion analysis of the Richards and modified models;
• Microscale physical pore network model supported by imaging technique;
• Microscale numerical experiments using Pore Network Model (PNM, solving Poiseuille form capillary flow equation in artificial pore network), Direct Numerical Simulation (DNS, solving Navier–Stoke equation in artificial beads package) and Lattice Boltzmann Method (LBM, solving discretised Boltzmann equation in virtual particles package).

7.1. The State-of-Art of the Continuum-Scale Experiments

Besides the early studies on observing dynamic nonequilibrium effects in SWRC and validating the conventional models, the one-dimensional soil column and core flooding tests are currently being used to investigate the influential factors in dynamic coefficient and the improved Richards or two-phase seepage models coupled with additional dynamic capillary pressure. The demonstrations of the instrumented soil column and core flooding tests are shown in Figure 9.

Many soil column and core flooding tests are available from many experimental works. As shown in Figure 9a–c, several common instruments are installed on the specimen. Compared to conventional pressure cell tests determining the capillary pressure/soil suction by external water back pressure under the ceramic plate and fluid content by accumulative in/outflow, the instrumented setup highly demands temporal and spatial data logging using pressure transducers and dielectric/electrical sensors. Upgrading the air-water to the two-phase testing system can be achieved by setting an extra reservoir for nonwetting phase fluid. Moreover, the tensiometers for soil suction ought to be upgraded to high-pressure transducers. The soil/rock specimen dimension can be shortened to less than 10 cm or extended to 100 cm or even up to 200 cm long.

As for the shorter measuring window, one or two pressure and moisture sensors are required to detect fluid content and pressure for a single Representative Elementary Volume (REV) of soil/rock specimen. Due to the ultra-low permeability of ceramic plate at the bottom of pressure cell, nonequilibrium data cannot be recorded as reliable. Nylon and semipermeable membranes are often utilised to replace the ceramic plate in standard tests in order to overcome this inability. Such membranes own unique properties in high permeability for the wetting phase fluid but impermeability for the nonwetting phase fluid. As for the longer measuring window, the fully saturated soil/rock specimen under the partially saturated zone can prevent the nonwetting phase fluid from breakthrough more like in situ conditions in the field. Moreover, with more sensors inserted along with the measuring profile, both REV-scale and vertical profile can be investigated. Nevertheless, compared to shorter tests, the larger one has more time consumption and labour cost in terms of experimental setup complexity.

The hydraulic boundary conditions are as same as the standard tests. It includes one-step, multistep in/outflow, fluctuated water table, constant infiltration/imbibition flux from top or bottom, etc. For example, a vibrating water head overflowing system (see Figure 9b) can implement most hydraulic pressure boundary conditions. Likewise, the hydraulic pump in Figure 9c can control the flow to achieve constant flux boundary conditions for infiltration and imbibition tests. All sensors are connected to the datalogger for consistent data recording, so the temporal variation of state variables can finally be studied for transient flow conditions.

The soil column tests for the air-water system are general in sensor selection and boundary control in many experimental works. Therefore, it leads to few comprehensive reviews of soil column setup in the literature but more focus on reviewing experimental observations [49,164]. In contrast, due to the complexity of setting up the two-phase testing system in terms of high pressure and temperature, various two-phase flow core flooding and carbon dioxide sequestrations were sufficiently reviewed by Sun, et al. [165] and Li, Luo, Li, Liu, Tan, Chen and Cai [50]. With an aim to investigate influential factors in dynamic capillary effects, those experimental approaches were often applied with the numerical solution of advanced theories.

7.2. The Influential Factors in Dynamic Nonequilibrium Effect

Stauffer [37] intuitively proposed the dynamic capillary pressure in Equation (36) based on the earlier experimental analysis, far before the theoretical derivation. Furthermore, Stauffer [37] presented an empirical relationship between intrinsic soil properties and dynamic coefficient:

$$\tau ={\alpha }_{stauffer}\frac{n\mu }{K{n}_{BC}}{\left(\frac{{\alpha }_{BC}}{\rho g}\right)}^2$$

(42)

where α_stauffer = the constant proposed by Stauffer [37] with a value of 0.1 for the air-water system, n = the porosity, K = the intrinsic permeability for the saturated soil, µ = the dynamic viscosity of soil water, ρ = the density of soil water, g = gravitational accelerator, and α_BC and n_BC = the fitting parameters of the SWRF in Equation (3). By far, the influential factors to dynamic effects have still been continuously studied by Das and Mirzaei [166] on saturation dependency, Abidoye and Das [34] on scale effect, Hanspal and Das [41] on temperature effect, Goel and O’Carroll [39] on fluid viscous effects, Mirzaei and Das [35] on local heterogeneity influence, Manthey, Hassanizadeh and Helmig [40] on permeability dependency, O’Carroll, Mumford, Abriola and Gerhard [122] on wettability dependency, and Mirzaei and Das [167] on hysteretic dynamic effects, etc.

The saturation dependency: The dynamic coefficient was originally defined as a constant independent of saturation. However, the linear relationship between dynamic coefficient and saturation was firstly found by O’Carroll, Phelan and Abriola [28] using multistep outflow experiments. Later, Sakaki, O’Carroll and Illangasekare [30] also explored the saturation dependency in dynamic effects on primary drainage and main imbibition. This work found that the dynamic coefficient increases with decreasing wetting phase saturation. Moreover, the linear relationship between them was updated to become a log-linear relationship. This log-linear relationship was later reconfirmed by another recent experimental contribution on hysteretic dynamic effects from Zhuang, Hassanizadeh, Qin and de Waal [47]. Meanwhile, Das and Mirzaei [166] used a 1D soil column set up to study the dynamic coefficient and found that dynamic coefficient is not a linear function of saturation but a nonlinear function in which dynamic coefficient can only be treated as a constant within high saturation 70–100%. When saturation is lower than 60%, the dynamic coefficient increases nonlinearly with saturation [166]. Conflictingly, the most recent soil column test results from Luo, Kong, Ji, Shen, Lu, Xin, Zhao, Li and Barry [161] showed the uniqueness of dynamic coefficient-saturation relationships on both drainage and imbibition for specific sand under a given period of water table fluctuation.

The scale effect: Abidoye and Das [34] applied dimensional analysis to nine parameters (gravity g, isotropic intrinsic permeability K, bubbling pressure ψ_AEV, the domain volume representing domain scale V, fluid density ρ, fluid viscosity μ, saturation S, porosity n, pore size distribution index of SWRF in Equation (3) n_BC), which are reported as essential variables in the determination of dynamic coefficient, to derive a nonlinear relationship between two dimensionless groups as

$${\prod }_1=\frac{\tau \sqrt{g}}{K^{0.25}{\lambda }_{AEV}}=a{\left[{\prod }_3\right]}^b=a{\left[\frac{V{K}^{-1.5}\rho \mu }{nS{n}_{BC}}\right]}^b$$

(43)

where ∏₁ = the first dimensionless groups, ∏₃ = the third dimensionless groups, and a and b are fitting coefficients. As for the experimental results of Das and Mirzaei [166], a = 9 × 10⁻¹⁴ and b = 1.31, other notations are given in the previous content. Prediction from Equation (43) shows good agreement with the experimental results not only from Das and Mirzaei [166] but also Bottero [168]. Hence, this might be the first mathematical form to quantify the dynamic coefficient impacted by nine other essential variables, including the domain scale effect. Bottero, Hassanizadeh, Kleingeld and Heimovaara [33] and Abidoye and Das [34] investigated nonequilibrium capillary effects at various scales, thereby concluding that the dynamic coefficient increases with the observation scale. The same conclusion was again validated by Abidoye and Das [169] using an artificial neural network modelling approach. Later, Goel, et al. [170] further studied the scale dependency of dynamic relative permeability. This work found no observing scale dependency of dynamic wetting phase relative permeability, but the dynamic nonwetting phase relative permeability slightly increases with domain size decrease. In addition, the location dependency was also unveiled in this work. According to the data from Goel, Abidoye, Chahar and Das [170], with the measuring zone moving from top to bottom, the dynamic wetting relative permeability decrease but the dynamic nonwetting phase relative permeability increase. As a result, Goel, Abidoye, Chahar and Das [170] concluded that the location dependency varies according to the location of the fluid injection point.

The temperature effect: Hanspal and Das [41] carried out a numerical simulation of unsteady capillary flow in porous media between 20 °C and 80 °C. Their results showed that dynamic coefficients were nonlinear functions of temperature and saturation, and the dynamic coefficient increases with a temperature increase [41]. Meanwhile, Civan [171] also investigated the temperature effect and reconfirmed the conclusion from Hanspal and Das [41]. However, in spite of this numerical exploration of the temperature effect, any experimental investigations on the temperature dependency were quite rare and less found in the literature by far.

The fluid viscosity dependency: Goel and O’Carroll [39] experimentally studied the variation of dynamic coefficient impacted by the viscosity of non-wetting phase fluids using a 1D sand column drainage test. In their study, three important points are mentioned: (1) there is a delay response of the tensiometer due to permeability of ceramic cup, which implies using a high conductivity tensiometer to reduce response postponement during the dynamic experiment; (2) dynamic coefficient decreases for non-wetting phase fluids having smaller viscosity; (3) their work is the first experimental study used to validate previous numerical experiments for viscosity effect and provided accurate data against some contradictory conclusions from numerical studies. The primary purpose of their work is to validate the conclusion that the dynamic coefficient can be enlarged with increasing the effective viscosity (µ_eff = µ_nS_n + µ_wS_w) proposed by Barenblatt, et al. [172]. Li, et al. [173] also drew a similar conclusion for the fractured tight reservoir that the dynamic coefficient is proportional to the effective relative viscosity (µ_ew = µ_eff/µ_w) defined by them. Other studies did not consider those newly defined terms. Instead, the viscosity ratio (µ_n/µ_w) and wetting phase viscosity were often applied. Joekar Niasar, Hassanizadeh and Dahle [38] concluded the stronger dynamic capillary effects with a larger viscosity ratio. Later, Abbasi, et al. [174] numerically explored the viscosity dependency of dynamic capillary effect and found the dynamic coefficient increasing with increasing wetting phase viscosity, which is consistent with Equation (42) experimentally determined by Stauffer [37]. Moreover, Goel, Abidoye, Chahar and Das [170] investigated the viscous effects on dynamic relative permeability, in which they found the dynamic coefficient increase with mobility ratio (M = K_rwµ_n/K_rnµ_w) decrease.

The local heterogeneity influence: Mirzaei and Das [35] conducted a numerical study on micro-scale heterogeneities influencing the dynamic multiphase flow in porous media. In this study, different distributions and intensities of micro-scale heterogeneities were generated in the solving domain to study dynamic coefficients changed by those two influential factors. This study demonstrates that the dynamic coefficient is dependent on flow conditions and domain geometry. To be concise, the dynamic coefficient also increases with the higher intensity of heterogeneity. Other numerical studies from Helmig, et al. [175] and Abidoye and Das [169] also drew similar conclusions.

The permeability dependency: Manthey, Hassanizadeh and Helmig [40] and Mirzaei and Das [35] numerically and experimentally studied the transient two-phase flow in homogeneous and heterogeneous porous media at the continuum scale. A common research finding revealed amongst those works was that the dynamic coefficient is reversely proportional to the intrinsic permeability of porous media. This conclusion highly agrees with Stauffer [37] in Equation (42). Later, the same conclusion was confirmed again by Li, et al. [176] in dealing with rocks of the deep reservoir with different permeabilities. They also recommended the non-negligibility of dynamic effects for low permeable rock (K < 9.87 × 10⁻¹⁴ m²). Furthermore, the fractures as local heterogeneity in rocks also affect the dynamic coefficient. Salimi and Bruining [177] developed an upscaling model to study the dynamic effects in fractured porous media. It revealed that the dynamic effect enlarged with a higher seepage speed of fracture flow. There was low efficiency of oil recovery by water flooding for this scenario. Later, Tang, Lu, Zhan, Wenqjie and Ma [56] revisited the same research objective using numerical simulation, concluding the dynamic coefficient of fractured porous media is higher than that of unfractured porous media. Similar findings were repeatedly confirmed in later studies for fractured tight reservoirs [50,163,173].

The wettability dependency: O’Carroll, Mumford, Abriola and Gerhard [122] implemented two-phase multistep outflow experiments with inverse modelling to study the wettability dependency of dynamic effects. In addition to experimental exploration, O’Carroll, Mumford, Abriola and Gerhard [122] derived a microscale capillary advancing equation (Equation (44)). It is in the same form of macroscale dynamic capillary pressure equation from Hassanizadeh, Celia and Dahle [22], using Washburn [145] equation coupled with the current development of dynamic interfacial physics:

$$\frac{dl}{dt}={\left(\frac{2}{r}\zeta +\frac{8\mu L}{r^2}\right)}^{-1}\left(\Delta P+\frac{2{\sigma }_s\cos \mathsf{\theta }}{r}\right)$$

(44)

where l = advancing distance, t = advancing time, r = capillary tube radius, ζ = a coefficient of contact line friction, μ = fluid dynamic viscosity, L = total length of capillary tube, ΔP = two-phase pressure difference, σ_s = surface tension, and θ = static contact angle. By comparing Equation (44) with Equation (36), it is obvious to see that when the macroscale dynamic capillary equation is applied to a single capillary tube, the dynamic coefficient can be seen as

$$\tau =\frac{2}{r}\zeta +\frac{8\mu L}{r^2}$$

(45)

in which the calculation of the coefficient of contact friction (ζ) can be checked in detail from O’Carroll, Mumford, Abriola and Gerhard [122], with other notations in Equation (44). O’Carroll, Mumford, Abriola and Gerhard [122] and Li, Li, Chen, Guo, Wang and Luo [162] draw common conclusions that dynamic effect might be negligible for intermediate wetting conditions (60° < θ < 130°) but are necessary for more water-wet systems. In contrast, via numerically simulating flow in a bundle of capillary tubes, Mumford and O’Carroll [178] draw reversed conclusion that the larger the contact angle, the stronger dynamic effects.

The hysteretic dynamic effects: Mirzaei and Das [167] experimentally study the dynamic coefficient variation for primary drainage and main imbibition using 1D soil column experiments. Their result confirmed that the hysteresis nature of SWRC also happens to dynamic SWRC, and the dynamic coefficient is different for various hydraulic loading paths. However, only primary drainage and main imbibition data are available in this study. Thus, it provides the first vision to look into the hysteresis nature of dynamic SWRC. A similar experiment was also carried out by Sakaki, O’Carroll and Illangasekare [30] for one drainage-imbibition cycle. In this work, the theories of dynamic capillary pressure and fluid phase redistribution were both studied from their experimental data. They found the differences in dynamic coefficient (τ) and redistribution time (τ_B) between the drainage and imbibition process. Beforehand, Chen [179] already partially investigated hysteresis in dynamic effects for primary drainage, secondary drainage, and the main imbibition using a small flow cell setup, concluding the differences in a dynamic coefficient. Moreover, they found higher dynamic coefficients for the main imbibition than primary and secondary imbibition. The same trends also occurred for dynamic effects in relative permeability. Zhuang, Hassanizadeh, Qin and de Waal [47] experimentally revisited the hysteretic dynamic capillary effects. They further investigated scanning drainage curves beyond prior studies merely on primary drainage and main imbibition paths. Specifically, Zhuang, Hassanizadeh, Qin and de Waal [47] reconfirmed the log-linear relationship between dynamic coefficient and saturation and further found the nonuniqueness of dynamic coefficient for various hydraulic loading paths. Recently, Luo, Kong, Ji, Shen, Lu, Xin, Zhao, Li and Barry [161] revisited the hysteretic dynamic effects using a soil column test under fluctuated water table dynamics. This work experimentally determined the flow rate-dependent hysteretic dynamic coefficients, decreasing with an increasing fluctuation rate [161]. In summary, the cycling wetting and drying paths affecting dynamic coefficient is worth to be deeply studied in the future. There is also a need of identifying if the concept of the dynamic coefficient is a practical and straightforward approach for quantifying the discrepancy between static and dynamic SWRC with hysteresis.

The pressure boundary conditions effects: The pressure boundary conditions usually include continuously fluctuating hydraulic head, one-step or multistep in/outflow. Using an instrumented soil column test and full-scale embankment model, Scheuermann, Montenegro, et al. [180] investigated the so defined “hydraulic ratcheting” effect in soil water retention behaviour. Using the Spatial Time Domain Technique (Spatial TDR) and tensiometers, Scheuermann, et al. [180] found the matric suctions lost at minimal water contents for transient infiltration. Moreover, cumulative water content storage was first explored after several cyclic imbibition and drainage processes. This “ratcheting” effect was later verified by numerical studies using multiphase Lattice Boltzmann simulation [181,182]. However, other similar soil column setups from Cartwright, et al. [183] and Cartwright [184] presented no “ratcheting” effect but only hysteresis in static scanning curves. Those works also found no dynamic effects for both drainage and imbibition. The most recent experimental work from Luo, Kong, et al. [161] verified the dynamic effects for soil column tests under the fluctuated water table but had no confirmation of the “ratcheting” effect.

In addition to the water table dynamics, multistep in/outflow tests were carried by Schultze, Ippisch, Huwe and Durner [25], O’Carroll, Phelan and Abriola [28], Chen, et al. [185], Hui, Changfu, Huafeng, Erlin, Huan, Fratta, Puppala and Munhunthan [31], etc. There were a common in terms of multi-stepwise accumulative outflow in those studies. However, those works delivered conflicting manifestations of dynamic soil water retention curves. The results from O’Carroll, Phelan and Abriola [28] and Chen, Wei, Yi and Ma [185] showed no significant multi-stepwise overshot and undershot of soil suction for single soil water content on drainage and imbibition, respectively. In contrast, Schultze, Ippisch, Huwe and Durner [25] presented apparent multi-stepwise dynamic nonequilibrium effects. The conflicting results for multistep in/outflow should be deeply investigated in the future as well.

Acoustic excitation effects: Regardless of the mainstream reviewed above, some special boundary conditions raised research interests and are worth to be studied. One contribution from Lo, et al. [186] navigated the investigations in the dynamic response of soil water retention behaviour to a novel orientation. Instead of exclusively diving into studying hydraulic boundary conditions, Lo, Yang, Hsu, Chen, Yeh and Hilpert [186] applied acoustic excitations to transient drainage tests on Ottawa sand. Based on their experimental findings, Lo, Yang, Hsu, Chen, Yeh and Hilpert [186] concluded that the acoustic excitations insignificantly influence the static soil water retention curves but have more apparent effects on dynamic ones with the draining flow rate increases. Additionally, they experimentally determined that the dynamic coefficient decreases with frequency increases [186]. The hypothesised mechanism proposed by Lo, Yang, Hsu, Chen, Yeh and Hilpert [186] refers to the dynamic contact angle of capillarity varied with excitation frequency. Nonetheless, this work is merely an outstanding commencement of studying transient two-phase seepage under complex environmental conditions and still needs theoretical and experimental development in multiscale.

Summary: Although early studies and current experimental works offer many insights into dynamic SWRC on the macroscale, each study merely focuses on a few influential factors affecting the dynamic coefficient. Moreover, there are still many conflicts of research findings amongst those studies. Therefore, it is worthful to reinvestigate those conflicting conclusions in the future. On the other hand, Sakaki, O’Carroll and Illangasekare [30] suggested that future investigation should be focused on identifying extreme conditions on which the negligibility of the dynamic term can be determined. Moreover, the hysteretic behaviour of the dynamic coefficient (τ) or redistribution time factor (τ_B) needs further investigation. Finally, the influential factors for dynamic capillary effects still demand more experimental and numerical efforts to enhance the clarity and reinforce the principles concluded in prior works.

7.3. The Validations of Advanced Theories against Experiments

The investigation of transient two-phase flow in porous media focused more on the influential factors in dynamic capillarity effects because of the versatility of the thermodynamic-based theory. Specifically, there is only a single demand of dynamic coefficient to replicate the nonequilibrium transient two-phase flow seepage using the simplified version of the thermodynamic-based theory. However, as with other theories mentioned above, additional research efforts were devoted to validating all advanced theories. The validations were achieved by comparison against many experimental results lately published. The extra constitutive parameters (e.g., dynamic coefficient, fluid redistribution time, etc.) could be determined using straightforward experimental approaches or inverse modelling.

The validation of the Richards model: To the authors’ best knowledge, validating the Richards model could be sourced since the earlier 1960s. The invalidity of this theory was found in terms of dynamic nonequilibrium effects. Since Gardner [119] developed the method to calculate the diffusivity and hydraulic conductivity inversely using pressure plate outflow data, Rawlins and Gardner [123] later found the velocity-dependent diffusivity. In another straightforward interpretation, if one unsaturated soil hydraulic diffusivity in the function of moisture or soil suction is adopted, this theory will fail to simulate unsaturated soil water seepage for a different seepage velocity. Liakopoulos [24] also validated this model against a series of transient seepage tests using soil column setup and eventually concluded the failure as well. Nevertheless, instead of merely ending on the seepage speed dependence, Liakopoulos [24] further deducted the cause of failure was due to the inability to model soil water inertial using the two-phase Darcy seepage equation. Due to lacking a physical basis for modelling nonequilibrium scenario, the inversion of the Richards model somehow turned into a fitting process to determine the hydraulic properties. Therefore, new theories of fluid redistribution [151] and dynamic soil suction [37] were experimentally explored and mathematically developed for modern theoretical validations.

The validation of fluid redistribution theory: The aforementioned fluid redistribution theory in Section 6 was initially developed in the 1970s by Barenblatt [151] in petroleum engineering. Due to the complexity of integrating wetting phase fluid redistribution term into both capillary pressure or dimensionless one as Leverett J function-saturation and relative permeability-saturation, very few works on theoretical verification could be found in the literature, except the theoretical developer. Barenblatt and Vinnichenko [45] and Barenblatt, Garcia Azorero, De Pablo and Vazquez [172] continued to correct and modify this theory for modelling oil-water displacement in porous media. Finally, they successfully verified their model against the counter current seepage tests implemented by Zhou, et al. [187]. Later, Schembre and Kovscek [29] also applied this theory to simulate spontaneous countercurrent imbibition. The experimental results of spontaneous imbibition in diatomite (high porosity and ultra-low permeability) from Le Guen and Kovscek [188] and Zhou, et al. [189] were selected to compare with the simulation outcomes. Due to adding the fluid redistribution time, Schembre and Kovscek [29] could successfully simulate unsteady-state imbibition. Moreover, the shifting between nonequilibrium and static constitutive relationships could be eliminated by involving this term. However, this theory has been less validated in the latest decade because of overlooking dynamic capillary pressure from the physical perspective. The dynamic capillarity effects were more selected for modelling two-phase flow seepage through rock specimens even in petroleum engineering works. In fact, Sakaki, O’Carroll and Illangasekare [30] already proved the transferability between the fluid redistribution time and dynamic capillary coefficient.

As for soil hydrology, Ross and Smettem [152] validated their soil water redistribution model against the co-author’s experimental results. Later, this kinematic equilibration term was integrated into the dual-fraction model developed by Diamantopoulos, Iden and Durner [153]. Diamantopoulos, Iden and Durner [153] validated their newly proposed theory against multistep outflow experimental data and compared it with the Richards, single kinematic equilibration, and dual-porosity models. According to their simulating performance, the dual-fraction model coupled with the soil water redistribution model could simultaneously simulate the smooth transition of hydraulic pressure head and abrupt transition of accumulative outflow between each pair of in/outflow steps. In addition to multistep in/outflow, Diamantopoulos, Durner, Iden, Weller and Vogel [164] further validated their model for various boundary conditions published in prior experimental works [25,49,130]. Compared to the Richards model’s performance, their dual-fraction model succeeded in modelling the nonequilibrium transitional behaviour in terms of full state variables. Although there is a lack of nonequilibrium soil suction, it is still a successful and outstanding soil hydrological theory for modelling the transient seepage in unsaturated soil under nonequilibrium conditions.

The validation of dual-porosity and dual-permeability theory: Other theories such as dual-porosity and dual-permeability indeed can be applied to simulate the nonequilibrium effects in unsaturated soil water flow. The validations of those theories can be sourced from many. For instance, one outstanding contribution is from Šimůnek, et al. [159]. However, the physical nature of those theories is to resolve the transient seepage in macroscale heterogeneous porous media (structured porous media), which is beyond the scope of this review. For interests in transient seepage in structured porous media (e.g., soil with aggregate and clay content, rock with fissures and fractures), comprehensive reviews from Šimůnek, Jarvis, Genuchten and Gärdenäs [149] and Jarvis [190] may help to advance the understanding.

The validation of the thermodynamic-based theory: By far, the mainstream of modelling transient two-phase flow in porous media is based on the comprehensive theory developed by Hassanizadeh and Gray [43] and the simplified version given by Hassanizadeh, Celia and Dahle [22]. This thermodynamic-based theory applies to both soil hydrology and petroleum engineering.

As for soil hydrology application, Sander, et al. [191] numerically solved the simplified version in Equation (41) coupled with a hysteresis model in one- and two-dimensional domains. Thereby, the unstable fingering flow can be replicated by this numerical model. As for petroleum engineering, Cao and Pop [192] derived the uniqueness of weak solutions for the two-phase seepage theory coupled with both dynamic capillarity and hysteresis. El-Amin [193] developed a numerical solution of nonequilibrium two-phase seepage theory using the Implicit Pressure Explicit Saturation (IMPES) method. Meanwhile, Skiftestad [3] also developed a numerical solution using IMPES to investigate Microbial Enhanced Oil Recovery (MEOR). The details of this numerical model in terms of theory, numerical scheme, boundary and initial conditions, as well as numerical stability analysis, were provided entirely in this work. Even though there was no effort to validate the model with any experiments, they successfully conducted a parameter study using this numerical model. Their numerical investigation found the scale dependence of dynamic effects in capillarity and significant impact on interfacial tension by increasing the concentration of microbes [3]. Nevertheless, as this work had no validation against the oil recovery test, Skiftestad [3] set their prospects for collaboration with other experimental contributions. Hence, the conclusion based on the numerical solution might only bring inspiration and ideas for future experimental works on MEOR. Moreover, other numerical studies numerically solving this theory also include the investigation of local heterogeneity from Mirzaei and Das [35], the temperature influence in carbon dioxide sequestration from Das, et al. [194] and the scale effect from Hou, et al. [195], etc. However, none of the above validated their numerical model with experiment results but only numerically solved the theory with enormous mathematical efforts. The numerical simulations based on this theory have somehow been more utilised as numerical experiments with parameters studies.

There are several validations of numerical solutions against experiments from recent studies. For instance, Fučík and Mikyška [196] compared their numerical solutions calculated by a fully implicit scheme against a set of fast drainage and imbibition experimental results from the Colorado School of Mines. Distinctively, they concluded less importance of dynamic effects for homogeneous porous media but more dynamic effects for highly heterogeneous porous media. Moreover, their numerical solutions did not highly agree with the capillary pressure measured in experiments. Later, Das and Mirzaei [166] conducted a series of transient silicon oil-water seepage experiments in homogeneous coarse and fine sand columns to study the saturation dependence of the dynamic coefficient. Having conducted those tests, they also numerically simulated those experiments in 3D cylindrical domains using this theory for two-phase conditions. This numerical solver had been substantially developed by Mirzaei and Das [35] as well as Hanspal and Das [41]. Both experimental and numerical results mutually overlapped and also confirmed the nonlinear relationship between dynamic coefficient and saturation.

Afterwards, a petroleum engineering study from Zhang, et al. [197] reconfirmed the importance of considering dynamic effect into a core flooding test in sandstone. The numerical solutions of conventional and dynamic capillarity theory were implemented using IMPES as well and compared with the saturation profile detected by the X-ray CT technique. Zhang, He, Jiao, Luan, Mo and Guo [197] concluded on a better agreement with experiments when considering dynamic capillarity in the numerical model. Moreover, this dynamic term could yield relative permeability, which highly agreed with the relative permeability experimentally determined. Then, another petroleum study from Ren, et al. [198] numerically solved the theories of Buckley and Leverett [8], Barenblatt [151], and Hassanizadeh and Gray [43]. Based on a Bayesian analysis examining the efficacy of these theories to experimental results of decane/pentane-brine seeping in sandstone, Ren, Rafiee, Aryana and Younis [198] concluded that the two-phase fluid redistribution theory agreed more with experimental observations for the higher viscosity ratios of 4.4 and 15, while the other two theories outperformed the two-phase fluid redistribution theory for the lowest viscosity ratio of 1.1.

Besides, many soil hydrology studies also validate and develop thermodynamic-based theory. As the thermodynamic-based theory developed by Hassanizadeh and Gray [43] was based on a Gibbes free energy for air-water interfaces in an earlier theoretical work of Hassanizadeh and Gray [148], it can be used to model both dynamic effects by Equation (39) and hysteresis by the interfacial area advancing model. Therefore, Zhuang, et al. [199] compared the hysteresis-coupled Richards model and the thermodynamic-based theory against their experimental results to investigate horizontal soil water redistribution in terms of hysteresis. It finally concluded that the interface advancing model outperformed the conventional one with hysteresis. For validating these two models again, these horizontal unsaturated soil water seepage experiments were carried out with other settings by Zhuang, et al. [200]. They found the conventional model with hysteresis only performed well for the initial condition of very dry, while the interfacial area model can simulate for all experimental settings by fitting out the parameter in the interfacial production term in thermodynamic-based theory. Later, with the prior experimental findings on the log-linear relationship between saturation and the dynamic coefficient [47], Zhuang, van Duijn and Hassanizadeh [48] modelled saturation overshoot during the infiltration through dry sand in consideration of saturation-dependent dynamic coefficient. Not only was analytical analysis given in detail regarding parameter studies for saturation overshoot, but Zhuang, van Duijn and Hassanizadeh [48] also validated the thermodynamic-based theory with prior experimental works from DiCarlo [201] and Fritz [202] on primary imbibition exclusively. The most updated work from Zhuang, Hassanizadeh, et al. [203] modelled spontaneous imbibition in the thin fibrous layer. In this work, Zhuang, et al. [203] concluded that spontaneous imbibition is a fast capillary-driven process. Moreover, this process can only be reproduced by the dynamic capillarity theory but failed to be modelled by the conventional theory.

Summary: Many studies have validated all newly proposed theories against experimental outcomes. Due to more complicated constitutive relationships by adding additional dynamic terms, the numerical solution seems the only method to solve the theory mathematically. Most studies on validating each theory concluded that the newly proposed theory outperformed and outcompeted the conventional theory of two-phase flow in porous media and the Richards model. However, some numerical solutions did not match the experimental observations very well. Furthermore, the theory has been expanded to simulate versatility (e.g., various types of drainage, imbibition, hysteresis, etc.), whereas the fingering flow replication might be still debatable, referring to numerical stability. As a result, the numerical solutions of advanced theories are still worth to be applied to the model validations with laboratory experiments and specific geophysical scenarios with thermal, chemical, and mechanical coupling.

7.4. The State-of-Art of the Micromodels

With the development of imaging techniques and a transparent micromodel, there is another excellent chance to study dynamic SWRC in visualization. The terminology of the pore-scale physical model was defined as the micromodel shown in Figure 10. It has been applied to many disciplines to investigate two-phase flow in porous media. In detail, it was initially designed and fabricated in a few centimetres with small pore conduits (<1 mm). The in/outlets are opened on the left and right fluid reservoirs for in/outflow boundary control. The boundary condition can be regulated by pressure or flux controllers. The high transparency of the artificial model allows observing and recording seeping flow dynamics at the pore-scale level using visualization techniques (e.g., microscope and digital camera).

Those microscale mechanisms of pore-scale flow have been studied by Chen and Wilkinson [206], Lenormand, et al. [207], Avraam and Payatakes [208], Pyrak-Nolte, et al. [209], etc. A comprehensive review of the history of developing micromodel and fabricating methods could be sourced from Karadimitriou and Hassanizadeh [210]. Due to its advancements, it has been used to explore the invalidities of conventional theory for steady-state flow conditions. However, since the dynamic nonequilibrium effects under transient flow conditions drew more attention in recent two decades, it has been more applied to study transient two-phase flow seepage with advanced theories in soil science and petroleum engineering.

Earlier studies of transient two-phase seepage using micromodel: Most earlier works on micromodel focused on steady-state flow, which is not the objective of this review. To the best of the authors’ knowledge, Tsakiroglou, Theodoropoulou, et al. [124] presented the earliest study on nonequilibrium effects using micromodel, followed by another continuous work on transient and steady-state flow [204]. Tsakiroglou, Theodoropoulou, et al. [124] carried out drainage experiments by percolating paraffin oil-water flow through a hydrophilic planar glass-etched pore network in their earlier work. They also conducted numerical simulations solving thermodynamic-based two-phase flow theory. The parameters in the continuum-scale model were determined by inversion analysis with a Bayesian estimator. In addition to the continuum-scale study, they also explored the microscopic percolation theory on these experiments. As a result, they determined the relation between macroscopic and microscopic parameters in percolation theory: dynamic coefficient to the capillary number and other scaling coefficients, etc. In this study, Tsakiroglou, Theodoropoulou and Karoutsos [124] concluded on the dynamic effects in capillarity and relative permeability depending on the capillary number (10⁻⁵–10⁻⁸), flow regime (capillary or viscous dominating transient flow), flow pattern (capillary or viscous finger), pore network size, etc. In addition, this flow pattern altered from more clustered patterns to shaper driving front with the capillary number increasing. Moreover, Tsakiroglou, Theodoropoulou and Karoutsos [124] determined the decreasing relation between dynamic coefficient and capillary number and finally emphasised on applying flow rate-dependent constitutive relationships for simulating transient two-phase seepage flow dominated by viscous effects.

The same conclusion on the relation between dynamic coefficient and the capillary number was reconfirmed again by a later work from Tsakiroglou, Avraam and Payatakes [204]. This continuous work further explored the wettability-dependence on the relative permeability-saturation relationship under transient flow conditions. Specifically, as for less hydrophilicity, the transient relative permeability exceeds the corresponding steady-state one for smaller capillary numbers, while tending to be smaller than the steady-state one for higher capillary numbers [204]. As for better hydrophilicity, the transient relative permeability increases with capillary number and the transient nonwetting phase relative permeability is higher than the steady-state one, yet the transient wetting phase relative permeability substantially decreases and turns to be smaller than the steady-state one for smaller capillary numbers [204]. This series of experiments shed light on dynamic nonequilibrium effects under transient two-phase flow seepage. However, it still merely investigated the flow pattern and regimes in terms of capillary number as well as wettability. The interfaces dynamics or specific interfacial areas were not investigated.

Recent studies of transient two-phase seepage using micromodel: One recently representative work on multiphase flow in a transparent micromodel was given by Karadimitriou, Hassanizadeh, Joekar-Niasar and Kleingeld [51]. Their work originally started from Karadimitriou and Hassanizadeh [210], in which a review of micromodels for two-phase flow studies is given. In this work, fabrication methods, including the Hele-Shaw cell, Lithography, and Wet-etching technique, as well as visualizing processes, such as the camera, camera coupling with a microscope, and laser-induced fluorescence, were reviewed to show their advantages and drawbacks. Karadimitriou, Joekar-Niasar, et al. [211] started with a glass-etched micromodel for a two-phase flow experiment based on this experimental method review. Their experimental result agreed with numerical pore network simulation developed by Joekar-Niasar and Hassanizadeh [150] for a single SWRC.

However, Karadimitriou, et al. [211] pointed out the disadvantages of deep reactive ion etched (DIRE) fabricated micromodel in three points: anisotropic of the DRIE model, sloping vertical wall, and rough boundary surfaces. Hence, Karadimitriou, Musterd, Kleingeld, Kreutzer, Hassanizadeh and Joekar-Niasar [205] used a soft lithography method to fabricate polydimethylsiloxane (PDMS) micromodel for a two-phase flow experiment. However, the PDMS micromodel has one critical drawback. The wettability of material changes with time when two-phase fluids are mixed into the porous matrix [210]. Therefore, they developed a reliable procedure for PDMS fabrication with a relatively stable wettability in this study [205]. Furthermore, this reliable micromodel experimentally proved the uniqueness of the capillary pressure-saturation-specific interfacial area (P_c-S_w-A_nw) at the static condition.

The most updated study using this physical model for the isothermal two-phase flow was presented by Karadimitriou, Hassanizadeh, Joekar-Niasar and Kleingeld [51]. It is a contributive study of two-phase flow under transient flow conditions using the PDMS micromodel. In this study, they physically demonstrated that the proposed constitutive surface was unique for a given capillary number. Moreover, one constitutive surface cannot be used for both transient and steady-state flow conditions [51]. By far, this PDMS micromodel was continuously implemented to compare with other numerical studies by Kunz, et al. [212], Nuske, et al. [213], Konangi, et al. [214], and Yiotis, et al. [215].

7.5. The State-of-Art of the Pore-Scale Simulations

Computational Fluid Dynamic (CFD) simulation also contributes to experimental studies on two-phase flow in porous media. To the authors’ best knowledge, there are mainly three numerical experiments: Pore Network Model (PNM), Direct Numerical Simulations (DNS) and Lattice Boltzmann Method (LBM). Each model has been well-developed, so these methods can be directly used to explore constitutive relationships under transient conditions by applying them in the micro-sized domain where pore conduits can be visually identified.

The Pore Network Model (PNM) application: The representative work on PNM of two-phase flow and validating thermodynamic basis theory using PNM were conducted by Joekar-Niasar and Hassanizadeh [150]. The PNM has two types: PNM for static conditions and the Dynamic Pore Network Model. The pore matrix can be artificially generated into different forms, such as a pure pore throat network or tubes coupled with large pores to study non-wetting phase trapping mechanisms. In an early study of PNM from Joekar-Niasar, et al. [216], the PNM was generated using the sphere-tube network, and the Poiseuille Law calculated the flow rate of each pore throat. After upscaling by averaging variables on a tube scale, they found the uniqueness of P_c-S_w-A_nw surface and k_rw-S_w-A_nw surface under static conditions. So this work confirmed that coupling specific interfacial area with the original constitutive curve to form a 3D parabolic surface can reduce the hysteresis between boundary curves (primary drying and wetting) [216].

Later, Joekar Niasar, et al. [217] improved PNM to own unstructured pore network features for studying SWRC and HCF under static drainage and imbibition. A significant finding is that non-wetting phase trapping will change the parabolic relationship between specific interfacial area and wetting phase saturation. In another way, the unique parabolic P_c-S_w-A_nw surface can only be generated from PNM when there is no tapping in PNM [217].

Since 2010, this PNM was upgraded to Dynamic Pore Network Model (DPNM) using Washburn [145] dynamic capillary flow equation for a single tube [38]. As a previous content review of dynamic interfacial physics, the Lucas–Washburn equation can only be applied for the laminar flow regime where capillary driving is balanced by viscosity. In this DPNM, the dynamic interfacial physics (dynamic interface or any free energy terms) was found never to be used in the determination of the driving force term. Instead, Joekar Niasar, Hassanizadeh and Dahle [38] geometrically derived a pore-scale capillary pressure-saturation relationship to provide tube-scaled dynamic capillary pressure in a Lucas–Washburn equation, still based on the static Young–Laplace law. This raises concern on how a pore network can generate a dynamic capillary flow using static capillary geometry on this pore scale. However, according to the causes of dynamic effects reviewed above, other causes also include local heterogeneity. Thus, it still could be used to validate the advanced theory proposed by Hassanizadeh and Gray [43]. It confirmed this theory that derivatives of interfacial area and saturation should appear in driving terms of the new Darcian capillary flow motion equation (Equation (38)) [218].

Joekar-Niasar and Hassanizadeh [150] studied dynamic capillary flow in angular pore network using DPNM and gave the same conclusion as Karadimitriou, Hassanizadeh, Joekar-Niasar and Kleingeld [51] that proposed constitutive surface (P_c-S_w-A_nw) exists for non-equilibrium drainage and imbibition but the difference between dynamic surface and static surface increases with both viscous forces and capillary number increase. Later, Sweijen, et al. [219] upgraded the angular DPNM to the tetrahedra DPNM using the Discrete Element Method (DEM). Then, they carried out the DPNM simulation through the genuine tetrahedra pore unit formed by an assembly of granular particles. The same as Joekar-Niasar and Hassanizadeh [150], Sweijen, Hassanizadeh, Chareyre and Zhuang [219] also applied the static capillarity into pore-scale modelling and more focused on generating macroscale dynamic effects by regional heterogeneity. Finally, they concluded that regional heterogeneity-induced fingering flow formed dynamic capillarity effects in macroscale.

A comprehensive literature review of PNM and DPNM is available from Joekar-Niasar and Hassanizadeh [220], in which the authors of DNPM already acknowledged neglecting dynamic contact angle and inertia effects. However, DPNM is still useful for analysing dynamic two-phase flow in porous media for a larger continuum scale (several REV in simulating domain). The DPNM costs less computational expenses than later introduced pore-scaled numerical simulation methods (DNS and LBM), which are only conducted on the microscale equal to or less than a single REV. Joekar-Niasar and Hassanizadeh [220] also confirmed a challenge of validating DPNM using currently available experiments. Compared DNPM with the other two methods, the governing equation for multiphase flow in DPNM is still Poiseuille-based. Therefore, the driving forces should be clearly identified. Otherwise, it is impossible to replicate dynamic two-phase flow completely without additional terms fully capturing dynamic capillary pressure at pore-scale.

The Lattice Boltzmann Method (LBM) application: Pore-scale simulation is a method in which two-phase fluids flow are simulated by solving either the Boltzmann equation (LBM) or the Navier–Stokes equation (DNS) in a domain with artificially packed beads. Compared to DNPM, both LBM and DNS are computationally expensive and applied to a small domain. The Boltzmann equation is a statistical equation used to describe the streaming and collision of fluid molecules in space. The Boltzmann equation is more fundamental and focuses on particle statistics than the Navier–Stokes equation describing momentum conservation with external force action for a single Control Volume (CV). The Boltzmann equation is

$$f^1\left(\mathbf{x}+d\mathbf{x},\mathbf{p}+d\mathbf{p},t+dt\right)d\mathbf{x}d\mathbf{p}=f^1(\mathbf{x},\mathbf{p},t)d\mathbf{x}d\mathbf{p}+\left[{\Gamma }^{+}-{\Gamma }^{-}\right]d\mathbf{x}d\mathbf{p}dt$$

(46)

where f¹(x,p,t) is the probability of finding one molecule with given position x and momentum p in time t, Γ⁺ is molecules invading into aiming portion of molecules, and Γ⁻ are molecules escaping out of a targeting region of molecules [221]. The last term of the right-hand side describes the particle collision, and the other two terms form particle streaming. This nonlinear partial differential equation has no analytical solution, so the continuum equation is discretised into lattices. Thus, this discretised Boltzmann equation can be solved numerically, thereby defined as the Lattice Boltzmann Method (LBM). Sukop [222] interpreted the Boltzmann equation and solutions for single and two-phase flow in detail. Initially, solving the discretised Boltzmann equation could only simulate the single-phase fluid in a domain. With the addition of three forces: an adhesive force between wetting phase particles and solid surface, a repulsive force between non-wetting phase particles and solid surface, and another repulsive force between two immiscible fluid phases’ particles, the interface can be generated between two fluids in LBM. Shan and Chen [223] invented this method of two-phase generation, so it is also named Shan–Chen Lattice Boltzmann Method (SCLBM).

Shan–Chen LBM has already demonstrated its powerful function in studying steady and unsteady two-phase flow in porous media. Sukop and Or [224] applied SCLBM to study fundamental multiphase flow physics, including surface tension determination, adsorption, capillary condensation, and capillary pressure-saturation of a complex pore network. This work discussed the limitation of SCLBM due to its thermodynamic basis (discrepancy between ideal Equation of State (EOS) and LBM EOS). Moreover, it gave an outlook on studying transient capillary flow using such a powerful tool. Pan, et al. [225] approximated the solution of the Boltzmann equation using SCLBM in a spherical package to determine SWRC and hysteresis curves. LBM results agreed with the experimental result generated from a two-phase displacement in the glass beads package. Ahrenholz, et al. [226] further studied hysteresis of SWRC using SCLBM. They compared PNM and DPNM to demonstrate the advantage of SCLBM, which allows for observing microscale phenomena. Regarding the P_c-S_w-A_nw constitutive relationship proposed by Hassanizadeh and Gray [148], Porter, et al. [227] applied SCLBM in a 2D domain and validated SCLBM against experimental results from a computed microtomography (CMT). Their LBM results did not only show good agreement with experiments, but they also confirmed the non-hysteretic constitutive surface. However, the uniqueness of the P_c-S_w-A_nw surface was later challenged by Galindo-Torres, et al. [228]. They demonstrated the non-uniqueness of the parabolic constitutive surface using SCLBM by applying two-phase displacement in different principal directions. Most of the recent SCLBM studies focus on experimental validating SCLBM and studying impact factors (pore geometry, viscosity of fluid, wettability, fracture features, etc.) influencing SWRC and HCF under steady-state flow conditions [229,230,231,232,233]. Only very few looked into the transient two-phase flow in porous media.

The SCLBM study investigating transient effects during wetting-drying cycles in unsaturated soil was given by Galindo-Torres, Scheuermann, Li, Pedroso and Williams [182]. Compared to other SCLBM studies neglecting gravity for less computational effort, this work considered gravity to define the hydraulic head by a linear fitting relationship between wetting fluid density and hydraulic head. A sinusoidal oscillation water table was generated to investigate the hydraulic ratcheting effect (moisture content accumulation under cycling hydraulic loading) by SCLBM. However, this study did not give any SWRC and HCF data under transient conditions but offered a relationship between oscillating hydraulic head and saturation vibration. Liu, et al. [234] and Liu, et al. [235] simulated a process of carbon dioxide (CO₂) sequestration in porous media using colour-fluid LBM and quantified the transient effect by Capillary number (Ca). This Ca-dependent behaviour agreed with the micromodel experimental result on constitutive surface varying with Ca from Karadimitriou, Hassanizadeh, Joekar-Niasar and Kleingeld [51]. However, they also did not provide any information on SWRC, HCF, and P_c-S_w-A_nw surface. Instead, they determine a unique relationship between non-wetting phase saturation and specific interfacial length (specific interfacial area in the 2D domain) under different Ca. So far, this new proposed relationship has not been validated by any experiments yet but a relationship given by colour-fluid LBM. Other dynamic multiphase flow studies using SCLBM only investigated capillary tubes rather than packing beads [236,237]. As the object was only a capillary tube, they only studied the dynamic contact angle varying with Ca. Therefore, there are few LBM studies on the transient effect of SWRC and HCF, whereas LBM can reproduce the dynamic capillary pressure using Shan and Chen [223] inter-particle potential method.

The Direct Numerical Simulation (DNS) application: The newest method for pore-scaled two-phase simulation is achieved by numerically solving the Navier–Stokes (N-S) equation with interfacial tracking by Volume of Fluid (VOF). This method is also defined as DNS. One representative work was given by Ferrari and Lunati [238]. The significance of this study is that a force resulting from surface tension is accounted into a body force portion of the N-S equation. In addition, this surface tension force is not determined by the static Young–Laplace equation but calculated from a dynamic capillary equation having Helmholtz free energy, which transforms back to the Young–Laplace equation when the infinitesimal increment of Helmholtz free energy is minimised to zero (following equilibrium approached). This utilization agrees with a thermodynamic basis theory developed by Hassanizadeh and Gray [148]. Compared to LBM using a virtual lattice unit and needing case-specific calibration, DNS coupling with fluid function accounting Helmholtz free energy can directly use physical parameters to simulate dynamic capillary pressure and inertia effect [238]. This requires less effort for model parameters calibration. Using this DNS model, Ferrari and Lunati [238] demonstrated the invalidity of Darcy law due to neglecting viscous effects and trapping.

Latter, Ferrari and Lunati [52] studied the inertial effect of the meniscus and demonstrated that different phase distributions could be achieved by varying inertial force during the fluid redistribution process. Despite the emphasis on considering non-equilibrium effects, Ferrari and Lunati [52] also concluded on which conditions the dynamic effects should be considered:

• if the characteristic time of the applied boundary flow is larger than redistribution time, Darcy’s law is still applicable;
• otherwise, dynamic should be considered.

This finding provides an approach to distinguish the conventional multiphase flow theory from advanced theory considering dynamics by comparing boundary flow conditions with fluid-phase redistribution time. Ferrari, et al. [239] validated this DNS model against experimental results from a Hele-Shaw cell packing cylindrical obstacles. Although uncertainties of spatial features induced some mismatches between experiments and simulations, good agreements could be achieved between them regarding upscaled state variables [239]. However, although Ferrari [240] developed this method, there has been no study on SWRC and HCF. Therefore, it is necessary to replicate this method using accessible commercial multiphase simulation software. This might offer another opportunity to investigate the dynamic effect in SWRC and HCF using numerical experiments.

Despite PNM, LBM, and DNS (VoF), other fluid dynamic simulation methods still exist for modelling transient two-phase flow in porous media. It includes another DNS with Level-Set (LS) and Smoothed Particle Hydrodynamics (SPH) methods. One representative work from Helland, et al. [241] applied DNS with the LS method to explore the uniqueness of capillary pressure-saturation relationship under the steady-state flow condition. However, they made no attempts to investigate transient two-phase flow conditions, which could be a future research interest for LS simulation. In contrast, a few contributions from Kunz, Zarikos, Karadimitriou, Huber, Nieken and Hassanizadeh [212], Sivanesapillai, Falkner, Hartmaier and Steeb [53], and Sivanesapillai and Steeb [242] indeed studied the dynamic interfaces of immiscible two-phase flow in porous media using SPH. In addition, a method named continuum surface force (CSF) was introduced into the SPH model to simulate the two-phase separation. Sivanesapillai, Falkner, Hartmaier and Steeb [53] comprehensively investigated the pore-scale phenomenon and REV-scale outputs with various fluid properties settings using the SPH-CSF method. For interest in conducting numerical experiments by this method, their works bring more insights in modelling development and multiphase physics for transient two-phase flow dominated by capillary-viscous effects.

8. Engineering Applications in Unsaturated Soil Mechanics

In comparison to soil water equilibration theory, the simplified thermodynamic-based theory could bring an additional chance to simulate soil suction loss-induced shear strength reduction during intensive rainfall. Moreover, as previous rainfall-trigged landslides modelled by the Richards model coupled with slope stability analysis (Ng and Shi 1998), the development of dynamic nonequilibrium theory contributes to expanding the predicting model of rainfall-induced landslides with higher temporal precision.

8.1. Observation of Transient Effects in Natural Slopes

Although several studies did not find dynamic nonequilibrium or transient effects under transient two-phase flow conditions [183,184,243,244,245], the literature cited in the previous content has already proven the existence in laboratory experiments. Durner, et al. [246] also partially reviewed and discussed dynamic/transient effects and concluded ambiguity on whether it is critical in the field because other conditions such as heterogeneity might manifest more significantly. Moreover, the earlier study on rainfall-induced slope stability from Ng, et al. [247] has not reported any observation of transient effects. However, it is still imprudent to conclude there are no such effects in situ because the observing scale is enlarged to the large scale for the field. Observing transient effects depends on rainfall intensity, initial soil water content, groundwater dynamics, in situ intrinsic permeability, functioning mechanism of sensors, temperature, humidity, etc.

The latest fieldwork of Soil Water Retention Curve (SWRC), which proved the existence of transient effects on imbibition SWRC in the field, was presented by Bordoni, Bittelli, Valentino, Chersich and Meisina [66]. Their work implemented in situ soil water content and suction tests using Time Domain Reflectometry (TDR) probes and tensiometers in Montuè and Centonara test sites in Italy. Bordoni, Bittelli, Valentino, Chersich and Meisina [66] found the nonequilibrium process due to fast infiltration during intense summer rainstorms (10–25 mm/h), in which soil suction was almost totally zero for a nearly constant soil water volumetric content of 20% at 0.4 and 0.6 m from the ground surface. Furthermore, the transient effects smeared out with a depth of 1.2 m under the ground surface. The soil, in which they investigated transient effects, has less content of sand and gravel (<20–30%) but high content of silt and clay in 50–60% and 20–30% respectively. It differs from the uniform sandy soil high frequently selected in laboratory experiments.

This in-situ test is extremely rare in demonstrating the dynamic nonequilibrium effects for the transient infiltration process. Therefore, it will foster experimental and numerical efforts on this research objective as the prospects for developing unsaturated soil mechanics. However, there are few applications of nonequilibrium transient two-phase seepage theories in geotechnical engineering. Only several attempts in Geotechnics are available so far and, therefore, are reviewed in this section for a modest spur to induce someone to come forward with more valuable contributions.

8.2. Transient Effects Coupled in Unsaturated Soil Effective Stress

Unsaturated soil is a unique material for which soil suction needs to be included, compared to saturated soil, where positive water pressure is the only neutral stress changing soil skeleton stress. Bishop [248] already gave the stress state variable for unsaturated soil, which is

$$\sigma \prime =(\sigma -u_a)+\chi (u_a-u_w)$$

(47)

where σ′ = effective stress of unsaturated soil, σ = total normal stress, u_a = pore air pressure, u_w = pore water pressure, and χ = Bishop factor depends on effective soil saturation (S_e). In most cases, pore air pressure is assumed at atmospheric pressure as a zero reference pressure so that the pore water pressure can be defined as the soil suction when it turns to be negative pore water pressure. Equation (47) has been experimentally validated by Fredlund and Morgenstern [59] and also derived in the textbook of Lu and Likos [1] with spherical packing assumption. Although there are some other forms of unsaturated soil effective stress considering solute suction (osmotic suction) [249,250], the Bishop effective stress form is still the most generic equation used in Geotechnics.

Based on the thermodynamic-based theory of transient air-water flow in unsaturated soil, Nikooee, et al. [251] further extrapolated the effective stress to

$$\begin{gathered} \sigma \prime =\sigma -u_a+{\sigma }_s=\sigma -u_a+\chi {\psi }_m+k^{nw}{a}^{nw} \\ {k}^{nw}=-\frac{\partial A^{nw}}{\partial n}{\Gamma }^{nw} \end{gathered}$$

(48)

where σ_s = suction stress proposed by Lu and Likos [60], ψ_m = static soil matric suction (equilibrium soil suction), k_nw = a material constant, a = specific immiscible fluids interfacial area, A^nw = Helmholtz free energy of interfaces, n = porosity, Γ^nw = interfacial mass density and others as same as notations in Equation (47).

With Equation (48), the effective stress of unsaturated soil under static or equilibrium conditions can be described for all hydraulic loading paths, notably including hysteresis. Moreover, after the parabolic surface was fitted into the constitutive relationship of soil suction-saturation-specific interfacial area, the suction stress (σ_s) can be determined, subsequently calculating the effective stress (σ′) as well. In this way, the conventional hysteretic Soil Water Retention Functions (SWRF) will not be required to give the Bishop χ factor. Nikooee, Habibagahi, Hassanizadeh and Ghahramani [251] successfully validated this derivation against the Soil Water Retention Curves (SWRC) of kaolin, silt, and mixtures of sand and clay samples. For performance in terms of suction stress-soil suction-saturation constitutive surfaces, the work from Nikooee, Habibagahi, Hassanizadeh and Ghahramani [251] deserves to be followed up on.

Later, Nikooee, et al. [252] additionally derived the effective stress considering dynamic nonequilibrium effects in the mathematical form:

$$\sigma \prime =\sigma -u_a+{\sigma }_s=\sigma -u_a+\chi ({\psi }_m-\tau \frac{\partial S}{\partial t})$$

(49)

where τ = dynamic coefficient, S = saturation, t = time, and the rest notations were already given in Equations (47) and (48). The formulation of Equation (49), when χ = S_e = (S − S_r)/(1 − S_r), was accepted for sandy soil in other hydro-mechanical coupling studies [253,254].

Nevertheless, Nikooee, Hassanizadeh and Habibagahi [252] made no efforts to experimentally validate this proposed unsaturated soil effective stress under transient flow conditions. Since validating effective stress usually demands direct shear coupled with the soil water retention test, there were no available experiments to test shear strength for transient flow conditions at that moment.

8.3. Transient Effects Coupled in Unsaturated Soil Shear Strength

To the authors’ best knowledge, the shear strength determination of unsaturated soil under transient flow conditions was investigated by Milatz, Törzs, Nikooee, Hassanizadeh and Grabe [255]. However, before they achieved the goal, a series of experiments were developed in the following steps.

Initially, Milatz and Grabe [254] set up a simple shear testing device to study coarse-grained soils’ mechanical behaviour under monotonous and cyclic mechanical loading conditions. This shear apparatus was suction-controlled using the conventional Axis Translation Technique (ATT), desaturating and saturating soil water in/out of soil specimen in an air-opened flow cell by vacuum controller. It was similar to the hanging column setup but set in a smaller flow cell. Instead of determining suction by suction control, a tensiometer centrally embedded into the porous cap disk could measure the soil suction after sensor insertion by capping the disk. Its advantage is high accuracy in the low suction range <10 kPa. However, due to the tensiometer’s measuring upper constraint, this apparatus could only detect soil suction <85 kPa. Thus, it can only be applied to coarse-grained soil above silt and clay. In addition, this direct shear test was still achieved as the standard direct shear but with a circular shear box instead of the square one. Furthermore, they succeeded in carrying out drained-constant suction and undrained-constant water content tests. Finally, Milatz and Grabe [254] concluded on hysteresis-induced dilatancy as well as contracting and suction-reinforcing soil specimens during cyclic shearing.

Later, Milatz, et al. [256] investigated the settlement behaviour of unsaturated granular soils by varying soil suction and water content. Using the same device from Milatz and Grabe [254], the shear loading was switched off, and the normal loading was preserved so that the suction-controlled oedometer could be achieved. Eventually, according to their experimental outcomes, Milatz, Törzs, et al. [256] replicated the wetting-induced soil collapse. This shearing and consolidating apparatus demonstrated its advantages in studying the hydro-mechanical coupling behaviour of unsaturated coarse-grain soil. However, this apparatus could not investigate transient flow conditions but only equilibrium conditions because of laying a ceramic porous disk under the soil specimen.

Milatz, Törzs and Grabe [253] upgraded their experiment by replacing the ceramic disk with a microporous filter membrane, preventing air breakthrough in the testing system while preserving the high permeability. This upgraded setup succeeded in producing dynamic effects in SWRC on primary drainage and main imbibition. Furthermore, by applying the suction stress (Equation (49)) concept proposed by Lu and Likos [60], Milatz, Törzs and Grabe [253] successfully estimated soil suction stress as

$${\sigma }_s=\chi ({\psi }_m-\tau \frac{\partial S}{\partial t})=S_e({\psi }_m-\tau \frac{\partial S}{\partial t})$$

(50)

in shear stress of unsaturated soil:

$$\tau =c\prime +(\sigma -u_a)\tan \varphi \prime +{\sigma }_s\tan \varphi \prime =c\prime +(\sigma -u_a)\tan \varphi \prime +\chi (u_a-u_w)\tan \varphi \prime$$

(51)

where τ = the shear stress, c′ = cohesion at zero suction, Φ′ = the friction angle for zero suction, and other notations are the same as above. Their work reconfirmed the nonequilibrium soil suction overshot equilibrium on primary drainage. Controversially, their nonequilibrium soil suction also overshot equilibrium one on the main imbibition. In fact, prior studies showed soil suction undershooting for the imbibition process. As a result, both Suction Stress Characteristic Curve (SSCC) for transient primary drainage and imbibition overshot the SSCC under low seepage velocities. Later, Milatz, Törzs, Nikooee, Hassanizadeh and Grabe [255] continued this experimental and theoretical work by expanding the hydraulic loading path to involve transient hysteretic scanning loops. Similar to the last work of Milatz, Törzs and Grabe [253], Milatz, Törzs, Nikooee, Hassanizadeh and Grabe [255] also found dynamic soil suction overshooting phenomenon for drainage, imbibition, and hysteretic scanning curves. According to those findings, Milatz, Törzs, Nikooee, Hassanizadeh and Grabe [255] finally estimated dynamic SSCC, which overshot equilibrium SSCC, even for dynamic hysteretic SSCC partially.

Compared to the theoretical work from Nikooee, Hassanizadeh and Habibagahi [252], Milatz, Törzs, Nikooee, Hassanizadeh and Grabe [255] indeed carried out experimental investigations on transient air-water seepage. They also estimated soil suction stress using available SSCC theory [60] coupled with dynamic capillarity theory [148]. It is a contributive commencement for studying the mechanical behaviour of unsaturated soil under transient flow conditions. Nevertheless, the soil suction contribution to the shear strength as suction stress is still a theoretical prediction without validation by shear tests. So far, this validating work is still open to those who have great interests and intend to tackle those experimental challenges.

8.4. Transient Effects Coupled in Unsaturated Soil Deformation

Due to the advancements and complexities in the thermodynamic-based theory, effective stress and shear strength characterization already have numerous difficulties theoretically and experimentally. Therefore, it is even more challenging for any attempt in deformation modelling. To the best of the authors’ knowledge, only one theoretical formulation with numerical solution from Zou, et al. [257] is available in the literature by far. Their work considered the dynamic capillarity effects in simulating transient two-phase flow in porous media. Meanwhile, they also modelled the drainage-induced settlement (shrinkage) by expanding the pore-elastic consolidation theory [258,259,260] to embrace transient two-phase flow simulators.

Table 6. The entire theoretical framework for modelling unsaturated soil elastic deformation in coupled with the thermodynamic-based theory of transient soil gas-water flow in unsaturated soil.

Seepage Equations	Mathematical Formulations
Governing equation for transient soil water flow in unsaturated soil considering soil skeleton deformation (Zou, Saad and Grondin [257])	$\begin{array}{l}\frac{\partial {\rho }_w{n}_e{S}_w}{\partial t}+\nabla \left[{\rho }_w\left(\frac{K_{r}K}{{\mu }_w}(\nabla P_c{}^{dyn}-{\rho }_{w}g)+n_e{S}_w\frac{\partial u}{\partial t}\right)\right]=0\\ \frac{\partial {\rho }_w{n}_e{S}_w}{\partial t}={\rho }_w{n}_e\frac{\partial S_w}{\partial P_c{}^{stat}}\frac{\partial P_c{}^{stat}}{\partial t}+{\rho }_w{S}_w\frac{\partial n_e}{\partial t}\end{array}$	(52)
Dynamic capillary pressure (Hassanizadeh and Gray [148])	$P_c{}^{dyn}-P_c{}^{stat}=(P_a-P_w)-P_c{}^{stat}=-\tau \frac{\partial S_w}{\partial t}{P}_a=0$	(53)
Effective porosity-volumetric strain coupling (Zou, Saad and Grondin [257])	$n_e=\frac{n_0+{\epsilon }_v}{1+{\epsilon }_v}$	(54)
Permeability-effective porosity coupling (Modified Kozeny-Carman Wang and Nackenhorst [260])	$K=K_0\frac{n_e{}^3}{{(1-n_e)}^2}\frac{{(1-n_0)}^2}{n_0{}^3}$	(55)
Dynamic coefficient-porosity and permeability coupling (Stauffer [37])	$\tau =\frac{{\alpha }_{stauffer}{n}_e{\mu }_w}{K{n}_{BC}}{\left(\frac{{\alpha }_{BC}}{{\rho }_{w}g}\right)}^2$	(56)
Soil water retenion fuction(Modified van Genuchten [14])	$S_e=\frac{S_w-S_r}{1-S_r}={\left[\frac{1}{1+{\alpha }_{VG}{(P_c{}^{stat})}^{n_{VG}}}\right]}^{m_{VG}}{m}_{VG}=1-\frac{1}{n_{VG}}$	(57)
Relative permeability-saturation model (Modified van Genuchten [14])	$K_r=S_e{}^{0.5}[1-{(1-S_e{}^{1/m_{VG}})}^{m_{VG}}]$	(58)
Brooks & Corey to modified van Genuchten fitting parameters transformation (Morel-Seytoux, et al. [261])	$\begin{array}{l}\frac{{\alpha }_{BC}}{{\rho }_{w}g}=\frac{1}{{\alpha }_{stauffer}}\frac{p+3}{2p(p-1)}\frac{127.8+8.1p+0.0092p^2}{55.6+7.4p+p^2}\\ p=3+\frac{2}{n_{VG}}\end{array}$	(59)
Mechanical equations	Mathematical formulations
Unsaturated soil effective stress (Biot and Willis [259])	$\begin{array}{l}{\sigma }_{ij}\prime ={\sigma }_{ij}-b\left[(1-S_w)P_a+S_w{P}_w\right]{\delta }_{ij}\\ P_a=0; b=1-\frac{{\kappa }_b}{{\kappa }_s}\le 1\end{array}$	(60)
The full pore-elastic stress-strain constitutive relationship (Biot [258])	$\begin{array}{l}{\sigma }_{ij}\prime ={\lambda }_b{\epsilon }_v{\delta }_{ij}+2G_b{\epsilon }_{ij}\\ {\epsilon }_v={\epsilon }_{xx}+{\epsilon }_{yy}+{\epsilon }_{zz}\\ {\epsilon }_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\\ {\lambda }_b=\frac{v_b{E}_b}{(1+v_b)(1-2v_b)},G_b=\frac{E_b}{2(1+v_b)}\end{array}$	(61)
Simplified pore-elastic constitutive relationship under isotropic loading condition (Zou, Saad and Grondin [257])	$\begin{array}{l}{\epsilon }_v=\frac{{\sigma }_{mean}}{{\kappa }_b}+S_w{P}_c{}^{dyn}\left(\frac{1}{{\kappa }_b}-\frac{1}{{\kappa }_s}\right)\\ {\kappa }_b=\frac{(3{\lambda }_b+2G_b)}{3}\\ {\sigma }_{mean}=\frac{{\sigma }_{xx}+{\sigma }_{yy}+{\sigma }_{zz}}{3}\end{array}$	(62)
Nonequilibrium soil suction increase-induced shrinkage by assuming rigid particles(Zou, Saad and Grondin [257])	$\begin{array}{l}{\epsilon }_v=S_w{P}_c{}^{dyn}\left(\frac{1}{{\kappa }_b}-\frac{1}{{\kappa }_s}\right){\kappa }_s\to +\infty \\ {\epsilon }_v=S_w\left(P_c{}^{stat}-\tau \frac{\partial S_w}{\partial t}\right)\left(\frac{1}{{\kappa }_b}\right)\end{array}$	(63)
Mechanical equilibrium equation	${\sigma }_{ij,j}+{\rho }_{b}b=0$	(64)

where contains the following two sets of parameters: Seepage parameters: ρ_w = the density of soil water, n_e = the effective porosity, S_w = the soil water saturation, t = the time, K = the intrinsic permeability, K_r = the relative permeability, µ_w = the dynamic viscosity of pore water; g = the gravational accelerator; P_w = the pore water pressure, P_a = the pore air pressure assumed zero for unsaturated soil, P_c^stat = the static capillary pressure (equilibrium soil matric suction), P_c^dyn = the dynamic capillary pressure (nonequilibrium soil matric suction), τ = the dynamic coefficient, u = the soil matrix displacement, n₀ = the initial effective porosity, ε_v = the volumetric strain, K₀ = the initial intrinsic permeability, α_stauffer = the Stauffer coefficient with a value of 0.1 for the air-water system, α_BC and n_BC = the Brooks and Corey fitting parameters in Equation (3) in Table 1, S_e = the effective saturation, S_r = the irreductable saturation, α_VG, n_VG and m_VG = the van Genuchten fitting parameters in Equation (4) in Table 1, and p = a coefficient for parameters transformation between Brooks and Corey and van Genuchten functions; Mechanical parameters: σ_ij′ = the effective stress, σ_ij = the total normal stress, b = the Biot’s coefficient, δ_ij = the Kronecker delta (unit diagonal matrix, dimension of matrix depending on 2D or 3D), κ_b = the drained bulk modulus of soil, κ_s = the bulk modulus of solid particles (assumed = +∞ for rigid particles), λ_b = the Lame’s moduli of dry porous media, ε_v = the volumetric strain, ε_xx, ε_yy, and ε_zz = the linear strain on x, y and z directions, G_b = the shear modulus, ε_ij = the strain tensor, u_i or u_j = the displacements on i or j directions (i or j = x, y, z), v_b = the soil Poisson ratio, E_b = Young’s elastic modulus of soil, σ_mean = the mean total normal stress for isotropic loading, σ_xx, σ_yy, and σ_zz = the total normal stress on x, y and z directions, P_c^dyn = the dynamic capillary pressure (nonequilibrium soil matric suction), P_c^stat = the static capillary pressure (equilibrium soil matric suction), τ = the dynamic coefficient, σ_ij,j = the entire stress tensor (dimension of matrix depending on 2D or 3D), ρ_b = the bulk density of soil, and ρ_bb = the total body force including the gravational force (ρ_bg), centrifuge force, Coriolis force for large scale groundwater and geological models, and other field-generated body forces, usually the gravational force exclusively adopted for geotechnical engineering.

describes the theoretical framework developed by Zou, Saad and Grondin [257]. For details of the derivation and assumptions of this theoretical framework, the work from Zou, Saad and Grondin [257] is worth following up on.

The numerical method applied by Zou, Saad and Grondin [257] was the Finite Element Method (FEM). Both transient two-phase seepage and mechanical governing equations were discretised in the spatial domain using FEM and in the temporal domain using an implicit time-stepping method, the backward Euler method. The governing equations were transformed into weak form by the Galerkin method. Zou, Saad, et al. [257] checked the numerical stability and convergence. According to Zou, Saad, et al. [257], the inversion analysis by the Levenberg-Marquardt (LM) algorithm was used for the parameter optimization for Kozeny-Carman and transient effects models. The hydro-mechanical coupling illustration is shown in Figure 11.

Zou, Saad and Grondin [257] successfully validated their transient seepage simulator against experimental results from Zhuang, Hassanizadeh, Qin and de Waal [47] to agree on the log-linear relationship between dynamic coefficient and saturation. In addition, they further carried out the hydro-mechanical coupling simulation with dynamic/transient effects for one-dimensional settlements under lateral constraint (K₀) conditions. As a result, Zou, Saad and Grondin [257] finally concluded that more significant settlements (shrinkage) could be generated by transient effects in comparison to calculating deformation without such effects.

Even though Zou, Saad and Grondin [257] theoretically and numerically contribute to the unsaturated soil hydro-mechanical coupling under transient flow conditions, their work still owns many limitations, which they already reflected. For instance, they expected to couple this model with vapour diffusion, fluid-phase transformation, and, even more significantly, pore matrix deformation rather than elastic settlement, perhaps pore-elastic-plastic theory. Still, it is a distinctive attempt to apply the nonequilibrium transient two-phase seepage theory in unsaturated soil mechanics. Therefore, any future work will be highly expected and will deserve to be followed up on.

8.5. Discussion on Practical Engineering Application

Last but not least, a few complexities and limitations in theoretical expansion could be discussed here regarding theoretical development. As the Soil Water Retention Curve (SWRC) highly depends on a pore matrix and capillarity, the significant deformation will turn the static SWRC into a Soil Water Retention Surface (SWRS), including porosity or void ratio [65,80]. Any attempts to add transient terms into such complicated SWRS will lead to theoretical extravagance and numerical challenges.

Furthermore, no experimental method exists to test and validate the hydro-mechanical coupling process under nonequilibrium transient conditions. With more sensors inserted into soil specimens for detecting transient effects, the mechanical behaviour will also be disturbed, consequentially leading to the inaccuracy of mechanical properties determination. So far, many non-destructive moisture-detecting methods (e.g., CT scan, Gamma-ray, etc.) could replace moisture sensors. However, the pore pressure sensor (e.g., tensiometer) is irreplaceable for detecting nonequilibrium soil matric suction. Thus, it is extraordinarily challenging to couple transient effects with mechanical deformation, referring to those aspects above.

Nonetheless, the assumption on an insignificant variation of pore matrix could be conserved for evaluating unsaturated soil shear strength when considering dynamic effects. It somehow opens a gate for applying stability analysis of earth profile in slope, dam, etc., based on the extreme equilibrium theory. Therefore, theoretical, experimental, and numerical efforts on this application might deserve to be more investigated in the future.

9. Conclusions

This comprehensive review is constructed for the transient two-phase flow in porous media and its engineering applications in Geotechnics. It commences from a literature review of conventional two-phase flow seepage in three sections: fundamental multiphase physics and soil water retention curve, steady-state and transient seepage theories with hydraulic properties, and experimental methods for relative hydraulic conductivity. Then, the constraints in those theories and experiments are challenged based on a review of experimental observations, which violate conventional model predictions and assumptions, followed by a multiphase physio-analysis of possible reasons at pore-scale. Therefore, the dynamic nonequilibrium effects, abbreviated as dynamic or transient effects, are defined for modern investigations. By far, four advanced theories have already been composed for modelling nonequilibrium transient two-phase flow in porous media. This work collects them all for numerical modellers’ interests. Moreover, the advantages and disadvantages of each theory are discussed regarding engineering applications. Modern experimental investigations based on these newly proposed theories have been executed in multiscale. This study further reviews the state-of-art of current studies in five sections: the experimental setup, influential factors in dynamic/transient nonequilibrium effects, and modelling validations, as well as micromodel and numerical approaches at the microscale. Last but not least, the engineering application of nonequilibrium transient seepage into the hydro-mechanical coupling of unsaturated soil is reviewed with discussions on their applicability. Including transient effects into effective stress, shear strength of unsaturated soil, and extreme equilibrium analysis should be practical to stability analysis of natural earth profiles in geotechnical engineering practices, including slope, dam, etc.

Based on this comprehensive review, conclusions can be summarised in the following points:

(1)

The conventional theory of transient two-phase flow in porous media is still valid for transient flow conditions only if the instantaneous equilibrium condition can be achieved by carefully controlling the boundary conditions of geotechnical tests.

(2)

There are still many research questions left for instantaneous equilibrium two-phase flow in porous media in terms of deformation coupled with hysteresis, the high suction range for soft soil, hysteretic hydraulic properties for problematic soil, etc.

(3)

The conventional experiments determine the soil water retention and hydraulic properties under the assumption of instantaneous equilibrium. However, when such an assumption is violated, the inversion analysis can still produce velocity-dependent hydraulic diffusivity. Lacking a physical basis, this scenario becomes a fitting process rather than a physical and experimental characterization.

(4)

Earlier experimental observations of dynamic nonequilibrium effects under transient two-phase flow conditions fostered the theoretical development for modelling such effects. Thus, it is expected to advance theories with more experimental findings in multiscale.

(5)

The four advanced theories have their advantages and disadvantages. The soil water redistribution model, dual-fraction model with soil water redistribution, and dual-porosity and dual-permeability model can simulate transient effects in multiscale. However, according to this listing sequence of models, the theoretical complexity increases with more parameters in the governing equations. Except that the soil water equilibration time can be straightforwardly determined using instrumented soil column test, other parameters can only be given by inverse modelling.

(6)

Compared to the previous three theories in soil hydrology, the thermodynamic-based theory can be applied to simulate both soil moisture dynamics and nonequilibrium soil suction under transient two-phase flow conditions. Therefore, it has a unique advantage for estimating the mechanical properties of unsaturated soil, whereas others neglect this critical application.

(7)

Modern continuum-scale experimental methods often incorporate soil suction and moisture sensors into a soil column with various hydraulic boundary conditions. In petroleum engineering, core flooding tests can be implemented with non-destructed measuring methods for fluid saturation (e.g., CT scan, Gamma-ray, outflow, etc.). Nevertheless, pore pressure transducers are irreplaceable for determining dynamic capillary pressure.

(8)

Recent continuum-scale experimental and numerical studies focus on investigating influential factors of dynamic nonequilibrium effects in terms of dynamic capillary coefficient. The influential factors mainly include properties of porous media, fluid properties, multiphase physical properties, etc. However, many conflicting conclusions were drawn amongst those studies for each influential factor, therefore, needing more experimental revisits. Furthermore, the hysteretic and extreme conditions for dynamic effects should be continuously investigated using experimental and numerical methods in the future.

(9)

All advanced theories have been validated successfully by many experimental studies on specific hydraulic loading paths. However, there are failure cases, which might be due to numerical solution, parameter selection, etc. So, it is still worth casting the numerical tools constantly in order to improve the accuracy of numerical solutions. Moreover, the fingering flow simulated by the thermodynamic-based theory considering hysteresis could be due to conditionally numerical instability, which is quite debatable for mathematicians and numerical modellers.

(10)

The pore-scale modelling transient two-phase flow displacement can be achieved by the physical micromodel, numerical Pore-Network Model (PNM), and Computational Fluid Dynamic (CFD) method. Each is a powerful tool to investigate transient two-phase flow from pore-scale up to Representative Elementary Volume (REV) scale or even several REVs. Still, the physical micromodel is not often available in every geotechnical laboratory and is not generic. On the other hand, the Dynamic PNM (DPNM) can simulate the dynamic effects with lower computational expenses than CFD. However, the interfacial dynamics was not counted at pore-scale. Therefore, only local heterogeneities induced dynamic/transient effects can be reproduced by DPNM.

(11)

CFD methods include multiphase Lattice Boltzmann Method (LBM), Direct Numerical Simulation (DNS) coupled with Volume of Fluid (VoF) or Level-Set (LS), and Smoothed Particle Hydrodynamics (SPH). Only a few numerical experiments using those methods already replicated dynamic effects in dynamic capillarity, capillary-viscous and inertial dominating flow conditions. Although those numerical methods have been implemented to study steady-state flow conditions in flow patterns and regimes, the contribution to dynamic effects under transient flow conditions are still rare and, therefore, strongly urge pursuit in the future.

(12)

The observation of dynamic nonequilibrium effects in natural slopes was firstly reported in the recent literature. It partially supports the application of transient effects for engineering practices. However, the in-situ data are still insufficient to prove the importance of considering dynamic effects in modelling transient two-phase flow seepage. It is expected to receive more contributions from the field or large-scale observations than smaller-scale findings using the pressure or flow cell tests.

(13)

The dynamic/transient effects have been coupled into unsaturated soil effective stress, soil suction characteristic stress, subsequently in shear strength, and pore-elastic consolidation. Those studies proposed the theoretical frameworks for hydro-mechanical coupling with transient effects and succeeded in simulating nonequilibrium transient water flow in unsaturated soil. However, the validation can only be given to the transient seepage, while the mechanical prediction cannot be verified by any experimental results yet.

(14)

There are many challenges in setting up an experimental apparatus to test the mechanical properties of unsaturated soil under transient flow conditions. An inevitable drawback is the irreplaceability of tensiometers in detecting nonequilibrium soil suction. The insertions of those sensors will cause mechanical properties variation in terms of deformation. Moreover, deformation-induced nonuniqueness of soil water retention curve will dramatically increase complexity for additional coupling with transient effects. Therefore, those reasons constrain hydro-mechanical investigation when transient effects are considered. However, it seems still applicable for shear strength estimation and extreme equilibrium analysis when the soil pore matrix varies insignificantly.

To conclude, there is no doubt that extreme environmental conditions are unnecessary in the practical geotechnical design when considering cost and efficiency. Still, it can be used to predict any potential failure for a warning system in the future in order to save lives and costs.