Soil hydraulic conductivity in the state of nonequilibrium

Hydraulic nonequilibrium in soil during water infiltration and drainage is a well‐known phenomenon. During infiltration, water initially invades easily accessible pores before it slowly redistributes towards some state of energetic minimum. In analogy, during drainage, easily drainable pores are emptied more rapidly than those blocked by bottlenecks. The consequence is that the water content is lagging behind the water potential and both state variables do not follow a unique water retention curve as typically assumed when applying Richards equation. Current models that account for nonequilibrium allow for the required decoupling of water content and water potential; however, they do not consider the consequences for the hydraulic conductivity. In this contribution, we present a physically based approach to estimate hydraulic conductivity during nonequilibrium, which depends on both water content and water potential during nonequilibrium conditions. This approach of a dynamic hydraulic conductivity function is demonstrated for an infiltration process into relatively dry soil and for a stepwise drainage and rewetting with decreasing and increasing water fluxes (i.e., multistep flux experiment). The new approach reproduces well‐known phenomena such as pressure overshoot and preferential flow across infiltration fronts using a unified concept for hydraulic conductivity. This was not possible with existing models assuming some fixed unsaturated conductivity function depending on either water content or water potential.

It can be described by a unique water retention curve θ(ψ) with water content θ and water potential ψ corresponding to capillary pressure. This situation is accomplished when the wetting and nonwetting phases are distributed according to the Young-Laplace equation with the same curvature of water-air interfaces everywhere while the air phase is at atmospheric pressure everywhere. In practice, however, this state is probably never reached because soils are open systems where the time required to establish this equilibrium is much longer than Vadose Zone J. 2023;22:e20238.
wileyonlinelibrary.com/journal/vzj2 1 of 12 https://doi.org/10.1002/vzj2.20238 the typical time scale of changes in water content as triggered by the atmospheric boundary condition. However, in modeling soil water dynamics using state-of-the-art physically based models such as Richards equation, we implicitly assume such an equilibrium by assuming some fixed θ(ψ) function that is considered to be a material property. In some cases, hysteresis is taken into account (Hannes et al., 2016), which is a well-known feature. It is typically explained by the heterogeneous porous structure that allows for the same water potential at different water contents due to phase entrapments (Poulovassilis, 1962).
Generally, it is widely accepted that the prerequisites for the application of Richards equation are that the porous system is rigid, nonswelling, and isotropic, and that liquid water flow is isothermal and air can move freely at atmospheric pressure. Another prerequisite is typically taken as grantednamely, that the constitutive relations θ(ψ) and the hydraulic conductivity function [K(θ) or K(ψ)] actually exist with or without hysteresis. It is well recognized that there is considerable uncertainty in how to measure and parameterize these relations and there is an enormous body of literature on the different methodological aspects. In contrast, it is rarely questioned whether moving water can, in principle, adhere to these constitutive relations and how to cope with the situation if not. This problem has been discussed in the framework of hydraulic nonequilibrium by many authors as cited below; however, all the proposed concepts refer to the nonequilibrium between θ and ψ (i.e., the water retention curve), but the implications for hydraulic conductivity were completely ignored so far. The objective of this paper is to discuss these implications and to develop an alternative to formulate unsaturated hydraulic conductivity as a function of θ and ψ during nonequilibrium conditions.

EXPERIMENTAL EVIDENCE
Since the pioneering work of Topp et al. (1967), followed by Smiles et al. (1971) and Vachaud et al. (1972), who were probably the first who systematically studied the nonuniqueness of θ(ψ), it became obvious that θ(ψ) obtained from drainage and infiltration experiments depend on the dynamics (i.e., the rate of change in water content). Hence, the relation between θ and ψ should rather be described by a possibility space than a unique function. The observed phenomena were later summarized under the term "dynamic effects" (Hassainzadeh et al., 2002). Diamantopoulos et al. (2012) reviewed the various phenomena that are attributed to the "dynamic nonequilibrium of water flow in porous media," as they termed it. This included effects of dynamic contact angles, air entrapment, and the relaxation of air-water interfaces. Typical experiments to demonstrate dynamic nonequilibrium are drainage and imbi-

