A novel double disc method to determine soil hydraulic properties from drainage experiments with tension gradients

Journal of Hydrology 615 (2022) 128625 Available online 29 October 2022 0022-1694/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Research papers A novel double disc method to determine soil hydraulic properties from drainage experiments with tension gradients D. Moret-Fern´ andez * , B. Latorre Departamento de Suelo y Agua, Estaci´ on Experimental de Aula Dei, Consejo Superior de Investigaciones Científicas (CSIC), PO Box 13034, 50080 Zaragoza, Spain ARTICLE INFO This manuscript was handled by Corrado Corradini, Editor-in-Chief, with the assistance of Renato Morbidelli, Associate Editor Keywords: Water retention curve Hydraulic properties Inverse analysis Hydraulic conductivity Upward infiltration Constant head method ABSTRACT Determination of the saturated hydraulic conductivity, K s , and the water retention curve, θ(h), is of paramount importance to characterize the hydraulic behavior of the vadose zone. Given the van Genuchten hydraulic model, defined by the residual, θ r , and saturated, θ s , volumetric water content and the α and n parameters, this work presents a new laboratory procedure to estimate K s , θ s , n and α for a drainage process, α dr , from the inverse analysis of successive drainage steady-states curves generated by a tension-gradient between the surface and the base of a soil column. To this end, a double disc system, one connected a bubbling tower and placed at the soil surface and the second one placed under the soil core, was employed. The second disc was connected to an airvacuum system. The experiment presented two parts: a first 1D downward infiltration at saturation on a dry soil column, followed by successive drainage steps. During the drainage process, the tension of the upper and lower discs varied between 0 and −5 cm, and from −5 to −100 cm, respectively. The soil sorptivity, S, and θ s were calculated from the 1D transient infiltration measure, K s was calculated by Darcy’s law, α dr and n were optimized from the inverse analysis of the steady-state curves under tension-gradient and α for a wetting process, α w , was calculated from previously obtained S, θ s , K s and n. Once K s estimated, α dr and n were optimized by minimizing the Q= |hb−hn|objective function, where h b and h n are the experimental and calculated tensions at the base of the soil core. Given a α dr value, the optimimum n was computed as the value that provides a minimum Q. By repeating this process for a sequence of α dr , different Q-isolines were obtained, one for each h b value, which crossing-point corresponded to the actual α dr and n values. The method was tested on 2.5 cm high columns of four different synthetic soils. Next, it was applied on an experimental sand column of 5 cm height and on 2.5 cm high columns filled with sieved loam, clay loam and clay soil. The estimated α dr and n were compared with corresponding values measured in the same soils with the pressure plate technique and α w was contrasted with the corresponding value calculated with an empirical hysteresis model. The method, which was fast (from 1 to 2 h) and easy to implement for small-scale experiments, was successfully applied to soil samples 2.5 cm high and allowed to explore a range of soil tensions from 0 to −100 cm. Overall, accurate estimates of θ s , K s , α dr and n were obtained in both synthetic and experimental soils. A significant relationship was also obtained between α w estimated from S and the corresponding value calculated from the hysteresis model. 1. Introduction Accurate determination of soil hydraulic properties is of paramount importance for correct simulation of soil water flow in the vadose zone. The soil hydraulic properties are defined by the water retention and the hydraulic conductivity functions. The soil water retention curve, θ(h), expresses the relationship between the volumetric water content, θ [L 3 ⋅L −3 ], and the matric potential, h [L] under static conditions. The hydraulic conductivity is defined as the ability of the fluid to pass through the pores of the material (Newby et al. 2009). Its value decreases with the soil water content: as the soil moisture decreases, the total number of water conductive pathways along which fluid can travel reduces, and hence the K(h) value. Accordingly, the K(h) is a measure of the increased impedance to water flow with decreasing moisture content (McCartney et al., 2007). These functions are characteristic for different types of soil. van Genuchten (1980)-Mualem (1976) described one of the most employed models to characterize θ(h) and K(h). In this case, θ(h) is defined by the saturated (θ s ) and residual (θ r ) volumetric water contents, * Corresponding author. E-mail address: [email protected] (D. Moret-Fern´ andez). Contents lists available at ScienceDirect Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol https://doi.org/10.1016/j.jhydrol.2022.128625 Received 2 June 2022; Received in revised form 22 August 2022; Accepted 19 September 2022 Journal of Hydrology 615 (2022) 128625 2 the empirical α and n factors, and the m parameter, commonly defined as m=1−1 n. θ r expresses the water content for which dθ/dh becomes zero (excluding the region near θ s which also has a zero gradient), n [−], which is related to pore-size distribution, defines the slope of θ(h) and α [L −1 ] is a scale factor that describes the shape of θ(h) near θ s . K(h) can be defined as a function of α , n and the saturated hydraulic conductivity, K s , defined this last one as the soil ability to transmit water at saturation conditions. The soil water retention function is subject to the hysteresis phenomenon, manifested as a difference between equilibrium curves of soil wetting and drying (hysteresis loop) (e.g. Hillel, 1998). This means that water content in the drying (or drainage) branch of water potential is larger than water content in the wetting branch for the same value of water potential. The causes of hysteresis can be the variation in contact angle at different wetting and drying cycles, entrapped air in a newly wetted soil, temperature, swelling and shrinking, and ink-bottle effect due to nonuniformity in shape and sizes of both individual pores and interconnected pore networks (Bachmann and van der Ploeg, 2002; Maqsoud et al., 2004). Laboratory methods to determine the soil hysteresis are both time consuming and complicated. To overcome this limitation, different physical and empirical models for predicting one branch of the water retention relation have been developed (Pham et al., 2003; Gebrenegus and Ghzzehei, 2011). To date, two families of methods to characterize the soil hydraulic properties under laboratory conditions are available: those based on direct measures and those founded on the inverse numerical analysis of Richard’s water flows. The saturated hydraulic conductivity, K s , can be directly measured with either the constant head or the falling-head methods (Klute and Dirksen,1986). Among the direct techniques to measure the water retention curve, the pressure membrane apparatus (Richards, 1941), which is considered as the reference laboratory method, determines θ(h) from measured θ and h pairs. Although this technique has been progressively improved by introducing different ways to estimate θ (Jones et al., 2005; Moret-Fern´ andez et al., 2012), the long time required to complete a measurement together with errors of apparatus in fine-textured soils (Solone et al., 2012) may limit its use. Another direct procedure is the evaporation method (Wind, 1968; Wendroth et al., 1993, among others), where the soil water retention and hydraulic conductivity data are obtained through monitoring of the water content (evaporation) and matric potential dynamics at various heights of a soil sample exposed to evaporation. Although this method can extend the range of pressures up close to the wilting point (Schindler et al., 2010), it is time consuming, it has to use relatively long soil cores (6–11 cm) and results inaccurate at near saturation conditions (Sarkar et al., 2019). To overcome this last limitation, the evaporation method is commonly complemented with the double disc or quasi unit-gradient percolation procedure, which currently is the most straightforward way to measure unsaturated hydraulic conductivity near saturation (Dirksen, 1999). In this case, water is supplied at the top of a soil column and drained at the bottom by suction; this requires the pressure head at the bottom to be adjusted to a value corresponding to the pressure head measured in the soil near the top (Sarkar et al., 2019). However, although this method allows estimating the hydraulic conductivity at near saturation conditions, it is time consuming and the applied tension is limited to from 0 to −10 cm (Sarkar et al., 2019). Currently, four different 1D transient methods can be used to estimate soil hydraulic properties by inverse analysis: evaporation and horizontal-, downwardand upward-infiltration processes. The main advantage of these techniques is the simultaneous estimate of θ(h) and K (h) (Moret-Fern´ andez et al., 2016). In the horizontal infiltration method (Shao and Horton, 1998), K s is measured by Darcy’s law and α and n parameters are estimated with an integral method that simulates the problem of water absorption into a horizontal soil column. To this end, a 20 cm-length transparent cylinder should be used. The 1D downward infiltration method based on the inverse analysis of a single transient cumulative infiltration at saturated conditions allows determining the soil sorptivity, S, and K s (Lassabatere et al., 2009; Latorre et al., 2015; Moret-Fern´ andez et al., 2021b), but not the β parameter (Latorre et al., 2018) defined Haverkamp et al. (1994a, 1994b). S, which is defined as a measure of the capacity of a medium to absorb or desorb liquid by capillarity (Philip, 1957), can be expressed as function of the initial (θ i ) and saturated (θ s ) water content, K s and the α and n parameters of the van Genuchten (1980) water retention model (Moret-Fern´ andez et al., 2017). β is an integral shape parameter, which can be expressed as function of θ i , θ s and the shape n parameter of the van Genuchten (1980) water retention model (Moret-Fern´ andez and Latorre, 2017). However, given the inverse analysis of a 1D transient infiltration allows only estimating two (S and K s ) of the three parameters of the Haverkamp et al. (1994a, 1994b) model (Latorre et al., 2018), α and n cannot be estimated from S, and thus from a transient infiltration curve (Simunek and van Genuchten, 1996). Among the wide variety of existing procedures to estimate the soil hydraulic parameters from the inverse analysis of an upward infiltration curves (Hudson et al., 1996; Young et al., 2002; Moret-Fern´ andez et al., 2016; Pe˜ na-Sancho et al., 2017), the most complete one is that presented by Latorre and Moret-Fern´ andez (2019) and Moret-Fern´ andez et al. (2021a), which allows determining four of the five van Genuchten (1980) parameters from the inverse analysis of a single infiltration curve. However, although this method is accurate, can work with short 5 cm-high cores and run with disturbed and undisturbed soils, the inverse analysis is time consuming (from 2 to 20 days) and requires high-performance computing. On the other hand, since the inverse analysis of transient infiltrations measured in stratified soils profile leads to erroneous estimates of soil hydraulic parameters (Moret-Fern´ andez et al. 2021b), the inverse analysis of a 1D upward infiltration should be limited to homogeneous soil profiles. However, soil heterogeneity is the rule rather than the exception under field conditions (Abou Najm et al., 2019). For example, the application of the Sequential Inverse Analysis (SIA) (Moret-Fern´ andez et al., 2021b) method on 20 infiltration experiments showed that in 85 % of the soils analyzed, the thickness of the homogeneous topsoil layer was less than 3 cm . This problem associated with soil layering could be minimized, for instance, by using shorter soil columns or analyzing steady-state infiltration curves. In summary, although there is a large variety of laboratory methods to estimate the soil hydraulic properties, most of them presents some of the following constraints: (i) long measurement or analysis time; (ii) limitations in the applied tension; (iii) high computing requirements; (iv) poor behavior in stratified soils; or (v) use of long soil columns. The objective of this work is to present a new laboratory method and a numerical procedure to estimate the soil hydraulic properties (S, K s , θ s , α and n), attenuating most of the above-mentioned limitations. In the proposed method, the soil hydraulic properties for a drainage process are estimated from the inverse analysis of successive steady-states drainage curves generated by applying a tension gradient between the surface and the base of a soil cylinder. Based on the quasi unit-gradient percolation procedure, the method employs a double disc system (one placed at the surface and the other one at the base of the soil column), where tensions in the upper and lower discs range between [0, −5] and [-5, −100] cm, respectively. This procedure generates a first downward infiltration on an initially dry soil column followed by successive steadystate drainages curves of increasing rates. In this case, S and θ s are measured from the 1D transient infiltration curve at saturation conditions, K s is measured by Darcy’s law when tensions at the surface and bottom of the soil column are equal to 0 cm, and n and α for a drainage process, α dr are estimated form the inverse analysis of the drainage curves. Finally, α for a wetting process, α w , is calculated from previously obtained S, θ s , K s and n. The method was first validated on four different synthetic soils columns of 2.5 cm high. The successive steady-state drainage experiment was also employed to check how quickly the infiltration curve stabilizes after the suction at the bottom increases. Making use of a new design of double disc plus air-vacuum system, the D. Moret-Fern´ andez and B. Latorre Journal of Hydrology 615 (2022) 128625 3 method was applied on four sieved soils columns with different textures (from sand to clay) and the estimated α dr and n were compared with corresponding values measured with the pressure plate method. The estimated α w was compared to the corresponding value calculated from α dr using an empirical hysteresis model. 2. Theory 2.1. Soil hydraulic functions and cumulative infiltration model The governing flow equation for one-dimensional isothermal Darcian flow in a variably saturated rigid porous medium is given by the modified form of the Richards equation. dθ dt =d dz (Kdh dz +K)(1) where θ is the volumetric water content (L 3 L -3 ), h is the soil–water pressure head (L), K is the hydraulic conductivity (L T −1 ), z is a vertical coordinate (L) positive upward, and t is time (T). Under steady-state conditions, Eq. (1) simplifies to. q=d dz (Kdh dz +K)(2) where q is the water flux density [L T −1 ]. For saturated soils condition, Eq. (2) reduces to the Darcy’s law (Lichtner et al., 1996). q=Ks dz dH (3) where and K s is the saturated hydraulic conductivity [L T −1 ] and H =h +z is the total head. The soil water retention and hydraulic conductivity functions can be described with the van Genuchten (1980) - Mualem (1976) models according to. θ(h) = θr+ (θs−θr)[1 1+ ( α h)n]m (4) K(h) = Ks{1− ( α h)n−1[1+ ( α h)n]−m}2 [1+ ( α h)n]m 2 (5) where θ s and θ r are the saturated and residual volumetric water content [L 3 L −3 ], respectively, α [L −1 ] and n [−] are a scale and shape parameter, respectively, and m =1–1/n. According to van Genuchten (1980), θ r is defined as the water content for which the gradient dθ/dh becomes zero. The soil water retention function is subject to a hysteresis phenomenon, manifested as a difference between equilibrium curves of soil wetting and drying (hysteresis loop) (e.g. Hillel, 1998). Given that the soil hysteresis affects only the α value (Haverkamp et al. 2002), a scale parameter for wetting, α w , and drainage, α dr , processes should be defined. The α w parameter can be predicted from drying water retention data using the I H index (Gebrenegus and Ghezzehei, 2011), IH=(rn−1 rn(n−1) 2n−1−1)1 n−1 −(rn−1 rn−rn2 2n−1)1 n−1 (6) where r = α dr / α w . In the absence of measured wetting and drying water retention data, I H is calculated as. I H =0.378 ln(n) (7). The 1D cumulative infiltration, I1D, under saturation conditions can be described by the Haverkamp et al. (1994a, 1994b) model as. 2(Ks−Ki)2 S2t=2 1−β (Ks−Ki)(I1D−Kit) S2−1 1−βln[1 βexp(2β(Ks−Ki)(I1D −Kit)/S2)+β−1 β] (8) expression that is valid for the whole infiltration time. S is the sorptivity (L T −0.5 ), K i (L T −1 ) is the hydraulic conductivity value corresponding to the initial, θ i , volumetric water contents (L 3 L -3 ) and β is an integral shape parameter (-). The β parameter varies between 0.6 and 1.7 (Lassabatere et al., 2009), however, because β has a negligible influence on 1D infiltration (Latorre et al., 2018), a constant β is commonly considered (Angulo-Jaramillo et al., 2019). Given the complexity to solve the implicit Eq. (8), the cumulative 1D infiltration for saturated soil conditions and negligible initial hydraulic conductivity can be approximated by the explicit 4-terms, 4 T, expansion (MoretFern´ andez et al., 2020). I1D4T(t) = St1 2+(2−β 3Ks)t+K2 s 9S(β2−β+1)t3 2+2(β−2)(β1)(1−2β) 135 K3 s S2t2 (9) which is valid for long infiltration times (between 2.000 and >50.000 s for coarse and fine soil textures, respectively). The soil sorptivity expressed as a function of van Genuchten (1980) model results (MoretFern´ andez et al., 2017). S2=(1−m)Ks(θs−θr) α wm∫ Θs Θi [Θs+Θ−2Θi]Θ1 2−1 m[(1−Θ1 m)−m+(1−Θ1 m)m−2]dΘ (10) 3. Material and methods The proposed 1D laboratory method consists of a downward 1D infiltration at saturation followed by a sequence of drainage curves obtained, these last ones, by applying increasing tension-gradients between the surface and the base of the soil core. This compound measure allows obtaining a transient infiltration curve followed by successive steady-state flows with increasing slopes (Fig. 1). In order to distinguish between the saturated hydraulic conductivity estimated from the transient cumulative infiltration to that calculated from the steady-state, from now on this term will be defined as K Hk and K s , respectively. 3.1. Estimation of K Hk , S and θ s from the 1D cumulative infiltration curve S and K Hk were calculated by applying the Sequential Infiltration Analysis (SIA) (Moret-Fern´ andez et al., 2021b) to the 1D transient cumulative infiltration curve, which runs from t =0 to t at which free drainage begins (Fig. 1). The SIA method estimates S and K Hk of the upper soil layer by fitting 4 T (Eq. (9)) to increasing time series, and computes the RMSE as a function of the number of samples. The optimal infiltration time, t o , is defined by the minimum RMSE, and the actual S and K Hk are the corresponding values calculated for an infiltration curve of time t o . The Levenberg-Marquardt (More, 1978) optimization algorithm was incorporated in the inverse analysis. A total of 50 increasing times ranged from 50 to 100 s to the whole available infiltration data were considered. The thickness of the soil surface layer was defined as the position of the wetting front advance (WFA) at time t o , which was calculated as (Lassabatere et al., 2009): WFA =I1D(t0) θs−θi (11) where I 1D (t 0 ) was calculated by applying the previously optimized K s , S and β =0.6 to. D. Moret-Fern´ andez and B. Latorre Journal of Hydrology 615 (2022) 128625 4 I1D(t0) = St 1 2 0+(2−β 3Ks)t0+K2 s 9S(β2−β+1)t 3 2 0+2(β−2)(β1)(1−2β) 135 K3 s S2t2 0 (12) More details about the SIA procedure can be found in MoretFern´ andez et al. (2021b). The θ s was calculated as. θs= (If/L)+θr(13) where L is the length of the soil column and I f denotes total infiltration at the time when the wetting front reaches the bottom of the soil column and the steady-state begins (Fig. 1). 3.2. Estimation of α dr , n and K s from the succssive steady-states curves As above mentioned, the drainage process consisted of applying successive tension gradients between the soil surface, h t , and the base, h b , of the soil core. To this end, the following procedure was used: 1. The first pair of tensions, defined by atmospheric conditions at the top and bottom (h t =h b =0 cm) of the soil column allowed calculating K s according to Eq. (3). For this specific case, dz =dH ⇒ K s =q, where q corresponded to the slope of the correspondig steady-state curve (Fig. 1). 2. Next, to allow a partial desaturation of the soil surface, and thus facilitate the entry of air to favor the drainage of the lower soil layers (Sarkar et al., 2019), a second set of tension of h t =h b =-5 cm was applied. From this time, incresing h b values, keeping h t =-5 cm, were selected. This process allows obtaining a seccuence of straight lines of increasing slope, which are related to the incresing h b values (Fig. 1). In this case, the q value calculated for each pair of h t -h b values correspondeded to the slope of the respective drainage sections. Given this analysis corresponds to a drainage process, the estimated α will be correspond to α dr . 3. Once K s and corresponding q were obtained, the optimal n and α dr values were calculated as follows: - Firstly, a numerical procedure to estimate the bottom tension of the soil column, h n , was developed. Given a soil column of height L and assuming K s , q, h t , α dr and n as known values, h n can be iteratively calculated from Eq. (2) as: h2=ht+dL(Ks q−1)→Kh2=Ks{1− ( α h2)n−1[1+ ( α h2)n]−m}2 [1+ ( α h2)n]m 2 z2=0+dL h3=h1+dL(Kh2 q−1)→Kh3z3=z2+dL hn=hn−1+dL(Khn−1 q−1)→Kh3zn=zn−1+dL where z n =L. •Next, taken K s , q and h t as known values, the Q= |hb−hn|objetive function was defined, where h n is the bottom soil tension caclulated for different combinations of α dr and n. Given a α dr value, the optimimum n was computed as the n value that gives a minimum value of Q. This optimization, which was carried out using the LevenbergMarquardt algorithm, was repeated for 1000 different α dr values ranged from 0.001 to 1 cm −1 . This analysis generated different Qisolines, one for each h b , where each line contains the α dr and the corresponding n that results a minimum Q value. The actual α dr and n values corresponded to the crossing-point of the different Q-isolines. •Once obtained the optimum α dr and n values, the soil tension profiles were calulated from the (h t , h 1 , ….h n ) and (0, z 2 ,…. z n ) relationships. The water content profiles were esitmated by applying the van Genuchten function (Eq. (4)) to the corresponding soil tension profile. 3.3. Numerically generated data This new method was first tested on infiltration-steady-state curves generated in synthetic soils. To this end, the theoretical θ s , K s , α and n Fig. 1. Example of a transient infiltration curve followed by successive steady-state drainages obtained after applying a tension gradient between the surface and the base of the synthetic loam soil core of 2.5 cm height. I f denotes total infiltration at the moment the wetting front reaches the bottom of the soil column. D. Moret-Fern´ andez and B. Latorre Journal of Hydrology 615 (2022) 128625 5 values of loamy sand (LS), loam (L) and clay (Cl) soils (Carsel and Parrish, 1988) were compared to the corresponing optimized ones. The infiltration curves were simulated with the HYDRUS-1D software (ˇ Simunek et al., 1999). The van Genuchten (1980) and Mualem (1976) models for water retention and unsaturated hydraulic conductivity functions were used (Table 1), where no hysteresis phenomenon was considered. Homogeneous soil columns of 2.5 cm height discretized in a 1-D mesh of 500 cells were employed. Previous numerical analysis demonstrated that, under this discretization, the solution was grid independent. No evaporation rate was considered. The synthetic infiltration-steady-state experiment was divided in two parts. The first one corresponded to a 1D infiltration with saturation conditions at the soil surface and free drainage at the bottom. The initial soil tension of the soil matrix was -10 7 cm, and the experiment lasted until water started to drain by the bottom and steady-state curve stabilized. This first experiment allowed calculating θ s , S and K Hk (see section 3.1). The second experiment consisted of generating successive drainage curves by applying increasing tension gradients between the top and bottom of the soil. In this case, an initial tension of zero cm was defined in the soil matrix. The first and second pairs of tensions were h t =h b =0 and h t =h b =-5 cm, respectively. These values are basad in tension ranges previously evaluated. Keeping h t =-5 cm, the remaining h b values were −8, −10, −15, −20, −25, −30 for LS and −10, −20, −40, −60, −80 and −120 cm for Cl and L soils, respectively. On the basis of previous synthetic experiments, the time intervals defined for each h b were 200 s for LS and 1000 s for Cl and L. The q values were calculated from the central part of each steady-state section, ensuring that the water flow was stabilized. K s , which was directly calculated by Darcy’s law (Eq. (3)) from the steady-state corresponding to free drainage, was compared to K Hk . α dr and n were estimated from the crossing point of the different Q-isolines obtained from the different combination of h t and h b (see section 3.2). Taken θ i as measurable data (in our cases θ i =θ r ), α for a wetting process, here denoted as α w , was calculated by applying K s and n (obtained from the drainage process) and the corresponding θ s and S (calculated from the 1D infiltration curve analysis) to Eq.(10). The synthetic steady-state experiment was also used to check the time needed by a soil to stabilize the water flux curve when a change in h b occurs. To this end, a h t =-5 cm and an initial soil matrix tension of −10 cm were selected. The experiment was performed on a Cl soil (with the lowest K s ), where h b changed from −60 cm to −80 cm. The time for steady-state stabilization was determined as the time after which the difference between the slopes of the steady-state and infiltration curve were less than 5 %. 3.4. Experimental data 3.4.1. Double disc tension gradients system The device used in this new methodology is based on the same principles as the unit-gradient experiment (Sarkar et al., 2019), where the soil core is placed between two tension disc infiltrometers, one located on the soil surface and the second under the core base. Three different parts can be distinguished: (i) a base disc, (ii) an air vacuum system, and (iii) a surface disc infiltrometer (Fig. 2). The base disc consists of a perforated base (5 cm internal diameter -i. d.- and 0.7 cm high) contained in an aluminum receptacle of 10 cm diameter (Fig. 2). The surface of the perforated base is covered with a dry nylon mesh of 10 μ m pore size, which is hermetically closed against the aluminum receptacle. After testing different nylon meshes (from 10 to 44 μ m pore size), we concluded that the 10 μ m mesh presented the best K s - h maximum relationship, with K s =0.0564 cm s −1 and h maximum = -100 cm. Except for very high permeable sand, this mesh was conductive enough to determine the hydraulic properties of most of the soils. For the particular case of high permeable sand, the unit-gradient method could perfectly be applied. A 2.5 cm-high soil column consisting of a dry soil sample contained on a stainless-steel cylinder (5 cmi.d.) is placed on the nylon mesh. A thin air pathway (≅1 mm thinness) is left between the external wall of the stainless-steel cylinder and the aluminum ring. This path will allow the nylon mesh to be saturated by adding water, if necessary. The bottom of the aluminum receptacle, which has a draining hole of 10 mm diameter, is hermetically attached to a vacuum chamber that will receive the water drained from the wetted soil column. This vacuum chamber is connected to a water column manometer ([0, −100] cm tension) and a vacuum system that generates a suction on the base of the soil core, once the soil has been saturated. Two open/closed valves are placed between the vacuum chamber and vacuum system: a first one (valve 1 in Fig. 2) to connect the two systems and a second one (valve 2) to allow for atmospheric conditions at the base of the soil column, when necessary. Previous experiments showed that, once the nylon mesh was saturated, the tension on the mesh surface was the same as that inside the vacuum chamber. The vacuum system consists of bubbler tower connected to a large and closed water tank located at 2 m on the ground. A long pipe, the end of which is 2 m below the water tank, allows water to flow from inside to the outside of the tank. The water falling through the pipe (valve 3 in Fig. 2), which is regulated by a flow valve, generates a suction inside the tank. The air inlet into the tank, which allows the water to escape, comes from the 110 cm long bubbling tower which is connected to the vacuum chamber. The bubbling tower has a movable air inlet pipe that allows the exact tension to be applied to the vacuum chamber. When valve 3 is open (Fig. 2), the fall of water outside of the tank generates a suction inside the vacuum chamber and bubbling tower. The suction in the vacuum system is set by the movable air inlet pipe of the bubbler tower, and the suction force is regulated by valve 3. The more valve 3 is open, more bubbling intensity there will be inside the bubbling tower. Infiltration water through the soil surface comes from a tension disc infiltrometer of 4.7 cm diameter, which connected to a 15 cm long bubbler tower, allows infiltrations at tensions ranging between 0 and −10 cm. The disc of the infiltrometer is covered by a nylon mesh of 44 μ m of pore size. The water reservoir of the disc infiltrometer is 30 cm high and 2.0 cm i.d. Compared to the 5 cm i.d stain-steel cylinder, the slightly smaller diameter of the upper disc leaves part of the soil surface open, ensuring soil ventilation. A crucial point of the double membrane method is ventilation of the soil, because without air access the soil water content cannot change with variable suction (Sarkar et al. 2019). Although the ventilation of soil in cores could be provided by lateral holes in the cylinders, this would induce some additional difficulties in the handling of the cores and would exclude saturated percolation (Sarkar et al. 2019). A ±0.5 psi differential pressure transducer (PT) (Microswitch, Honeywell), connected to a datalogger (CR1000, Campbell Scientist Inc.), is installed at the bottom of the disc infiltrometer water-supply reservoir (Casey and Derby, 2002). The time evolution of the water level drop inside the infiltrometer reservoir was visible at real time through the screen of a laptop connected to the datalogger. This makes it possible to detect when the steady-state has been reached. 3.4.2. Setup of double disc tension-gradient system Before starting any measurements, the stainless-steel cylinder was Table 1 Theoretical values of initial (θ i ), saturated (θ s ) and residual (θ r ) water content, α and n parameters of the van Genuchten (1980) water retention curve, saturated hydraulic conductivity (Ks), sorptivity (S) of the synthetic soils. θ i θ r θ s α n K s S cm 3 cm −3 cm −1 cm s −1 cm s −0.5 Loamy sand 0.057 0.057 0.41 0.124 2.28 4.05 10 - 3 0.1025 Loam 0.079 0.078 0.43 0.036 1.56 2.88 10 - 4 0.0367 Clay 0.213 0.068 0.38 0.008 1.09 5.55 10 - 5 0.0076 D. Moret-Fern´ andez and B. Latorre Journal of Hydrology 615 (2022) 128625 6 primary filled with sieved and air-dried soil. The soil core, which was covered with a removable plastic lid, was weighted to obtain the initial gravimetric water content, W i , which, for our particular case, corresponded with the residual gravimetric water content, W r . To setup the system, the nylon mesh of the base disc must be initially dry. This allows air to escape from the soil matrix as infiltration progresses. Placing the soil core plus lid in an inverted position (with the lid in contact with the ground), the bottom disc plus the vacuum chamber set, also in an inverted position, was placed and fixed on the soil cylinder. Once the soil was in contact with the nylon mesh, the set was again inverted leaving the lid on the top. This process allows the soil core to be placed on the bottom disc without loss of material. The top lid was then removed and the soil sample was ready to be wetted. At this time, valve 2 was opened (Fig. 2) to allow atmospheric conditions at the base of the soil core. Next, the disc infiltrometer, set at 0 cm tension, was placed on the surface of the soil core and infiltration at saturation conditions started. The infiltration lasted until free drainage through the base of the core began, and a steady-state curve could be visible on the laptop screen (Fig. 1). At this moment, valves 1 and 2 opened and closed, respectively, and a −5 cm tension was applied to the bubbler towers of the disc infiltrometer and the vacuum system, respectively. To generate a suction inside the tank of the vacuum system, and hence on the vacuum chamber, valve 3 opened and the water started to flow down through the 2 m pipe. The bubbling frequency inside the vacuum system was regulated with valve 3. This new pair of tensions lasted during 5–10 min, until a steady-state line was also visible on the laptop screen. In principle, the infiltration process at saturation should be enough to saturate the nylon mesh, needed to generate a vacuum inside the vacuum chamber. However, in some cases, the no-uniform advance of the wetting front may produce an irregular nylon mesh saturation, which prevents creating a complete vacuum inside the vacuum chamber. Once the soil core has been saturated, this problem can be solved by adding some water between the external wall of the cylinder and the aluminum ring. This process, which allows saturating the nylon mesh, activates the vacuum process without altering the infiltration curve. Keeping constant the tension in the disc infiltrometer (h t =-5) cm, increasing h b were next applied on the vacuum chamber by submerging the mobile tube of the vacuum system bubbler tower. At the end of the experiment an 1D cumulative infiltration curve followed by a sequence of drainage steady-state lines were obtained (Fig. 1). Once finished the last stead-state measure, the soil sample was again saturated by opening valve 2 and leaving h t =0 cm. Next, the infiltrometer was removed from the soil surface and the soil sample was poured into a jar. To this end, valves 1 and 2 were closed and the vacuum chamber disconnected for the vacuum system. Following, a jar (with higher diameter than the soil core) was placed upside down on the soil sample, and the pot plus vacuum chamber, aluminum receptacle and cylinder set were turned over so that the soil fell down by gravity into the pot. Because of valves 1 and 2 are closed, the vacuum inside the vacuum chamber prevents the water inside it falls down on the soil core when the system is inverted. This procedure allowed retrieving the saturated soil sample, which was subsequently weighted, dried at 105 ◦C during 24 h, and weighted again to obtain the gravimetric saturated water content, W s . The soil bulk density ( ρ b ) was calculated as the quotient between the core volume and the dry-weight of the soil. The volumetric saturated water content obtained gravimetrically, θ s_w , and θ r were calculated as the product between ρ b and the corresponding W s and W r . 3.4.3. Experimental soils The method was applied of four sieved soils of different texture: sand, loam, clay loam and clay (Table 2). While loam and clay loam soils were sieved at 2 mm, a 0.25 mm sieving was applied to the clay. A 2.5 cm high soil column was employed in loam, clay loam and clay soils. This short soil column was selected in order to increase the homogeneity of the core. However, because the high permeability of the sand, a 5 cm length cylinder was applied to this material. The pair of h t -h b employed in the experimental soils are summarized in Table 2. The soil organic Fig. 2. Scheme of the double disc plus air-vacuum system. D. Moret-Fern´ andez and B. Latorre Journal of Hydrology 615 (2022) 128625 7 carbon content was measured with an improved chromic-acid digestion and spectrophotometric procedure (Heanes, 1984), and the electrical conductivity of saturated paste extract was measured using a standard method (Rhoades, 1982). As described in synthetic soils (section 3.3.), θ s , S and K Hk were calculated from the inverse analysis of the 1D cumulative infiltration curve, K s from the steady-state at free drainage, α dr and n from the crossing point of the different Q-isolines obtained from the different combination of h t and h b (see section 3.2) and α w by applying K s , n, θ s , θ i =θ r and S to Eq. (8). The measured K s was compared to K Hk and the optimized α dr and n were compared with those values calculated with a pressure cell method, PC (Moret-Fern´ andez et al. 2012). The volumetric water content in PC was measured at air-dried soil conditions, at saturation and at pressure heads of –5, –15, –30, –100 and –500 cm. More detalis of the PC experiment can be found in Moret-Fern´ andez et al. (2012). Finally, the estimated α w was compared to the corresponding value calculated by applying the estimated n and α dr to the I H index (Eqs. 6 and 7). 4. Results and discussion 4.1. Synthetic soils The experiment carried out on the synthetic clay soil to check the time required to stabilize the water flow after changing h b from −60 to −80 cm showed that the steady-state was reached after about 150 s (Fig. 3). This short stabilization time could be explained because of the short column used in the experiment, the high tensions applied at the base of the soil core, which accelerate the draining process, and probably also because only a small part of the soil profile takes part in the drainage process. Taking into account that time stabilization decreases with increasing the soil permeability, these results suggest that the presented method is relatively fast to implement: from several minutes in sand to several hours for clay. Fig. 4 shows that, after excluding the extremes within the different steady-states (between black dotted lines), the q values calculated from the drainage curve contained between the solid and dashed blue lines correspond to a zone of straight lines within each steady-state. Although a good relationship was obtained between theoretical (Table 1) and optimized S obtained from the analysis of the 1D cumulative infiltration curve (Fig. 5b), the differences between K Hk and K s increased in soils of finer textures (Fig. 5a). These differences should be attributed to the constant β =0.6 employed in the inverse analysis, whose value does not correspond to those of fine-textured soils. As reported by Lassabatere et al. (2009), the β value increases in soils with finer texture, however, the small influence of β on the 1D cumulative infiltration (Latorre et al., 2018) prevents this parameter from being optimized from 1D experiments. In summary, given β is an unknown parameter, the inverse analysis of a 1D infiltration curve allows only an approximation of K. A significant relationship with p less than 0.1 and slope close to one was obtained between the theoretical θ s (Table 1) and the corresponding value estimated from the cumulative infiltration curve (Fig. 5c). These results indicate that, while K Hk should not be used to optimize α and n from a multiple-tension drainage process, the consistent estimations of S and θ s allows these parameters to be employed to approach α w . Fig. 6 shows the Q-isolines estimated for LS, L and Cl soils. Each of the points of each Q-isoline is the value of n that minimizes the objective function Q within each α value. In this case, the 1000 values of α swept within the [0.001, 1] cm −1 interval used in each Q-isoline covers all α Table 2 Textural properties, organic carbon contents and electrical conductivity, EC, of the different experimental soils and pair of ht-hb values used in the experiments. Sand Silt Clay Organic carbon EC (1:5) Pair of h t -h b values g kg −1 dS m −1 (25 ◦C) cm Sand 1000 – – – 0.11 0–0, 5–5, 5–10, 5–15, 5–20 Loam 280 470 250 11.7 0.99 0–0, 5–10, 5–20, 5–30, 5–40, 5–60, 5–80 Clay Loam 205 497 298 19.9 0.43 0–0, 5–10, 5–20, 5–30, 5–40, 5–60 Clay 151 344 465 12.4 0.52 0–0, 5–10, 5–20, 5–30, 5–40, 5–60, 5–80, 5–100 Fig. 3. Example of time required by a synthetic clay soil to stabilize the water flow after modifying h b from −60 to −80 cm. Red continuous line denotes the projection of the regression line between time 1000 and 1500 s. Fig. 4. Example of steady-state sections (between solid and dashed blue lines) used to calculate the different q for the different h b applied to a synthetic loam soil. Black lines represent the times at which top and bottom tensions were changed. D. Moret-Fern´ andez and B. Latorre Journal of Hydrology 615 (2022) 128625 8 Fig. 5. Relationship between theoretical (a) K s , (b) S and (c) θ s values (Table 1) and the corresponding K Hk , S and θ s calculated from the analysis of the 1D transient cumulative infiltration curve and the four synthetic soils. Fig. 6. Q-isolines at the different soil tensions estimated for the three synthetic soils. Dotted grey lines and red cross denotes the optimized (Opt) and theoretical (Theor) α and n values, respectively. Numbers in legend denote the applied tension at the bottom of the soil core and S is the line for α -n relationship obtained by applying the measured K s , θ s , θ s , n and S to Eq.(10). D. Moret-Fern´ andez and B. Latorre Journal of Hydrology 615 (2022) 128625 9 range for soils. It can observe that in the three analyzed synthetic soils all Q-isolines cross in a single point (dotted grey lines) that corresponds to the theoretical α dr and n values (red cross). This crossover point was calculated as the point with the smallest standard deviation between the optimized n values for each α and the different applied tensions. Overall, a robust relationship between theoretical and optimized α dr and n was found (Fig. 7). An also consistent relationship was obtained between theoretical K s and the corresponding one estimated by Darcy’s law (Fig. 7). The α and n relationship calculated by applying to Eq.(8) the measured θ i and calculated θ s , S and K s , shows a kind of exponential line (dashed red line in Fig. 6) that crosses on the optimum n and α values. Because of the synthetic soils were defined without any hysteresis phenomenon, α w = α dr . Fig. 8 shows an example of the soil tension and water content profiles simulated in a loam soil during the drainage process. In this case, it can observe that water content and soil tension decrease in depth and with increasing h b . However, while important differences in θ and h profiles between two successive h b is observed at low values of h b , these differences decrease with increasing values of h b . This behavior may suggest that, at high h b , the method can lose effectivity in determining h b and, therefore, in estimating n and α . Thus, these results indicate that, depending on the soils, there is a maximum h b threshold (low for sand and higher for clay textures) above which the method begins to loss effectiveness. 4.2. Experimental soils The experimental design developed to measure the soil hydraulic properties from a tension-gradient percolation process resulted to be accurate, simple, inexpensive, fast and easy to implement. Accurate because the tension is directly measured with a water column manometer; simple and inexpensive because the vacuum system does not employ any electronic mechanism and only a datalogger attached to a single pressure transducer is needed; and quick and easy to implement because the complete experiment, including 8–10 tensions, takes no more than two hours and the soil sample is not touched or moved at any time. On the other hand, compared to other methods, the initially dry soil samples employed in the experiment allows obtaining additional information about the wetting processes, and therefore, to potentially characterize the soil hysteresis phenomenon. Furthermore, compared to the Sarkar et al. (2019) method, where the suction is generated by a water column drop, the employed air-vacuum system removes the head losses caused by the water flow through the 2 m length column, and also eliminates the negative effect of the possible entry of air bubbles inside the pipe. Finally, this design is perfectly compatible with the quasi unitgradient percolation method. Fig. 9 shows the experimental infiltration-steady-state curves obtained for the four sieved soils. Note that the greater drainage of water observed in sand should be attributed to the double-length column used in this subtract. The highest n value commonly observed in sand made that only up to 30 cm of tension could be applied (Table 2). On the other hand, the percolation curve with a greater slope at higher tensions observed in clay could be attributed, as will be seen later, to the lower α values of these kind of soils, which makes almost all the soil pores remain filled with water at high tensions. Under these conditions, the slope of the drainage curve becomes proportional to the applied soil tension (Eq. (3)). However, as the soil core base pores empty due to high suction (i.e. loam and clay loam soils with higher α values), K(h) decreases and thus the slope of the drainage curve. The selected experimental soils were no-saline and the corresponding bulk density varied from 1.18 to 1.51 g cm −3 for clay and sand, respectively (Table 3). A significant relationship (Fig. 10) was obtained between the θ s and θ s_w , where on average θ s_w was 2 % smaller than the corresponding obtained from the cumulative infiltration curve. Although both parameters gave comparable results, given the ease of estimating θ s , which does not require any soil handling, from now on, α w will be calculate using θ s . Except for clay and loam soil columns, the WFA values estimated for the sand and clay loam soils were close to the actual cylinder lengths (Table 3). These results suggest that the analyzed sand and clay loam soil columns were, overall, homogeneous. For the particular case of the loam soil column, the short WFA value should be attributed to a heterogeneity in the soil column, probably due to a soil compaction during the water infiltration experiment due to the disc infiltrometer weight. This made that K s was one order of magnitude lower than K Hk (Table 3). In contrast, the poorer estimates of WFA observed in clay could be attributed to the fact that Haverkamp et al. (1994a, 1994b) model does not work well in clay soils (Lassabatere et al., 2009). This caused Eq.(7) did not fit well with the experimental curve, forcing the SIA procedure to stop before the wetting front reached the soil base. Although a significant relationship between K s and K Hk was found (Fig. 11a), on average, K Hk (0.011 cm s −1 ) was one order of magnitude higher than K s (0.0029 cm s −1 ). Again, these large differences should be attributed to the high K Hk obtained in loam and clay (Table 3) soils. Although K s for sand was three times lower than K Hk , both values were within the same order of magnitude. These differences could be attributed to the high permeability of the sand, which, together with the short column, probably prevented the infiltration time from being long enough for a more accurate estimate of K Hk . The S values estimated for the different soils were within the same order de magnitude of the corresponding values measured for the same soils using an upward infiltration method (Moret-Fern´ andez and Latorre, 2017). In summary, these results showed that, although Haverkamp et al. (1994a, 1994b) model allows good approaches of K, the dependence of K Hk on the β parameter and on the soil homogeneity suggests that this parameter should not be employed in the inverse analysis of the drainage curves. In contrast, the robustness of the SIA procedure to estimate WFA and S indicates that the measured 1D infiltration curve provides interesting information to check the homogeneity of the soil column and valuable data about the soil wetting process, respectively. Given that K s is directly estimated by Darcy’s law, only α dr and n need to be optimized. Similarly to what was observed in synthetic soils, the Qisolines analysis showed in all the experimental soils a single crossingpoint, whose optimum α dr -n values (dotted grey line in Fig. 12) varied with the type of soils: the lowest and highest α dr -n values corresponded Fig. 7. Relationship between theoretical α , n and K s , and the corresponding values calculated for the four different synthetic soils from the steady-state lines obtained after applying increasing tension gradients between the surface and the base of the synthetic core. D. Moret-Fern´ andez and B. Latorre