Modeling dynamic non-equilibrium water flow observations under various boundary conditions

Modeling dynamic non-equilibrium water flow observations under various boundary conditions

Journal of Hydrology xxx (2015) xxx–xxx Contents lists available at ScienceDirect Journal of Hydrology journal homepage: www.elsevier.com/locate/jhy...

1MB Sizes 0 Downloads 55 Views

Journal of Hydrology xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Modeling dynamic non-equilibrium water flow observations under various boundary conditions E. Diamantopoulos a,⇑, W. Durner a, S.C. Iden a, U. Weller b, H.-J. Vogel c a b c

Technische Universität Braunschweig, Institut für Geoökologie, Germany UGT GmbH, Freising, Germany Helmholtz-Zentrum für Umweltforschung UFZ, Halle, Germany

a r t i c l e

i n f o

Article history: Received 23 March 2015 Received in revised form 19 July 2015 Accepted 23 July 2015 Available online xxxx This manuscript was handled by Corrado Corradini, Editor-in-Chief, with the assistance of Xunhong Chen, Associate Editor Keywords: Dynamic effects Non-equilibrium water flow Capillary pressure saturation relationship Dynamic Modeling

s u m m a r y Experimental evidence of variably-saturated flow in soils often shows a non-equilibrium between water content and water potential. Such phenomena, summarized by the term ‘‘dynamic effects’’ often occur in laboratory experiments specifically designed for the identification of soil hydraulic properties which are a mandatory input for water flow simulations at the field scale using the Richards equation. Depending on boundary conditions, dynamic non-equilibrium leads either to a relaxation of pressure head for constant water content or vice versa. Diamantopoulos et al. (2012) presented a model which describes dynamic effects in the case of multistep outflow experiments by assuming that the total water content is portioned into two fractions: one fraction which is in instantaneous equilibrium with the pressure head, and a second fraction for which the equilibration of the water content with the pressure is described by a first-order kinetics which leads to a temporal shift in water content at constant pressure head. In this work we applied this non-equilibrium model to four different observations of dynamic non-equilibrium effects presented in the literature and demonstrate that it successfully describes (1) non-equilibrium in water absorption experiments which results in a non-uniqueness of the diffusivity vs water content function, (2) the flow-rate dependence of soil hydraulic properties, (3) dynamic effects in multistep outflow and (4) dynamic effects in multistep flux experiments. We conclude the article by a discussion of possible macroscopic explanations, the consideration of dynamic effects for modeling water flow at the field scale, the uniqueness of the estimated non-equilibrium parameters among different experiments, and further theoretical considerations concerning dynamic effects. Ó 2015 Published by Elsevier B.V.

1. Introduction Accurate knowledge of the soil water retention curve (SWRC) and the hydraulic conductivity curve (HCC) is mandatory for modeling variably-saturated flow in porous media. The SWRC relates the water content of a soil to its pressure head at the scale of some representative elementary volume and is considered to be a characteristic relationship of a porous medium which reflects its equivalent pore-size distribution (Kutilek and Nielsen, 1994). It is usually assumed that a SWRC determined under static conditions can be used for modeling transient flow processes in the soil and vice versa. However, various experimental and theoretical studies suggest that there is a discrepancy between the SWRC estimated under static and transient conditions. This phenomenon has been associated with a dynamic non-equilibrium (DNE) between water

⇑ Corresponding author.

content and pressure head. Hassanizadeh et al. (2002) and Diamantopoulos and Durner (2012) give comprehensive reviews of experimental observations, physical reasons and approaches to model non-equilibrium flow in porous media. As pointed out by Hassanizadeh et al. (2002), early experiments in which DNE was observed were mainly focused on testing the validity of the diffusion equation. Rawlins and Gardner (1963) tested the validity of the diffusion equation for unsaturated flow of soil water by conducting absorption experiments in horizontal soil columns. They derived the diffusivity vs water content relationship D(h), found that it was not unique for the same porous medium, and concluded that this non-uniqueness could actually have resulted from either a non-uniqueness of SWRC, or HCC, or both. Topp et al. (1967) compared water retention curves which were estimated by equilibrium, steady-state, and transient methods for an identical soil column and found that for drainage more water was withheld by the soil matrix in the case of transient

E-mail address: [email protected] (E. Diamantopoulos). http://dx.doi.org/10.1016/j.jhydrol.2015.07.032 0022-1694/Ó 2015 Published by Elsevier B.V.

Please cite this article in press as: Diamantopoulos, E., et al. Modeling dynamic non-equilibrium water flow observations under various boundary conditions. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.07.032

2

E. Diamantopoulos et al. / Journal of Hydrology xxx (2015) xxx–xxx

