A simple hysteretic constitutive model for unsaturated flow

In this paper we present a constitutive model to describe unsaturated flow that considers the hysteresis phenomena. This constitutive model provides simple mathematical expressions for both saturation and hydraulic conductivity curves, and a relationship between permeability and porosity. The model is based on the assumption that the porous media can be represented by a bundle of capillary tubes with throats or"ink-bottles"and a fractal pore size distribution. Under these hypotheses, hysteretic curves are obtained for saturation and relative hydraulic conductivity in terms of pressure head. However, a non-hysteretic relationship is obtained when relative hydraulic conductivity is expressed as a function of saturation. The proposed relationship between permeability and porosity is similar to the well-known Kozeny-Carman equation but depends on the fractal dimension. The performance of the constitutive model is tested against different sets of experimental data and previous models. In all of the cases the proposed expressions fit fairly well the experimental data and predicts values of permeability and hydraulic conductivity better than others models.

In this section, we derive closed-form analytical expressions for saturation and hydraulic 71 conductivity curves. First, we present the pore geometry of the proposed model and we derive 72 some hydraulic properties which are valid for a single pore. Then, assuming a cylindrical 73 REV of porous media with a fractal pore size distribution, we obtain expressions for porosity, 74 saturated hydraulic conductivity, and saturation and relative hydraulic conductivity curves. The porous media is represented by a bundle of constrictive capillary tubes. Each pore 77 is conceptualized as a cylindrical tube of radius r and length L with periodically throats 78 represented by a segment of the tube with a smaller radius, as illustrated in Fig. 1. 79 Assuming that the pore geometry has a wavelength λ and that the length of the tube L 80 contains an integer number N of wavelengths, the pore radius along the tube can be described  83 where a is the radial factor (0 ≤ a ≤ 1), c is the length factor of the pore throat (0 ≤ c ≤ 1) 84 and n = 0, 1,...,N − 1. The parameter a represents the ratio in which the radius is reduced, 85 and the parameter c represents the fraction of λ with a narrow neck. Note that if c = 1o r 86 c = 0 we obtain a straight tube with radii ar or r , respectively.

87
Based on the above assumptions, the volume of a single tube can be calculated by inte-88 grating the cross-sectional area over the length L as follows: 92 f v is a factor that varies between 0 and 1, and quantifies the reduction in pore volume due to  where ρ is the water density, g gravity, µ water viscosity and h the head drop across the 100 tube.

101
Substituting Eq. (1)inEq.(4) yields: 103 105 f k is a factor that quantifies the volumetric flow rate reduction due to pore throats and varies 106 between 0 and 1. Figure 2b shows the variation of f k as a function of the radial factor a and 107 the length factor c. As it can be expected, low values of a drastically reduce the volumetric 108 flow rate of the tube.

109
If we compare Fig. 2a

116
To derive the expressions for saturation and hydraulic conductivity, we consider as a REV a 117 straight circular cylinder of radius R and length L. The choice of the REV geometry is based 118 on the shape of soil samples commonly used in laboratory tests. Other geometries, such as 119 rectangular REV, can also be considered by introducing minor changes in model derivation.

120
The pore structure of the REV is represented by a bundle of constrictive tubes (as described 121 in the previous section) with a fractal pore size distribution. We also assume that the pore 122 radius r varies from a minimum value r min to a maximum value r max .

134
The porosity φ of the REV can be computed from its definition: Replacing Eqs.
(2)and(8) into Eq. (9), the porosity of the REV can be expressed as: is the porosity of the REV considering straight tubes (i.e., a = 1).

141
The volumetric flow rate Q at REV scale can be obtained by integrating all the pores 142 volumetric flow rates given by Eq. (5) over the entire range of pore sizes: On the other hand, on the basis of Darcy's law (1856), the volumetric flow rate through the 145 REV can be expressed as:

147
where K s is the saturated hydraulic conductivity. Combining Eqs. (12)and (13), an expression 148 for K s is obtained: allowing us to represent media with high porosity, low permeability and low specific surface 157 area. Our model is also able to describe media which have the same porosity but different 158 permeabilities. For example, clay and sand soils have typically similar porosities, but their 159 hydraulic conductivities differ by several orders of magnitude (e.g., Carsel and Parrish 1988).

160
For most porous media, r min /r max ≃ 10 −2 (Yu and Li 2001); then, we can assume that 161 r min ≪ r max . Under the above assumption, the terms r 2−D min and r 4−D min in Eqs. (11)and (15) 162 can be considered negligible. Then, combining the resulting expressions, we can obtain the 163 following relationship between K s and φ:

179
where σ is the surface tension of the water and β the contact angle.

180
To obtain the main drying saturation curve, we consider that the REV is initially fully 181 saturated and is drained by a pressure head h. We assume that a tube becomes fully desaturated 182 if the radius of the pore throat ar is greater than the radius r h gi venbyEq.(18). Then it is 183 reasonable to also assume that pores with radii r between r min and r h /a will remain fully 184 saturated. Therefore, according to Eqs.
(2)a n d( 8), the drying saturation curve S d e can be 185 computed by: 208 where K r is the relative hydraulic conductivity which is a dimensionless function of h and 209 varies between 0 and 1.

210
Combining Eqs. (23)a n d( 24), and using Eqs. (5), (8)a n d ( 14), we obtain the relative 211 hydraulic conductivity for the drying process:  Similarly, the main wetting relative hydraulic conductivity curve K w r (h) can be derived 216 by integrating Eq. (23) over the range of saturated pores (r min ≤ r ≤ r h ): Note that saturation and relative hydraulic conductivity expressions for both drying and  In order to test the proposed relationship between permeability and porosity for different

260
where C KC is a parameter that depends on the specific internal surface area, an empirical 261 geometrical parameter and the tortuosity.

262
For each type of soil, Eqs. (30)a n d( 32) are fitted to measured data by minimizing the 263 root-mean-square deviation (RMSD):

265
where k calc and k dat correspond to the calculated and measured permeabilities, respectively. A 266 logarithmic scale was considered because of the wide range of variation for the permeability 267 values. The fitted parameters for Eqs. (30)a n d( 32) are listed in Table 1 as well as their 268 respective RMSD. It can be noted that, for all soils, the RMSD of the proposed model is 269 smaller than the ones from the KC equation. Figure 3 shows that the proposed relationship   Table 2. It can be 279 noticed that the proposed model shows a significant improvement over the one of Assouline 280 for the Gilat sandy loam (see Fig. 4b). Table 2 Table 3). Then, the fractal dimension D and the radial factor a have been 291 estimated by minimizing the RMSD between calculated and measured values of both drying 292 and wetting curves using an exhaustive search method. Table 3 shows the model parameters    Table 3. The prediction of the wetting curve from the 304 drying curve could be an additional advantage of the proposed model that needs to be verified 305 with a more exhaustive analysis and additional experimental data.