Core Ideas
• Soils are rarely in hydraulic equilibrium.
• We show consequences for their effective hydraulic conductivity. • We present a physically based concept how to better describe the unsaturated conductivity function. • The new approach describes pressure overshoot across fronts and the emergence of preferential during infiltration.
bition experiments performed at different gradients in water potential or in water flow, as done by Wildenschild et al. (2001), who observed more pronounced dynamic effects for coarse-textured materials. Most experiments are pressure controlled, and the generally observed feature is that the change in water content is lagging behind the change in water potential. The general understanding is that during drainage or imbibition, the water within a soil volume does not instantaneously distribute itself within the pore structure according to some state of energetic minimum. This is due to the structural complexity of the pore system, which allows some parts of the pore system to equilibrate faster than others and due to the fact that some of the available energy is required to move the interfaces within the porous medium. Thus, there is some characteristic time scale for equilibration that depends on soil structural properties.
Another interesting feature is observed in experiments that are flux controlled while θ and ψ can freely adapt to the imposed water flux. Stonestrom and Akstin (1994) measured the course of the water potential during infiltration into different materials from fine-textured soil to glass beads. They observed a clear overshoot in water potential during passage of the infiltration front and could exclude effects from entrapped air or particle rearrangement. Interestingly, this effect was considerably more pronounced if the infiltration rate was increased stepwise as compared with infiltration into initially dry materials.
This work was mainly recognized by researchers working on fingering flow where saturation and pressure overshoot are well known phenomena (Glass et al., 1989). However, Stonestrom and Akstin (1994) studied also fine-textured soil where fingering is not observed. The pressure overshoot was also found in multistep-flux experiments (Weller et al., 2011), which were originally designed to measure unsaturated hydraulic conductivity. In a follow-up paper, Weller and Vogel (2012) investigated the θ-ψ trajectory during such transient, flux-controlled experiments and demonstrated the overshoot of θ and ψ across moving water fronts. This clearly indicated that the same hydraulic conductivity is obtained with different combinations of θ and ψ. Consequently, it seems inappropriate to assume some fixed hydraulic conductivity function, whether it is K(θ) or K(ψ).
Despite all this experimental evidence and theoretical considerations, it remains unclear which processes under which conditions are most relevant. Diamantopoulos et al. (2012) concluded that "a generally valid definition of DNE (dynamic nonequilibrium) thus seems difficult and we propose to think of DNE as hydraulic state observations that cannot be described by a continuum model (Richards equation or a twophase flow model) with a set of unique SHPs (soil hydraulic properties)." This is certainly a valid conclusion. To proceed, we should either decouple water content and water potential in modeling, as we propose in this paper, or we stick to the hope that all these effects are of limited relevance in practice so that we can go on with the assumption of hydraulic equilibrium.
It is still an open question how important these nonequilibrium effects actually are if we want to predict the water dynamics at the scale of soil profiles. Can they be neglected for practical applications of Richards equation? Or are they the main reason why Richards equation often fails to predict water dynamics in natural soils? In fact, as recently pointed out by Vogel (2019), prediction of infiltration and drainage of field and lysimeter experiments are often not convincing even when based on calibrated parameter sets.