experiments compared with the static or steady-state experiments. Thus, Topp et al. (1967) have shown for the first time that the retention curve depends on the flow conditions. Because the experiments were conducted for a monotonous drainage branch, the observed non-uniqueness was not attributed to the effect of hysteresis in the SWRC. Similar results were also presented by Schultze et al. (1999) who investigated DNE effects in transient multistep-outflow (MSO) and continuous-outflow (CO) experiments conducted to identify soil hydraulic properties (SHP). In a MSO experiment, an initially water-saturated soil column is drained by a stepwise decrease in pressure at the bottom. The SHP are determined by monitoring the outflow across the bottom and the pressure head in the soil and solving an inverse problem using the Richards equation (RE). In a CO experiment, the pressure at the bottom is decreased continuously and both water content and pressure head within the soil are monitored using two pairs of tensiometers and TDR probes. Schultze et al. (1999) observed that in the case of the MSO experiments the pressure head in the soil column reached hydrostatic equilibrium relatively quickly after a change in pressure at the bottom while outflow from the column continued for hours. In the CO experiments, the quicker the change in the lower boundary condition, the more water was withheld in the soil at the same pressure head thus confirming the observations of Topp et al. (1967). The observations in MSO and CO experiments indicate a DNE between water content and pressure head and, if the data are evaluated, lead to a non-uniqueness in the SWRC. These results suggest that a SWRC reflecting equilibrium conditions cannot necessarily be determined uniquely from transient experiments. Poulovassilis (1974) tested the uniqueness of the water retention curve by performing a different series of experiments. He percolated water at a constant rate through soil columns which were packed with a sandy material and equipped with five tensiometers. After achieving unit-gradient conditions, the column was sealed on both ends and left at rest in a horizontal position. Poulovassilis (1974) observed a relaxation to higher pressure heads at all tensiometers which confirms the findings of Topp et al. (1967). In a second series of experiments, Poulovassilis (1974) continued a vertical steady-state flow for 3 h after macroscopic steady-state conditions had been established. In this case, the relaxation of pressure head was not observable after bringing the column to a horizontal position. These two series of experiments indicate that water was redistributed even during steady state flow. Weller et al. (2011) and Weller and Vogel (2012) presented the Multistep Flux (MSF) method for the direct measurement of the HCC. In the MSF method, the infiltration at the upper boundary of a soil column is controlled using a sprinkler and the pressure at the bottom is regulated to a value which corresponds to a pressure head measured near the surface of the column. The objective is to achieve unit-gradient conditions for which the water flux density equals the hydraulic conductivity. The authors observed that immediately after lowering the flux at the top the pressure head in the soil decreased. However, shortly afterward, the pressure head increased again while the water content of the column remained almost constant. This undershoot in pressure head contradicts the equilibrium assumption between water content and pressure head. It resembles the observations of Poulovassilis (1974) and has been confirmed very recently in the study of Kumahor et al. (2005). To sum up, we can distinguish between four different observations of dynamic NE. The first one is the non-uniqueness of the diffusivity vs water content function during imbibition as presented by Rawlins and Gardner (1963). The second observation is the flow-rate dependence of the water retention curve as observed by Topp et al. (1967) and in the CO experiments of Schultze et al. (1999). The third observation is the dynamic NE observed in

MSO experiments like those of Schultze et al. (1999) and the last observation is the dynamic NE observed in the MSF experiments of Weller et al. (2011) and Weller and Vogel (2012). It seems that in the case of experiments with controlled pressure head boundary conditions dynamic NE appears as a relaxation in the system-averaged water content while the pressure head in the soil column indicates hydrostatic equilibrium. For experiments with flux boundary conditions like the MSF method or the experiments of Poulovassilis (1974), DNE appears as a relaxation of the pressure head while the flux density and macroscopic water content distribution appear static. In this article we investigate whether it is possible to describe all the four above mentioned expressions of DNE with a variably-saturated flow model extended to account for DNE between local water content and pressure head. Diamantopoulos et al. (2012) presented a simple model for describing DNE in the case of MSO experiments which we refer to as Dual-NE. It assumes two fractions of water in the porous medium, one fraction in instantaneous equilibrium with the local pressure head, and the second fraction for which the equilibration of water content is time dependent. So far, we have only demonstrated that it is possible to describe DNE effects in MSO experiments (Diamantopoulos et al., 2012) with the Dual-NE model. In this article, we show that the model is able to describe the other expressions of DNE discussed above. To the best of our knowledge there is no model which has been tested against various expressions of DNE flow so far. We provide results of simulating the experiments leading to DNE under different boundary conditions and provide a thorough discussion of the experimental and modeling considerations concerning DNE water flow in unsaturated soils. 2. Theory 2.1. Richards equation and the nonlinear diffusion equation for water flow One-dimensional variably-saturated water flow is described by the 1D-Richards equation (Richards, 1931)

   @h @ @h K ¼ þ cosðuÞ @t @x @x

ð1Þ

where h [L3 L3] is the water content, t [T] is time, x [L] is the spatial coordinate, positive upwards, h [L] is the pressure head, K [L T1] is the unsaturated hydraulic conductivity and u [–] is 0° for vertical flow and 90° for horizontal flow. For horizontal flow, the RE can be written as a nonlinear diffusion equation

  @h @ @h DðhÞ ¼ @t @x @x

ð2Þ

where D(h) [L2 T1] is the diffusivity as a function of the water content defined as:

DðhÞ ¼ KðhÞ

@h @h

ð3Þ

2.2. Dual-fraction non-equilibrium model (Dual-NE) The Dual-NE model of Diamantopoulos et al. (2012) assumes two fractions of water in the porous medium, one fraction in instantaneous equilibrium with the local pressure head according to a given water retention curve, and the second fraction (fne) for which the equilibration of the water content is reached only slowly after a change in local pressure and thus is treated as being time-dependent. Water flow for the equilibrium fraction is described by the RE (1) and water flow for the non-equilibrium