MODELING APPROACHES
The currently available modeling approaches for hydraulic nonequilibrium introduce an additional dynamic term to Richards equation. Thereby, the dynamics of the water content are retarded and strive towards equilibrium according to the equilibrium water retention curve by first-order kinetics. This is, in principle, the same for the dynamic capillary pressure proposed by Hassanizadeh and Gray (1993) and the model introduced by Ross and Smettem (2000). The latter was extended by Diamantopoulos et al. (2012) towards a dualdomain approach to account for the observation that some part of infiltrating or draining water rapidly adapts to the change in water potential while another part is lagging behind.
In this way, a number of different experiments (multistep outflow and multistep flux) could be successfully described by Diamantopoulos et al. (2015) after parameter calibration. This also included the overshoot in water potential during flux-controlled experiments. In all these model concepts, nonequilibrium is introduced in terms of the water retention characteristic θ(ψ), whereas hydraulic conductivity is still assumed to be some well-defined function of either θ or ψ. In fact, Ross and Smettem (2000) assumed a fixed K(ψ) function to reproduce preferential flow phenomena for one-dimensional infiltration simulations. In contrast, Diamantopoulos et al. (2015) used a fixed K(θ) function where θ was the total water content including the nonequilibrium domain. In this way, the overshoot in water potential across infiltration and drainage fronts could be reproduced because hydraulic conductivity is changing along with the changing nonequilibrium water content. By choosing K(ψ), it is possible to reproduce preferential flow because with increasing water potential the conductivity gets high while the water content is still relatively low. In contrast, by choosing K(θ), it is possible to reproduce the observed pressure overshoot across infiltration and drainage fronts under flux controlled conditions because the increasing water content behind the front leads to an increase of hydraulic conductivity that is compensated by a decrease in water potential to adapt to the imposed flux.
Considering the water distribution at the pore scale, it can be shown that both concepts, K(ψ) and K(θ), are physically not consistent: the first means that we assume the same conductivity for a given water potential even though the water content is different. The latter means that we assume the same conductivity at a given water content even though this water resides in different pore size classes at different water potentials. Hence, it is crucial to investigate the change in hydraulic conductivity during nonequilibrium. In the following, we propose a physically based concept, how conductivity can be estimated for states of nonequilibrium.

HYDRAULIC CONDUCTIVITY IN THE STAT OF NONEQUILIBRIUM: THEORY
We start with the basic assumption that the pore structure is rigid and that for any water content there is a unique state of minimum energy that defines the water potential in equilibrium. Hence, an equilibrium retention curve exists in principle which, however, might be hardly reached in nature due to nonequilibrium and hysteresis. Then, we consider water infiltration into soil that is initially at hydrostatic equilibrium marked by some point on the equilibrium retention curve at matric potential ψ i and water content θ i (see Figure 1). The arriving infiltration front increases the matric potential to ψ f and the water content to θ f . During nonequilibrium, water content is lagging behind and, thus, θ f is lower than the equilibration water content at ψ f so that θ i < θ f < θ feq (ψ f ) as illustrated in Figure 1.
This transient hydraulic state characterized by ψ f and θ f is expected to depend on the rate of change in water saturation. If this is very low the infiltrating water may distribute according to the equilibration curve. If it is high, it is expected to deviate from the equilibrium curve to some extent. This is reflected in the modeling approaches of Hassainzadeh et al. (2002) and Ross and Smettem (2000). The question which has not been F I G U R E 1 Scaling of the equilibrium retention curve (solid red line) to describe a state of nonequilibrium defined by ψ f , θ f . The red dashed line represents the scaled curve according to Equations 1 and 2 with the blue section describing the nonequilibrium part of the actual state. ψ i , initial water potential; θ i , corresponding water content at equilibrium; θ feq , equilibrium water content at ψ f addressed so far is: what is the hydraulic conductivity at such a state of nonequilibrium? As explained above, it cannot be taken from some fixed K(θ) or K(ψ) function.
For infiltration, K(θ) would underestimate conductivity because the actual water potential at θ f is higher than the equilibrium potential for this water content so that larger pores contribute to water flow. In contrast, if we choose some fixed K(ψ) function, the conductivity is overestimated because the equilibrium water content at ψ f is higher as compared with θ f . Below, the conductivity for the transient hydraulic state ψ f , θ f is estimated based on a plausible distribution of the water phase within pores of different size. We assume a rigid pore structure so that the overall pore size distribution does not change. To characterize the hydraulic state at ψ f , θ f , we need to define a transient water retention curve connecting the state ψ f , θ f with the initial state ψ i , θ i . In this way the water distribution in the interval [ψ f , θ f ] and herewith the distribution of pore sizes that additionally contribute to flow can be described.
To do so, we follow the concept of Kool and Parker (1987) to define some virtual values for the residual water content, θ * r , and the saturated water content, θ * s to scale the water retention curve in the interval [ψ i , ψ f ] based on the shape parameters (n, α) of the equilibrium curve: and where S(ψ f ) and S(ψ i ) are the water saturations according to the equilibrium curve at water potentials ψ f and ψ i , respectively. For the equilibrium retention curve we chose a parameterization according to (van Genuchten, 1980): with m = 1 -1/n . An example for n = 1.6, α -0.006 cm −1 , θ s = 0.44, and θ r = 0.05 is shown in Figure 1.
It should be noted that this approach represents just one possible way to describe the nonequilibrium part of the water retention curve. It is chosen because of its operational simplicity. It implicitly assumes that all available pores that can be filled with water in the interval [ψ i , ψ f ] are filled at the same speed irrespective their size. The radii of these pores can be calculated based on the Young-Laplace equation: = 2σ cos α∕ψ (4) with surface tension σ and assuming perfect wettability with the contact angle α = 0˚. It should be noted that the velocity of capillary filling is expected to be proportional to the radius of the capillary according to the Washburn equation (Washburn, 1921). However, to reduce complexity, this is not considered here.
To estimate the relative hydraulic conductivity, we integrate over the water filled pore space following the concept of (Mualem, 1976): with the tortuosity factor τ = 0.5. The second term in the nominator reflects the contribution of the nonequilibrium part to the relative conductivity which is scaled by the factor a i according to the reduced saturation in the interval [ψ i , ψ f ] given by The relative conductivity K r calculated in this way for the state ψ f , θ f is substantially higher as compared to the value obtained for K r (θ f ) in hydraulic equilibrium because of the fact that larger pores are involved in water flow. In contrast, this conductivity is substantially lower than K r (ψ f ) in equilibrium because the water saturation in the state of nonequilibrium is substantially lower. The values are given in Table 1.
Next, we consider drainage starting from some point on the equilibrium curve. It should be noted that under natural T A B L E 1 Relative hydraulic conductivities in the nonequilibrium state (ψ f , θ f ) as illustrated in Figures 1 and 2 calculated based on the classical equilibrium functions K r (θ) and K r (ψ) and compared with the nonequilibrium approach K r (θ, ψ) according to Equations 5 and 9 and the alternative approach for drainage * r (θ, ψ) according to Equation 11 Approach Note. ψi, initial water potential; θi, initial water content.