Please cite this article in press as: Diamantopoulos, E., et al. Modeling dynamic non-equilibrium water flow observations under various boundary conditions. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.07.032

E. Diamantopoulos et al. / Journal of Hydrology xxx (2015) xxx–xxx

fraction is modeled using the approach of Ross and Smettem (2000). If we assume that the equilibration of the pressure heads between the two regions is instantaneous (at least relative to the movement of water in the main flow direction), the main equation describing water flow becomes (Diamantopoulos et al., 2012):

   @heq @hne @ @h ð1  f ne Þ þ f ne ¼ þ cosðuÞ K @z @t @t @t

ð4Þ

where heq [L3 L3] is the water content in the equilibrium region and hne[L3L3] is the water content in the non-equilibrium region. Ross and Smettem (2000) assume that h and h are decoupled and that the term @h@tne of Eq. (4) can be defined through the differential equation @h@tne ffi gðhne ; heq Þ where gðhne ; heq Þ is a known function of the actual and equilibrium water contents. Please note that the heq represents the equilibrium water content given by the (hydrostatic) water retention curve. Ross and Smettem (2000) assumed that

3. Application 3.1. Diffusivity vs water content relationship Eq. (2) is valid for describing horizontal unsaturated water flow only if the function D(h) is a unique function for a monotonous imbibition or drainage. Early studies were focused on testing whether this assumption holds. Rawlins and Gardner (1963) performed simple horizontal water absorption experiments by monitoring the evolution of the water content inside the soil column by using a gamma ray densitometer. Their experimental setup consisted of a 62 cm long, horizontal rectangular soil column and water was applied at a pressure head of 5.3 cm. From the measured water content vs distance profiles they calculated diffusivity vs water content function based on the method of Bruce and Klute (1956) who proposed to calculate D(h) from the soil misture vs distance profiles h(x) for a time t according to:

ðh h Þ

gðhne ; heq Þ ¼ eq s ne , where s describes the equilibration time of water content in the non-equilibrium region. 2.3. Parameterization of hydraulic properties The numerical solution of Eqs. (1), (2) and (4) requires to parameterize the SHP. We used the van Genuchten–Mualem (VGM) model (van Genuchten, 1980) for all simulations. The SWRC h(h) and HCC K(h) are given by the equations