F I G U R E 2
State of nonequilibrium during drainage characterized by ψ f , θ f . The blue line indicates a possible pathway starting from the initial state ψ i , θ i obtained from scaling of the retention curve (red dotted line) in analogy to Figure 1. Note: ψ i , θ i , is the initial water potential and water content at equilibrium; θ feq , equilibrium water content at ψ f ; ψ feq , equilibrium water potential at θ f conditions, drainage processes are expected to be less abrupt as infiltration could be. However, the multiflux experiments of Weller et al. (2011) clearly showed that a decrease in irrigation leads to a fast decrease in water potential accompanied by a pronounced dynamic effect with pressure overshoot and with the water content lagging behind the water potential. In analogy to infiltration, this nonequilibrium is again described by some transient hydraulic state ψ f , θ f deviating from the equilibrium curve as illustrated in Figure 2. A possible transient water retention curve linking the initial state ψ i , θ i with the state at the drainage front ψ f , θ f is obtained by scaling of the retention curve as was done for infiltration. This leads to virtual values for θ * r and θ * s (details are provided in the Appendix): and In this way, we obtain the blue segment for the nonequilibrium part in Figure 2.
To estimate the hydraulic conductivity, we again integrate over the water filled pore space following the concept of Mualem (1976). For water potentials smaller than ψ f , the entire pore space is water filled. The water volume in the interval between the actual water content θ f and the equilibrium water content θ freq is expected to be contained in pores that drain in the interval between ψ i and ψ f . Hence, the relative conductivity can be calculated in analogy to what was proposed for infiltration as where the second term in the nominator is scaled by the factor to account for the reduced volume of water within the integration interval and to preserve the total water capacity. This approach implicitly assumes that pore of all sizes that can be drained in the interval [ψ i , ψ f ] contribute in the same way to the water loss between θ i and θ f . Strikingly, the relative conductivity K r,drainage (θ f , ψ f ) calculated by Equation 9 for the example shown in Figure 2 is not only higher than K r (ψ f ) as expected but also higher than K r (θ f ). This is because in the transient state (ψ f , θ f ), larger pores that are in the interval [ψ i , ψ eq ] contribute to water flow in contrast to K r (θ f ) at equilibrium where these pores are completely drained. However, the assumption that pores in the interval between ψ i and ψ f contribute to conductivity according to their equivalent pore diameter is not in the same way justified for drainage as compared with infiltration. During infiltration, the water volume that is in nonequilibrium (i.e., in the interval between ψ i and ψ f ) has infiltrated pores that are easily accessible. In contrast, during drainage, the water volume in the interval between ψ i and ψ f is contained in pores that can be drained only slowly. Hence, these pores can be assumed to be less accessible or only connected via bottlenecks in the size range of the equivalent pore diameter of ψ f .
To take these considerations into account, we assume that the water-filled pores in the interval between ψ i and ψ f are arranged in sequence so that the harmonic mean of the pore diameter should be a suitable estimate for the hydraulically effective pore size. This leads to the solution that the integration over the nonequilibrium part of the water filled pore space is not done from ψ i and ψ f, but only from ψ* to ψ f with ψ* = (ψ i + ψ f )/2. Thus, Equation 9 is slightly modified to * This approach leads to reduced hydraulic conductivity during nonequilibrium, which lies in between the values obtained for K r (θ) and K r (ψ) ( Table 1).
Based on Equations 5 and 11, the relative conductivities for the entire range of possible states during nonequilibrium can be calculated. The results depend on the assumption of the initial state ψ i , θ i prior to infiltration or drainage. For infiltration, θ i should be in the interval [0, θ f ], whereas for drainage it is in the interval [θ s , θ f ]. As a physical interpretation, the jump from θ i to θ f represents the jump in water content across the drainage or the imbibition front. To illustrate the effect of the size of this jump, Figure 3 shows the results for three different scenarios. For imbibition the initial water content is calculated by bθ f , whereas for drainage it is obtained by θ f + (1 -b)(θ s -θ f ) with 0 ≤ b ≤ 1. With increasing b, the jump across the infiltration or drainage front is decreasing.
There are some important general features to note. At any state of nonequilibrium during imbibition, the hydraulic conductivity is increasing when approaching the equilibrium water content during an equilibration period (i.e., parallel to the water content axis at potential ψ f as implemented in the model of Ross & Smettem, 2000). For drainage processes, this is vice versa. Hence, for some constant flow rate, the water potential is expected to decrease during imbibition and it is expected to increase during drainage. Generally, the ψθ-trajectory is expected to follow the iso-conductivity lines as indicated in Figure 3. For drainage and imbibition, this explains the pressure overshoot as observed in flux-controlled experiments where water content and water potential can adapt freely, as shown experimentally (Stonestrom & Akstin, 1994;Weller & Vogel, 2012;Weller et al., 2011). The lower right panel in Figure 3 illustrates the result for the assumption made by Equation 9 (i.e., all pores contribute to flow according to their equivalent diameter). This is considerably different compared with the more realistic assumption of less connected pores according to Equation 11. This is also confirmed by the experimental findings of Weller and Vogel (2012), who observed an increasing water content during equilibration in a drainage process while Equation 9 would predict a decreasing water saturation.

IMPLEMENTATION
The proposed concept provides a physically based approach to estimate hydraulic conductivity in the state of hydraulic nonequilibrium. The implementation of this concept for numerical simulation of water dynamics based on Richards equation would require calculation of the conductivity fields under transient conditions while the structure of these fields change depending on the history and the gradient in water content and water potential across drainage and infiltration fronts.
For field situations, this could be limited to infiltration events where steeper gradients in water potential may occur. This is much less the case for drainage processes, which are typically very slow. For laboratory experiments such as multistep outflow, however, also sharp drainage fronts need to be considered.
One option for implementation would be to recalculate the conductivity fields for any infiltration (or drainage) front. This would come along with a considerable computational burden. As an alternative, precalculated conductivity fields could be used for a limited number of pressure jumps in combination with suitable interpolations. These concepts still need to be developed.
To demonstrate the effect of choosing a dynamic function K(ψ,θ) instead of some static K(ψ) or K(θ), we implemented a simplified approach. As we have seen (Table 1), realistic values for K(ψ,θ) during drainage and imbibition range between the values of K(ψ) or K(θ). Hence, as an approximation we propose a weighted geometric average between these conductivities: and p ∊ [0,1]. The results are illustrated for different values of p in Figure 4c-e in comparison with those for a fixed K(ψ) or K(θ) (Figures 4a,b). Obviously, the distortion of the unsaturated conductivity field looks similar to the physical based approach illustrated in Figure 3. However, in contrast with that, there is a smooth transition from imbibition (below the water retention curve) to drainage (above the water retention curve). A probable explanation for the kink in the conductivity field at the transition between infiltration and drainage ( Figure 3) is the same as we already discussed above for the difference in the effective conductivity during infiltration and drainage. During infiltration, the larger pores that are invaded by water are those that are easily accessible and well connected, so that K(ψ) is closer to the realistic value. In contrast,  Figures 1 and 2 and the corresponding initial values on the equilibrium curve with rhombs for infiltration and triangles for drainage. The white line corresponds to the equilibrium retention curve, and the color code represents log(K r ), where K r is relative conductivity. Note that the area above the equilibrium retention curve represents drainage and the area below this curve represents infiltration during drainage the larger pores that are still filled with water are those which are less connected. Hence, K(θ) might be closer to the realistic value. As an approximation we suggest an inverse weighting with p = q for imbibition and p = 1 − q for drainage with q ∊[0,1]: An example for q = 0.8 is shown in Figure 4f. This resembles closely the behavior shown in Figure 3.