( hðhÞ ¼

Se ¼

 m hr þ ðhs  hr Þ  1 þ jahjn ; h<0 hs ;

h  hr hs  hr

h  m i2 KðSe Þ ¼ K s  Sle  1  1  Se1=m

hP0

ð5Þ

ð6Þ

ð7Þ

where hs and hr [L3 L3] are saturated and residual water contents, respectively, a [L1], n [–], m [–] and l [–] are shape parameters, m ¼ 1  1n ; n > 1, and Se [–] is effective saturation. The hydraulic conductivity was calculated as a function of weighted average of the water content in the equilibrium and non-equilibrium region given by (1  fne)heq + fnehne.

3

Dðhj Þ ¼

  Z hj 1 dx  x dh 2t dh hj hi

ð9Þ

where hi is the initial water content and hj is the water content at which the value of D is calculated for time equal t. Fig. 1a shows the D(h) relationships estimated by Rawlins and Gardner (1963). We assume that, as diffusivity D decreases at any given water content as time t increases, the water content vs pressure head relationship is not unique. We performed numerical simulations for the absorption of water into a horizontal soil column with the RE and the Dual-NE model which reflect the experimental design of the study of Rawlins and Gardner (1963). We assumed a horizontal soil column 100 cm long filled with a sandy material. Initially the pressure head was equal to 750 cm and at t = 0 we applied a constant head boundary condition equal to 20 cm at one end of our simulation domain. For the non-equilibrium simulations we assumed that the equilibration time of water content is 2 h (s = 2h) and the fraction of the non-equilibrium site is 0.3 (fne = 0.3). These values correspond to values often estimated from laboratory experiments (Diamantopoulos et al., 2012). The D(h) function was calculated from the x vs t profiles for five different times (12, 30, 60, 120 and 180 min) according to Eq. (9). Results are shown in Fig. 1b. As the RE assumes instantaneous local equilibrium between water content and pressure head, it predicts that the diffusivity vs water content relationship is independent of the time of calculation.

2.4. Inverse modeling The data from the different laboratory experiments treated in this article were analyzed by inverse modeling. The objective function to be minimized for determining the vector of unknown p is given by the weighted-least-squares formulation parameters ~

Oð~ pÞ ¼

N X 2 wi ri ð~ pÞ

ð8Þ

i¼1

where N is the number of data points in the objective function, ri are the residuals, i.e. the differences between the observed and the model-predicted data, and wi are weights which reflect the reliability of the individual measurements. Iterative minimization of Eq. (8) with respect to the parameter vector ~ p was achieved with the SCE-UA global search scheme of Duan et al. (1992). Numerical simulations with both the RE and the Dual-NE model were done with the Hydrus-1D code (Šimu˚nek et al., 2008) which was modified according to Diamantopoulos et al. (2012). The details of the simulations including initial and boundary conditions are described in Section 3.

Fig. 1a. Diffusivity vs water content relationship estimated from the experimental water content vs distance profiles of a horizontal adsorption experiment (Rawlins and Gardner, 1963).

Please cite this article in press as: Diamantopoulos, E., et al. Modeling dynamic non-equilibrium water flow observations under various boundary conditions. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.07.032

4

E. Diamantopoulos et al. / Journal of Hydrology xxx (2015) xxx–xxx

Fig. 1b. Diffusivity vs water content relationship estimated from the simulated water content vs distance profiles for a hypothetical horizontal adsorption experiment. Simulations were performed with both the Richards equation and the Dual-NE model.

Fig. 2a. Measured and fitted outflow and pressure head data for a slow drainage experiment (Schultze et al., 1999). The applied pressure head at the bottom of the soil column changed smoothly from 16.7 cm to 60 cm in 192 h. The Richards equation was used for fitting the data.

Conversely, the Dual-NE model predicts that for a time equal to 12 min the diffusivity D is smaller than predicted by the RE because of the lower amount of water which has infiltrated into the soil column. This is the case because only a fraction of the water is in instantaneous equilibrium. As time increases both the second and third term of Eq. (9) increase but they resulted in lower value of D for the same water content. This is exactly what was observed by Rawlins and Gardner (1963) and can be explained by the fact that water moves through the more conducting pores of the sandy soil at the start of infiltration. Afterward, a redistribution of water occurs behind the wetting front which results in the increase of the water content for the same pressure head value. 3.2. Flow-rate dependence of the SWRC The second observation of DNE in water flow is expressed as a flow-rate dependence of the SWRC (Topp et al., 1967). Schultze et al. (1999) presented experiments which were explicitly designed for investigating this flow-rate dependence. They used undisturbed soil columns with a volume of 1 L filled with a sandy soil. The soil column was equipped with two pairs of tensiometers and TDR sensors in different depths. They conducted a series of drainage-imbibition cycles by changing the applied pressure head at the lower end of the soil column in a smooth way. The slowest drainage experiment had a duration of 192 h and the fastest drainage experiment had a duration of 3 h. The in-situ measured SWRC obtained from assigning the in situ water contents to the in situ pressure head data for the slow and the fast experiments can be seen in Fig. 2b. It becomes obvious that more water was retained by the soil matrix for the faster experiment in comparison to the slower one leading to a non-uniqueness of the SWRC. Firstly, we fitted the RE to the slowest experiment. This implies to assume that under these experimental conditions the SWRC determined by inverse modeling represents an equilibrium relationship. We estimated the SHPs by minimizing Eq. (8) which included the pressure head and cumulative outflow data. The results are depicted in Fig. 2a. The RE describes both cumulative outflow and pressure head data very well and the estimated SWRC passes very well through the in-situ SWRC observed during the slow experiment. However, it does not match the SWRC of the fastest experiment (Fig. 2b). In a second simulation, we treated the SHPs estimated from the first inverse simulation as constants and fitted the Dual-NE model to the data of the fastest experiment by optimizing only the non-equilibrium parameters s and fne. The fitted values were 0.2 ± 0.8 h for s and 0.87 ± 0.25 for fne, respectively. The results showed that the Dual-NE model can describe

Fig. 2b. In-situ measured water retention curves for a slow (192 h) and a fast (3 h) drainage experiment conducted for the same soil column. The SHPs estimated from the slow experiment (blue) cannot predict the measured in-situ water retention curve for the fast experiment (duration 3 h, red). However, using the estimated equilibrium SHP (blue) and fitting the fast experiment with the Dual-NE model with two non-equilibrium parameters, we were able to describe the dynamic retention curve for the fast drainage experiment. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the observed in-situ SWRC data by using the equilibrium SWRC determined from the slowest experiment as input to the Dual-NE model. To illustrate this more clearly, we conducted numerical simulations using the RE and the Dual-NE model for a hypothetical experiment similar to those of Schultze et al. (1999). We assumed an initially saturated soil column, 10 cm in length, packed with a sandy soil. As initial condition we assumed a hydrostatic pressure head distribution with a head equal to 0 cm at the lower boundary. In the first case (Fig. 3a), the pressure head at the lower boundary was smoothly changed from 0 cm to 80 cm in 500 h, then remained at 80 cm for 20 h and was then raised again to 0 cm within another 500 h. In the second case (Fig. 3b), we modified the speed of change of the lower boundary condition: it changed from 0 cm to 80 cm in 20 h, stayed equal to 80 cm for 60 h and was then raised again to 0 cm within another 20 h. Fig. 3a depicts the simulated ‘‘in-situ’’ SWRCs for the RE and for the Dual-NE model with different values of the non-equilibrium parameter s for the slow experiment. The results show that the ‘‘in-situ’’ curves are all identical for the four cases and equal to the equilibrium curve used for the simulations. However, in the case of the fast experiment, we observe a deviation from equilibrium (‘‘apparent hysteresis’’) of the in situ SWRCs. For the simulations with the RE the drainage and imbibition curves are identical since capillary hysteresis was not accounted for. However in the

Please cite this article in press as: Diamantopoulos, E., et al. Modeling dynamic non-equilibrium water flow observations under various boundary conditions. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.07.032

E. Diamantopoulos et al. / Journal of Hydrology xxx (2015) xxx–xxx

5

Fig. 4. This figure shows measured and fitted outflow and pressure head data for a MSO experiment for a loamy sand soil (Diamantopoulos et al., 2012). The experimental data show that after each pressure change at the lower boundary, the equilibration of the pressure head data in the column is quicker than that of the cumulative outflow data. The blue line illustrates that this cannot be described by the Richards equation since it requests simultaneous equilibration between pressure head and water content (or cumulative outflow). Accordingly, any attempt to describe the data must fail. On the contrary the Dual-NE model describes very well both the pressure head and outflow data. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Forward simulations were conducted for two basic drainage-imbibition experiments. In the first case (top) the pressure head at the lower boundary was changed smoothly from 0 to 80 cm in 500 h. After an equilibration time of 24 h it changed from 80 cm to 0 cm until 1000 h. In the second case (bottom) the pressure head at the lower boundary was smoothly changed from 0 to 80 cm in 20 h, remained there for 60 h, and then was raised again to 0 within 20 h. Forward simulations were conducted with the Richards equation and the Dual-NE model for various combinations of non-equilibrium parameters. The results show that in the first case (slow experiment, top) the in-situ measured water retention curves are identical for all models. Contrary to that, the fast drainage-imbibition experiment leads to an apparent hysteresis for the Dual-NE simulations. Specifically, in the case of drainage more water is withheld by the soil whereas in the case of imbibition less water is stored in the soil. The Richards equation predicts non-hysteretic flow as expected since water content and pressure head are always in equilibrium.

case of Dual-NE simulations, the water content is greater for the same pressure head as s increases in the case of drainage and smaller in the case of imbibition. This indicates that part of the hysteresis observed in dynamic laboratory experiments can be partly attributed to DNE effects.

remaining discrepancies between the observed and simulated outflow data can be attributed to error in the parameterization of the SHP or a possible water-content dependence of the rate parameter s which is not analyzed further in this article. The fitted non-equilibrium parameters were 3.2 ± 0.6 h for s and 0.78 ± 0.04 for fne, respectively.

3.4. MSF experiments Weller et al. (2011) presented a method for the direct measurement of the unsaturated hydraulic conductivity using a unit-gradient percolation experiment in which a water flux is applied at the top and a suction is applied at the bottom of a soil column. The pressure head at the bottom is adjusted to a value corresponding to the pressure head measured in the soil. The authors named this experimental technique the Multistep-Flux method (MSF) since it applies a sequence of different water fluxes at the top. The idea behind the MSF method is that after achieving steady-state conditions and a unit gradient in the hydraulic potential, the hydraulic conductivity equals the flux density.

3.3. MSO experiment MSO experiments are nowadays a standard method for the determination of SHP in the laboratory (Hopmans et al., 2002). Fig. 4 presents experimental results of an MSO experiment conducted on an undisturbed loamy sand (Diamantopoulos et al., 2012). The applied pressure heads at the lower boundary can be seen in Fig. 4. After a sudden change in pressure head at the lower boundary, the pressure head in the soil equilibrates relatively quickly and reaches a value reflecting hydrostatic conditions whereas the outflow dynamics are separated into two stages. Immediately after the pressure change at the bottom, a relatively large amount of water leaves the soil sample quickly. This first phase is followed by a relatively slow drainage of water during which the pressure head in the soil does not change any more. The RE cannot describe these dynamics because it assumes instantaneous local equilibrium between water content and pressure head. As shown in Fig. 4, the Dual-NE model describes both the pressure head and the cumulative outflow data very well and, more specifically, is able to describe a continuous outflow even for a static pressure head distribution in the soil column. The

Fig. 5. Experimental and simulated data from the study of Weller et al. (2011). After a stepwise change of the applied flux at the upper boundary, the pressure head drops as expected (drainage) but then it rises again to a higher value. Richards’ equation predicts that after the establishment of steady-state flux with constant water content in the soil, the pressure head remains constant as well and corresponds to the one obtained from the retention curve. This contradicts the experimental findings. In contrast, the Dual-NE model predicts the observed slow equilibration of the pressure head data even for constant water content.

Please cite this article in press as: Diamantopoulos, E., et al. Modeling dynamic non-equilibrium water flow observations under various boundary conditions. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.07.032

6

E. Diamantopoulos et al. / Journal of Hydrology xxx (2015) xxx–xxx

Fig. 5 presents experimental results for an MSF experiment conducted on a 10 cm long undisturbed sandy soil column. The experimental data are taken from Weller et al. (2011). The initially applied flux density was equal to 1.4 cm/h and it was subsequently reduced by a factor of two in each step thus causing a decrease in water content. The experimental results showed that after an applied flux at the upper boundary the pressure head drops quickly as is expected for a drainage experiment. However, the pressure head increases afterward although the column-averaged water content remains constant with time. One possible explanation is that after the reduction of the flux at the upper boundary, a part of the larger pores drains easily (equilibrium domain) and the water content and the pressure head adapt instantaneously to the applied flux. This corresponds to the phase of the pressure head undershoot. Afterward, water from the non-equilibrium domain becomes slowly connected to the main flow because of a limited pore continuity. This water is redistributed toward smaller pores which leads to the macroscopic increase in the pressure head observed in the experiment. Fig. 5 shows the results of fitting a numerical solution of the RE to the experimental data. The RE cannot describe the observed pressure head undershoot since it assumes that under unit-gradient conditions, the pressure head and the water content should be in equilibrium. In contrast to the RE, the Dual-NE model describes the observed pressure head data very well. 4. Discussion Macroscopically, all observations presented above can be described by the following mechanism. Even in well sorted materials there are some highly conducting pore paths in which the flow is initiated after the boundary condition is changed. Behind the wetting or drying front, however, there is a second phase-equilibration of a relative small amount of water between pores of different sizes within the REV. This second phase of equilibration is not instantaneous and it is responsible for the time-varying and thus non-unique behavior of the h vs h relationship at the REV scale. This explains why the moisture vs horizontal distance profiles change their shape differently than it is predicted by the diffusion theory and leads to the non-uniqueness of the diffusivity vs water content function (Section 3.1). It can also explain the rate-dependence of the SHP (Section 3.2), the two-stage dynamics of the outflow in the MSO experiments (Section 3.3) and finally the pressure head undershoot in the MSF experiments (Section 3.4). Although the Dual-NE model can describe all these different expressions of DNE, we admit that the exact pore-scale processes which are responsible for these macroscopic effects cannot be studied by macroscopic experiments, measurements and models. These processes have until now not been well identified and DNE in variably-saturated water content has been associated with (i) air–water interface equilibration (Barenblatt, 1971; Sakaki et al., 2010), (ii) water entrapment and pore water blockage (Topp et al., 1967), (iii) air-entry effects (Wildenschild et al., 2001; Diamantopoulos and Durner, 2012), (iv) gravity flow, (v) air entrapment (Schultze et al., 1999), (vi) dynamic contact angle (Friedmann, 1999) and (vii) fluid properties (Das et al., 2007; Goel and O’Carroll, 2011) among others. Given the various physical explanations for DNE in variably-saturated flow it appears questionable that a physically-based model accounting for only a single cause is able to describe experimental observations of DNE (Diamantopoulos et al., 2012; Diamantopoulos and Durner, 2012). In soil physics, the model most often applied to DNE data is the dynamic-capillary pressure model of Hassanizadeh and Gray (1993) which was derived on a thermodynamic basis. The model states that the dynamic capillary pressure in a two-phase flow system is different from the

capillary pressure under static conditions and that both are related by a function accounting for the temporal change of saturation. If the air pressure in the porous medium is equal to the atmospheric pressure, this function reads

hdyn ¼ hstat þ

s @h eqg @t

ð10Þ

where hdyn and hstat are the pressure heads measured under dynamic or static conditions [L], respectively, s [M L1 T1] is a material coefficient which may depend on saturation (Hassanizadeh et al., 2002), on texture and on sample size, e is the porosity [L1 L1], q is density of water [M L3]and g is the acceleration of gravity [L T2]. Barenblatt (1971) and Barenblatt and Gil’man (1987) presented a model for describing dynamic effects in two-phase flow. They distinguish between an effective water saturation and an actual saturation S [–] which can be measured in a sample and which tends toward the effective saturation during equilibration (Silin and Patzek, 2004). This model is very similar to the one of Ross and Smettem (2000), because both models assume that water saturation (or water content) changes according to a first-order kinetics. Hassanizadeh et al. (2002) point out that the identification of the equilibration time in the Barenblatt model is difficult because the model should be viewed as an empirical model which cannot be derived from first principles. The same holds for the model of Ross and Smettem (2000). Above we have demonstrated that the Dual-NE model is capable of describing all effects of dynamic nonequilibrium treated in this article. This is not the case for the other models discussed above and will be explained now. For the model of Hassanizadeh and Gray (1993), the dynamic pressure head in Eq. (10) equals the static one after achieving steady-state flow conditions. Therefore, the relaxation of pressure head observed in the MSF experiments and the experiments of Poulovassilis (1974) cannot be described by this model. Another consequence of the model of Hassanizadeh and Gray (1993), is that soil hydraulic properties estimated from steady-state experiments should be identical to those determined from static experiments. This was experimentally confirmed by Topp et al. (1967) who found that SWRCs determined by steady-state and static experiments were almost identical in the case of drainage. However, Topp et al. (1967) report that ‘‘the time to establish steady state varied from 300 to 600 min’’ for the coarse sand used in their study (particle-size distribution: 1–0.5 mm, 0.1%; 0.5–0.25 mm, 28.8%; 0.25–0.10 mm, 68.7%; < 0.1 mm, 2.4%). This implies that waiting times of 5–10 h to reach steady-state flow in this coarse medium may have masked most of the second-stage equilibration and effectively hidden the effects of DNE on the estimated SWRC. For the MSF experiments, the model of Ross and Smettem (2000) predicts a decrease in pressure head which is much stronger than observed when the flux is changed at the top. The reason for this is that immediately after changing the flux, all water is assumed to be in the non-equilibrium domain. This results in a very small effective water capacity which causes a strong decrease in pressure head. In contrast to this, the Dual-NE model assumes that some part of the water is in the equilibrium domain after the flux changes. This dampens the decrease in pressure head while sustaining the relaxation of the pressure head due to the presence of water in the non-equilibrium domain. For this reason we believe that no single domain non-equilibrium model can describe the pressure head dynamics observed in the steady-state MSF experiments discussed in Section 3.4. To sum up, single-domain non-equilibrium models like those of Hassanizadeh and Gray (1993), Barenblatt (1971) and Ross and Smettem (2000) can describe some of the often observed DNE effects in soil physics, but not all of them. At least some

Please cite this article in press as: Diamantopoulos, E., et al. Modeling dynamic non-equilibrium water flow observations under various boundary conditions. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.07.032

E. Diamantopoulos et al. / Journal of Hydrology xxx (2015) xxx–xxx

7

slow flow conditions is most appropriate for modeling water flow in the field. If one needs to conduct a quicker experiment, it should be analyzed with a flow model accounting for non-equilibrium. Only then is it possible to determine the correct SWRC which can be used for predicting water flow in the field with the RE. However, non-equilibrium water flow can occur in cases where water flow is fast such as in the case of rapid infiltration conditions or irrigation. 5. Conclusions

Fig. 6. Measured and fitted outflow and pressure head data for a MSO experiment for the same soil column as in Fig. 5. The experimental data show that after each pressure change at the lower boundary, the equilibration of the pressure head data in the column is quicker than that of the cumulative outflow data. The blue line illustrates that this cannot be described by the Richards equation. Fitting the DualNE model describes the experimental results very well, but with a significantly lower value of the s parameter (s = 3h, red dashed line). Keeping s parameter constant to the value fitted from the MSF experiment (s = 32.2h, drawn red line) gives a slightly reduced quality of the fit but describes the experimental data still reasonably well. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

phenomena, require dual-domain approaches like the Dual-NE model which describe an equilibration of water in two stages. An important issue in modeling DNE is to test whether multiple experiments causing different flow dynamics lead to similar estimates of equilibrium and non-equilibrium parameters. This will answer the important question whether DNE parameters are a function of physical or physicochemical soil properties like surface properties (e.g. wettability) or soil structural features. We are aware only of the study of Weller et al. (2011) in which MSO and MSF experiments were conducted on the same soil column. The results of the MSF experiments were already discussed in Section 3.4 and presented in Fig. 5. The estimated parameters for the non-equilibrium model were 30.2 ± 2.8 h for parameter s and 0.36 ± 0.01 for parameter fne. Fig. 6 shows the fitted RE and Dual-NE model to the MSO experimental results for the same soil column. Again, the RE cannot describe the second phase equilibration of the cumulative outflow data contrary to the fit obtained with the Dual-NE model. The fitted non-equilibrium parameters were equal to 3.3 ± 0.3 h for s and 0.31 ± 0.01 for fne. In a second step, we fitted again the DNE model but kept the parameter s fixed to 30.2 h as estimated from the MSF experiment. The results can be seen in Fig. 6 which shows that the Dual-NE model describes very well the experimental MSO results and most importantly with identical SHPs and the same values of the non-equilibrium parameters (fitted fne was equal to 0.35 ± 0.02). This indicates the possibility to describe both MSO and MSF experiments with the same set of equilibrium and non-equilibrium parameters. We have to notice that it is more likely that the non-equilibrium parameters should not stay constant but they probable be a function of saturation or of the pressure head. Additionally, the non-equilibrium parameters can be estimated from structural analysis at the pore scale. However, more experimental work is needed on this direction. Obviously, DNE is present in almost all soil physical laboratory experiments. The major difference between water dynamics in the field and water dynamics in laboratory experiments are that (i) stepwise changes in pressure head boundary conditions do not occur in the field and that (ii) laboratory experiments often accelerate flow processes to enable the determination of SHPs because of time constraints. To illustrate the problem, consider the two retention curves estimated by Schultze et al. (1999) which are presented in Fig. 2b. One may assume that the SWRC estimated under

We have applied the Dual-NE model of Diamantopoulos et al. (2012) to various experimental observations of DNE in variably-saturated flow problems. The model conceptualizes that the water content is divided into two fractions. In the first fraction the equilibration of the water content and the pressure head is instantaneous whereas in the second fraction this equilibration is slow and is therefore described by a first-order kinetics which leads to a temporal shift in water content even if the pressure head is constant. The model requires two more parameters than the RE: one describes the partitioning of the total water into the fractions and the second parameter describes the time-scale of the equilibration of the water content. This conceptual approach reflects our understanding that due to the complex structure at the pore scale a part of the pore volume can easily adapt its saturation in response to external forcing while another part cannot. The Dual-NE approach is a first order approximation of this general phenomenon. We have demonstrated that the Dual-NE model describes various DNE effects reported in the literature very well: (i) the non-uniqueness of the diffusivity vs water content function during absorption experiments, (ii) the flow-rate dependence of estimated SHPs, (iii) the two stage equilibration of cumulative outflow data in MSO experiments and (iv) the pressure head undershoot which has been observed during steady-state MSF experiments. Future research will focus on DNE water flow under different boundary conditions and on whether different experiments conducted on the same soil sample yield similar SHPs and non-equilibrium parameters. Finally, it should be analyzed whether DNE parameters can be related to soil textural and structural properties and surface properties like wettability and contact angle. Acknowledgments This study was financially supported within the DFG Research Group FOR 1083 ‘‘Multi-Scale Interfaces in Unsaturated Soil’’ (MUSIS), Grant DU283/10-1. References Barenblatt, G.I., 1971. Filtration of two non-mixing fluids in a homogeneous porous medium. Izvestiya (Bulletin), USSR Academy of Sciences, Mekhanika Zhidkosti i Gasa (Mechanics of Fluid and Gas), No. 5, pp. 144–151 (in Russian); Fluid Dynamics, 6, No. 5, pp. 857–864 (English Translation). Barenblatt, G.I., Gil’man, A.A., 1987. Nonequilibrium counterflow capillary impregnation. J. Eng. Phys. 52 (3), 335–339. Bruce, R.R., Klute, A., 1956. The measurement of soil moisture diffusivity. Soil Sci. Soc. Am. Proc. 20, 458–462. Das, B.D., Gauldie, R., Mirzaei, M., 2007. Dynamic effects for two-phase flow in porous media: fluid property effects. AIChE J. 53 (10), 2505–2520. Diamantopoulos, E., Durner, W., 2012. Dynamic nonequilibrium of water flow in porous media: a review. Vadose Zone J. 11 (3), (verified 14.02.14). Diamantopoulos, E., Iden, S.C., Durner, W, 2012. Inverse modeling of dynamic nonequilibrium in water flow with an effective approach. Water Resour. Res. 48 (3), n/a–n/a. http://dx.doi.org/10.1029/2011WR010717 (verified 14.02.14). Duan, Q., Sorooshian, S., Gupta, V., 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour. Res. 28, 1015–1031. Friedmann, S., 1999. Dynamic contact angle explanation of flow rate-dependent saturation-pressure relationships during transient liquid flow in unsaturated porous media. J. Adhes. Sci. Technol. 13 (12), 1495–1518.

Please cite this article in press as: Diamantopoulos, E., et al. Modeling dynamic non-equilibrium water flow observations under various boundary conditions. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.07.032

8

E. Diamantopoulos et al. / Journal of Hydrology xxx (2015) xxx–xxx

Goel, G., O’Carroll, D.M., 2011. Experimental investigation of nonequilibrium capillarity effects: Fluid viscosity effects. Water Resour. Res. 47 (9), n/a–n/a. http://dx.doi.org/10.1029/2010WR009861 (verified 14.02.14). Hassanizadeh, S.M., Gray, W.G., 1993. Thermodynamic basis of capillary pressure in porous media. Water Resour. Res. 29 (10), 3389–3405. Hassanizadeh, S.M., Celia, M.A., Dahle, H.K., 2002. Dynamic effects in capillary pressure-saturation relationships and its impact on unsaturated flow. Vadose Zone J. 1, 38–57. Hopmans, J.W., Simunek, J., Romano, N., Durner, W., 2002. Simultaneous determination of water transmission and retention properties – inverse methods. In: Dane, J.H., Topp, G.C. (Eds.), Methods of Soil Analysis, Part 4. SSSA Book Ser. 5. SSSA, Madison, WI, pp. 963–1008. Kumahor, S., de Rooij, G., Schlüter, S., Vogel, H.-J., 2005. Water flow and solute transport in unsaturated sand—a comprehensive experimental approach. Vadose Zone J. 14. http://dx.doi.org/10.2136/vzj2014.08.0105. Kutilek, M., Nielsen, D., 1994. Soil Hydrology. Catena Verlag, Cremlingen-Destedt, Germany. Poulovassilis, A., 1974. The uniqueness of the moisture characteristics. J. Soil Sci. 25, 27–33. Rawlins, S.L., Gardner, W.H., 1963. A test of the validity of the diffusion equation for unsaturated flow of soil water. Soil Sci. Soc. Am. Proc. 27, 507–511. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics (College Park, MD). 1 (5), 318–333, (verified 22.01.14). Ross, P.J., Smettem, K.R.J., 2000. A simple treatment of physical nonequilibrium water flow in soils. Soil Sci. Soc. Am. J. 64, 1926–1930.

Sakaki, T., O’Carroll, D.M., Illangasekare, T.H., 2010. Direct quantification of dynamic effects in capillary pressure for drainage-wetting cycles. Vadose Zone J. 9 (2), 424. Schultze, B., Ippisch, O., Huwe, B., Durner, W. 1999. Dynamic nonequilibrium during unsaturated water flow. In: Proc. Int. Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media, 1997, University of California, Riverside, CA. pp. 877–892. Silin, D., Patzek, T., 2004. On Barenblatt’s model of spontaneous countercurrent imbibition. Transp. Porous Media 54, 297–322. Šimu˚nek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of the HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7 (2), 587–600. Topp, G.C., Klute, A., Peters, D.B., 1967. Comparison of water content-pressure head data obtained by equilibrium, steady state, and unsteady-state methods. Soil Sci. Soc. Am. Proc. 31, 312–314. Van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898. Weller, U., Vogel, H.-J., 2012. Conductivity and hydraulic nonequilibrium across drainage and infiltration fronts. Vadose Zone J. 11 (3), (verified 12.06.14). Weller, U., Ippisch, O., Köhne, M., Vogel, H.-J, 2011. Direct measurement of unsaturated conductivity including hydraulic nonequilibrium and hysteresis. Vadose Zone J. 10 (2), 654, (verified 12.06.14). Wildenschild, D., Hopmans, J.W., Simunek, J., 2001. Flow rate dependence of soil hydraulic characteristics. Soil Sci. Soc. Am. J. 65, 35–48.

Please cite this article in press as: Diamantopoulos, E., et al. Modeling dynamic non-equilibrium water flow observations under various boundary conditions. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.07.032