NUMERICAL SIMULATIONS
We used the nonequilibrium model proposed by Diamantopoulos et al. (2015) to demonstrate the effect of the chosen model for hydraulic conductivity on the dynamics of infiltration and drainage. It is based on the approach of Ross and Smettem (2000) but distinguishes an equilibrium domain and a nonequilibrium domain. This accounts for the observation that during transient conditions along a drainage or a wetting front a part of the water content is in instantaneous equilibrium according to the equilibrium water retention curve while another part is not. This behavior can be observed for example in multistep-outflow experiments (Diamantopoulos et al., 2012). The Richards equation is modified to account for the non-equilibrium domain as follows: (1 − ne ) ∂θ eq (14) where f ne ∊[0,1] denotes the nonequilibrium fraction. The expected equilibrium water content at the given water potential, θ eq , is reached only partly (i.e., by the fraction 1 − f ne ), whereas the actual water content, θ ne , is in nonequilibrium. The hydraulic conductivity is now depending on both ψ and θ according to Equation 12. Following Ross and Smettem  Figures 1 and 2, and the corresponding initial values on the equilibrium curve with rhombs for infiltration and triangles for drainage. The white line corresponds to the equilibrium retention curve, the color code represents hydraulic conductivity, log(K r ) (2000), the dynamics towards equilibration is described by a simple linear function of the difference between the actual water content and that at equilibrium. The dynamics towards equilibrium is described by the rate parameter γ[T]: We conducted numerical simulations of infiltration and drainage based on the same van Genuchten parametrization of the equilibrium water retention curve as used above in the section describing the theory. The saturated hydraulic conductivity was K s = 140 cm d −1 , and the tortuosity factor in Equation 11 was chosen as τ = 0.5. The nonequilibrium fraction was set to f ne = 0.5 and the equilibration time was γ = 15 d. All parameters used are summarized in Table 2. Note. n and α, van Genuchten parameters; θ s , saturated water content; θ r residual water content; K s , saturated hydraulic conductivity; τ, tortuosity parameter; f ne , nonequilibrium fraction; γ, rate parameter for equilibration.
First we simulated a typical multistep-flux experiment where the flux at the upper boundary is changed stepwise. The results for different weighing of K(ψ) and K(θ) using the parameter q (Equation 13) is shown in Figure 5. With increasing q, the weight is more on K(ψ) for drainage and more on K(θ) for infiltration. The overshoot in water potential can be observed for both drainage and infiltration, but in opposite directions, which corresponds to what is typically observed in flux-controlled experiments (Diamantopoulos et al., 2015). The overshoot is more pronounced with increasing weight on K(θ) and cannot be reproduced by choosing a fixed K(ψ) function. This case is represented by q = 1 during infiltration. The difference between q = 1 and the standard Richards equation [fixed K(ψ)] during increasing flux steps in Figure 5 is due to the dynamic nonequilibrium model (DNE) where the water content reaction is delayed, leading to a deviation of the hydraulic state from the equilibrium water retention curve. The differences in the dynamics of water content are minor for different values of q. In the DNE model, the time to reach equilibrium is prolonged as compared with standard Richards equation.
In another example, we simulated a heavy rainfall event (55 cm d −1 ) on a relatively dry soil at an initial water potential of h = −1.000 cm. The depth profiles of hydraulic state variables after 2.4 h of infiltration are shown in Figure 6. As already shown by Ross and Smettem (2000), who used a static K(ψ) function, the water infiltrates more rapidly at lower water contents reflecting preferential flow initiated by nonequilibrium. More specifically, the water at the infiltration front enters larger pores at higher water potentials but lower water contents than when assuming standard Richards equation. Consequently, the infiltration front is moving faster. It should be noted that for values of q smaller than q ≈ 0.6, the water potentials become positive (results not shown). This is because with decreasing values of q, the relative weight of K(θ) is increasing and because θ is lagging behind, the highest achievable hydraulic conductivity is considerably lower than the static value of K s . Hence, for low values of q, the reduced K(ψ,θ) cannot reach the flow rate imposed at the upper boundary and the water potential becomes positive. The implications for natural boundary conditions is that for a heavy rain on a relative dry soil, the maximum infiltration rate might be well below the saturated conductivity. Figure 7 depicts the trajectories of the hydraulic state variables in the θ − ψ space for q = 0.8 during the multistep-flux F I G U R E 6 Depth profiles of the hydraulic state variables h and θ after t = 2.4 h of infiltration at a flow rate of 55 cm d −1 for the same parameter settings as shown in Figure 5. The dynamic nonequilibrium model (DNE) using the hydraulic conductivity function K(ψ, θ) and different weighing factors q is compared with standard Richards equation assuming fixed hydraulic properties according to a vanGenuchtem-Mualem (VGM) F I G U R E 7 Trajectories of water content (θ) and water potential (ψ) (symbols) for weighing factor q = 0.8. Left: the white line represents the equilibrium water retention curve, the black line the trajectory during infiltration during both simulated experiments, and the gray line is the trajectory during decreasing flux steps in the multistep-flux experiment. The color code describes the unsaturated hydraulic conductivities which corresponds to Figure 4d. Right: the same trajectories (dashed lines) relative to the static water retention curve ψ(θ) (blue) and hydraulic conductivity function K(θ) (red) experiment and during infiltration (Figure 7) at the depth of 5 cm. It shows clearly how the state of nonequilibrium is reached after each sudden change in flow rate. This is followed by an equilibration phase where the trajectory follows an isoline in the hydraulic conductivity field. If some static function of K(ψ) or K(θ) is assumed, as is typically done, this equilibration would be parallel to the θ axis or to the ψ axis, respectively, as illustrated in Figures 4a and 4b. There is no pressure overshoot possible when assuming K(ψ), and there is no preferential flow during nonequilibrium when assuming K(θ). The proposed concept that the hydraulic conductivity depends on both state variables during nonequilibrium is physically justified and reproduces all the observed phenomena. It should be noted that the trajectory shown in Figure 7 also illustrates that nonequilibrium leads to a hysteresis in the dynamic water retention curve. This suggests that the observed phenomenon of hysteresis could also be explained by hydraulic nonequilibrium. However, this will be the subject of future research.
The approach proposed here adheres to the Richards equation although the state variables are decoupled. This leads to the critical question whether, in this case, the state variables can be averaged within the discretization volume for the numerical solution of the Richards equation or if this violates the underlying physics. We think that for small-scale heterogeneities of the pore structure in the size range of millimeters to centimeters, pressure balance within the water phase is sufficiently fast so that a mean potential and mean water content within this volume can be reasonably specified. This is why we think the proposed approach is justified even though the state variables are not in equilibrium with respect to some static water retention curve corresponding to the energetic minimum. This situation changes when it comes to macropore flow produced along structures at a larger scale. Here, water moves at high potentials in the macropores while the soil matrix between the macropores is still relatively dry. In this case, the meaning of an average water potential is not obvious, and multidomain models (Gerke & van Genuchten, 1993;Jarvis et al., 1991) might be more suitable because they represent the existing potential differences. However, the modeling approach proposed here shows that even for small-scale heterogeneities, the velocity of infiltration and drainage fronts can be considerably higher than predicted by the classical Richards equation.
Finally, a practical question is how to deal with highly transient situations when equilibrium is not reached prior to a subsequent infiltration or drainage event. For the numerical simulation presented here, such an equilibration was actually reached and drainage and infiltration processes could be clearly distinguished. Specifically, if the actual water content was below the equilibrium water content at a given potential, the soil was in infiltration mode, and if it was above the equilibrium water content, it was in drainage mode. This clear separation is, however, not obvious for highly transient situations. For this case, the concept presented here might still be applicable, whereas further assumptions on the distribution of water within the pore space might be appropriate.

CONCLUSIONS
Hydraulic nonequilibrium during infiltration and drainage is a well-known phenomenon in soil hydrology. However, nonequilibrium effects pertain not only to the soil water retention function but also to the unsaturated hydraulic conductivity function. Depending on the steepness of the drainage or infiltration front, this can lead to substantial changes in the water distribution in the soil profile as compared with predictions by standard Richards equation. Current models that account for nonequilibrium allow for the required decoupling of water content and water potential; however, they do not consider the consequences for the hydraulic conductivity.
In this contribution, we present an extended nonequilibrium approach that accounts for the effect of hydraulic nonequilibrium on the unsaturated hydraulic conductivity. Based on physical considerations of the distribution of water within pores of different size, the unsaturated conductivity was derived for any combination of water content and water potential. A simplified version of this approach was built into an existing nonequilibrium model. As a result, the well-known phenomena of pressure overshoot and preferential flow across infiltration and drainage front could be reproduced. This was not possible by existing models of hydraulic nonequilibrium that do not account for the effect of nonequilibrium on hydraulic conductivity. A challenge for future work will be how to implement the developed concept without the simplification chosen here and to describe the transition between drainage and irrigation when the soil is not in hydraulic equilibrium at this transition point.

A C K N O W L E D G M E N T S
We thank the Deutsche Forschungsgemeinschaft (DFG) for financial support under grant VO 566/18-1, WO 1781/4-1, and GE 990/13-1 "Vadose Zone Modelling of water flow in hillslope soil (VAMOS)."

C O N F L I C T O F I N T E R E S T
The authors declare no conflicts of interest.