A general mathematical framework for modelling soil-water retention behaviour

A new constitutive model for the soil-water retention behaviour of unsaturated soils is proposed, able to reproduce the main drying and wetting paths, the cyclic retention behaviour and its dependence on the specific volume. The most significant aspect is the inclusion of the evolution, with the specific volume, of the degree of saturation when suction tends to zero in wetting paths considering the presence of entrapped air bubbles. The model is used to reproduce with success the drying/wetting cycles of two Pearl


Introduction
The ability to reproduce soil-water retention behaviour, usually through the definition of the relationship between the degree of saturation, S r , and suction, s, is one of the key aspects when modelling the behaviour of unsaturated soils. This relationship is the known soilwater retention curve, which should not be called as characteristic curve because it depends on soil's void ratio or specific volume. Several laboratory tests, performed in the last two decades, shown that the soilwater retention behaviour is quite complex, presenting hysteretic behaviour between drying and wetting paths ( [1], [2], [3], [4]), and a dependence on the void ratio ( [5], [6], [7]) and the temperature ( [5], [8]).
The hysteretic response of soil-water retention behaviour is a consequence of several factors. Lu & Likos [9] pointed out three possible reasons for the hysteretic response: x The difference in the contact angle between wetting and drying (liquid angle between the solid particles and the surface of tension); x Geometrical effects due to the geometric configuration of the voids; x The entrapment of air bubbles during the wetting stages.
Due to the coupling, established through the void ratio, between the soil-water retention behaviour and the mechanical behaviour of the solid skeleton, in a deformable soil, the aspect of the retention curve is always dependent on the value of the void ratio. Romero [5] has shown a translation of the main drying curves in the direction of higher suctions with the decrease of the void ratio, and Tarantino [6] and Vanapalli et al. [7] presented experimental results showing a dependency of the retention curves on the initial void ratio. These results also indicate that this dependency decreases for high suction values.
The most important constitutive models for the soilwater retention behaviour can be distributed by four different groups: single mathematical function, models incorporating hysteretic behaviour, models incorporating the dependence on the specific volume and models with both the dependence on the specific volume and hysteretic behaviour. Although the existing possibilities, in this paper, only the models in isothermal conditions formulated with unimodal functions are analysed.
The initial formulations ( [10], [11], [12]) belong to the family of reversible non-linear models assuming that they are established by fitting mathematical functions to the experimental data. All the functions assume incompressible porous medium and tend to describe the experimental results without taking into account the mechanical and hydraulic history. No interpretation is given, except for the air-entry suction value and the residual saturation, which are derived from geometrical constructions. Sillers & Fredlund [13] and Sillers et al. [14] provide a detailed discussion on the mathematical attributes of the different curves proposed, and Leong & Rahardjo [15] made an exhaustive review of the most important equations. The previous equations, although useful in practical applications, only give a basic definition of soil-water retention behaviour, ignoring the existence of the capillary hysteresis and the dependence on the specific volume of the retention behaviour.
Several mathematical models have been proposed to reproduce the hysteresis of the soil-water retention behaviour ( [16], [17], [18], [19], [20]). Pham et al. [21] made an extensive review of the existing hysteretic constitutive models. Generally, the possible states of saturation are bounded by the main drying and wetting curves, which are respectively the curve of highest possible degrees of saturation upon drying and the curve of lowest possible degrees of saturation upon wetting. These are the curves corresponding to the first drying starting from the fully saturated state and to the first wetting starting from the driest possible state (different denominations can be found for these two branches of the hysteresis domain). Inside the domain limited by the main drying and wetting curves, the behaviour can be considered reversible (if the main curves are not reached) or not, being called scanning lines when reversible.
Several mathematical models have been proposed to reproduce the dependency on the specific volume of the soil-water retention behaviour ( [22], [23]). While the dependency on the void ratio is of particular interest here, the failure to reproduce the capillary hysteresis is inherited from the non-linear fully reversible models. Recently, several models have been proposed taking into account both the hysteretic behaviour and the dependency on the specific volume ( [24], [6], [25], [26]).
Finally, Sakaki et al. [27] extended the concept of soil-water retention model into negative suction values to deal with the entrapped air. The possibility of having unsaturated states with ‫ݏ‬ = 0 at the end of a wetting stage was also observed in micromechanical models of spheres' arrangements by Pereira [28]. Very few of the existing models consider the possibility of having entrapped air when ‫ݏ‬ → 0. This is one aspect that is taken into account in the model proposed here.
Therefore, the main objective of this paper is to present a soil-water retention model incorporating the hysteretic behaviour, the dependency on the void ratio and the possibility of having an unsaturated state when suction tends to zero in a wetting stage. The proposed model is restricted to isothermal conditions because more studies are needed to characterise better the evolution of the main wetting and intermediate stages with temperature.

General mathematical framework
In this section, a general mathematical framework able to reproduce the soil-water retention behaviour is presented. In this framework, different mathematical functions can be adopted for the main drying and wetting paths to reproduce the evolution of the degree of saturation, ܵ ୰ , with suction, ‫,ݏ‬ and specific volume, ‫.ݒ‬ The hysteresis, which is the different response between drying and wetting paths at intermediate states, is taken into account. Figure 1 schematically represents the key elements of the soil-water retention framework at constant ‫.ݒ‬ The red and blue lines represent the main drying, ܵ r d , and main wetting, ܵ r w , curves, which are function of ‫ݏ‬ and ‫.ݒ‬ These curves define the limits of the admissible values of ‫,ݏ(‬ ܵ r ) for each value of ‫.ݒ‬ The green line represents a drying path starting from the intermediate point A that converges to ܵ r d with the increase of ‫.ݏ‬ The brown line represents a wetting path starting from the intermediate point B that converges to ܵ r w with the decrease of ‫.ݏ‬ The main drying path is a sequence of drying states that starts from a fully saturated state, and the main wetting path is sequence of wetting states that starts from a completely (as possible) dry state. The main drying and wetting paths are not characteristic curves in the ‫,ݏ(‬ ܵ r ) representation space because they depend on the value of ‫ݒ‬ ( [5], [29]). Additionally, in a highly compressible soil, drying and wetting stages originate significant variations of ‫.ݒ‬ The methodology presented by Sugii et al. [30] is used here to define the main drying, ܵ r d , and wetting, ܵ r w , curves (red and blue curves of Figure 1). Both are defined as isolines of the same ‫ݒ‬ in the ‫,ݏ(‬ ܵ r ) representation space. The admissible ‫,ݏ(‬ ܵ r ) paths (green and brown curves of Figure 1) are limited by the two functions ܵ r d and ܵ r w already defined. In a compressible soil, the ‫,ݏ(‬ ܵ r ) path crosses several main curves obtained at constant ‫.ݒ‬ When an inversion of the suction direction occurs in a process controlled by ‫,ݏ‬ the values of the pair ‫,ݏ(‬ ܵ r ) start to converge to the main drying and wetting curves if suction increases (d‫ݏ‬ > 0) or decreases (d‫ݏ‬ < 0), respectively. When an inversion in the direction of the degree of saturation occurs in a process controlled by ܵ r , the values of the pair ‫,ݏ(‬ ܵ r ) start to converge to the main drying and wetting curves, respectively if dܵ r < 0 or dܵ r > 0. The convergence rate varies with the distance of the actual state to the corresponding main curve.
The degree of saturation increment, dܵ r , is here defined as a function of the suction increment, d‫,ݏ‬ and the specific volume increment, d‫,ݒ‬ given by  where ܵ r d and ܵ r w are the values at the main drying and wetting curves at current ‫ݏ‬ and ‫,ݒ‬ respectively, ߦ r the relative position of ܵ r between S r d and S r w , and ߞ d and ߞ w constants of the model that define how fast the current state ‫,ݏ(‬ ܵ r ) converges to the corresponding main curve.
In this formulation, the convergence to the main curves is exclusively due to ds, being assumed that when ds = 0, dS r is the value that maintains ξ r constant. Therefore, the evolution of S r due to dv is directly dependent on the evolution of S r d and S r w with dv. The main drying curve is defined between ܵ r = 1 (saturated state), when ‫ݏ‬ = 0, and ܵ r = ܵ res , when ‫ݏ‬ increases to its maximum value. The variable ܵ res is defined here as the residual degree of saturation, which is the degree of saturation measured for the greatest ‫ݏ‬ value imposed (it could be a relatively low suction value), considering ܵ r = 0 for the soil dried in an oven, at say 105°C. In purely granular soils, a zero value of ܵ res is expected because it is possible to extract all the water from the soil, achieving a completely dry soil. In clayey soils, ܵ res > 0 can be conjectured because the electrochemical forces do not allow the evaporation of the water molecules adsorbed on the surface of the particles at ambient temperatures, even for very low relative humidity values. The main wetting curve is defined between ܵ sat , when ‫ݏ‬ = 0, and ܵ r = ܵ res . The term ܵ sat takes into account the possibility of having an unsaturated state with zero suction.

Soil-water retention model
In this section, a soil-water retention mathematical model (WRM) is presented taking into account the framework presented in Section 2.

Main drying and wetting curves
The equation presented by van Genuchten [11] (the equation most used in the literature of the soil-water retention models) was assumed as starting point to represent both main curves. The evolution of the main drying curve, ܵ r d , with ‫ݏ‬ and ‫,ݒ‬ is defined as where ݊ d and ݉ d are constants and ܽ d ‫)ݒ(‬ a function that makes the main drying curve dependent on ‫.ݒ‬ The evolution of the main wetting curve, ܵ r w , with ‫ݏ‬ and ‫,ݒ‬ is given by where ܵ sat represents the degree of saturation at the end of a main wetting path ‫ݏ(‬ = 0), ݊ w and ݉ w are constants and ܽ w ‫)ݒ(‬ a function that makes the main wetting curve, together with ܵ sat , dependent on ‫.ݒ‬

Evolution of function ‫)ݒ(ܽ‬
Functions ܽ d and ܽ w control the position of ܵ r d and ܵ r w , respectively, in the ‫ݏ‬ axis. With the decrease of ‫,ݒ‬ the values of the degree of saturation of the main curves increase ( [5], [28], [29]). The following equation is proposed to describe the evolution of ܽ d and ܽ w with ‫ݒ‬ where ݅ = {d,w}, ܽ ଵ the maximum value of function ܽ (when ‫ݒ‬ → 1), ܽ ஶ the minimum value of ܽ (when ‫ݒ‬ → ∞) and ߙ , ߚ and ߛ constants. The present formulation moves the main drying and wetting curves into the higher suction values with the decrease of ‫.ݒ‬

Evolution of function ܵ sat ‫)ݒ(‬
The existence of ܵ sat < 1 at zero suction and its dependence on ‫ݒ‬ is a consequence of the occurrence of isolated air bubbles, which cannot be expelled from the soil on wetting, or when the air phase is connected. Its value could converge to 1 with the increase of ‫ݒ‬ ( [31]), assuming that permeability increases with increasing specific volume, favouring bubbles expulsion, or decrease with the increase of ‫,ݒ‬ when the air phase is assumed to be interconnected. Considering only the first case, the evolution of ܵ sat with ‫ݒ‬ is here defined by where ܵ ଵ is the minimum value of ܵ sat (when ‫ݒ‬ → 1), and ߚ S and ߛ S constants. Figure 2 presents a tridimensional representation of the main drying and wetting surfaces in the (s, v, S ୰ ) representation space, with the values indicated there. The main drying surface is represented with some degree of transparency. As shown in the plots on the right side of this figure, for fixed values of ‫,ݒ‬ the main drying curve is always higher than the main wetting curve, and the decrease of ‫ݒ‬ increases the value of ‫ݏ‬ for the same value of ܵ r .

Model constants
In its most general form, this model has 21 constants: Main wetting curve: Hysteretic behaviour: ߞ d , ߞ w .
Several simplifications can be done to reduce the number of constants. For example: x When dealing with soil particles without clayey minerals, it can be assumed ܵ res = 0. x When the evolution of ܵ sat with ‫ݒ‬ is not relevant, or there is no data available to calibrate this part of the model, the function ܵ sat becomes a constant, which can be equal to 1. In this case, ܵ ଵ = 1 and ߙ S , ߚ S and ߛ S are not relevant, reducing by four the number of constants.
x When a fixed ratio between ܽ w and ܽ d is assumed for simplification, only the evolution of ܽ d with ‫ݒ‬ needs to be characterised. In this case ߦ a = ܽ w ܽ d ⁄ is constant and ߙ w , ߚ w , ߛ w , ܽ ଵ w and ܽ ஶ w are not relevant, reducing by four the number of constants.
x When only data from drying paths that started in saturated stages are available, the constants that control the behaviour between main curves and the main wetting curve cannot be obtained. x Finally, when van Genuchten's equation is adopted considering only the main drying equation without dependency on ‫,ݒ‬ then ܽ d = ܽ ଵ d = ܽ ஶ d .
In Pereira [28], the WRM was computed using the Euler method and the fourth-order Runge-Kutta method, together with all partial derivatives necessary to implement the model.

Simulation of laboratory results
The WRM was applied to reproduce the soil-water retention behaviour of two Pearl clay samples prepared with different initial void ratios, by directly imposing d‫ݏ‬ and d‫ݒ‬ ( [28]). Pearl clay (50% silt and 50% clay, with the grain diameter less than 5μm) has liquid limit of 49%, plasticity index of 22% and specific gravity of 2.71. The mineral composition of the clay, determined using X-ray diffraction, identifies kaolinite and few expansive clay minerals. Two unsaturated laboratory tests were selected from the work presented by Sun et al. [31], characterised by the existence of wetting/drying cycles at constant mean net stress, p net , and significant different initial values of ‫.ݒ‬ The samples were prepared by static compaction in a steel mold using dry material and the intended water content or initial suction. Suction was applied using axis translation technique. The main characteristics of the samples tested are in Table 1 and the laboratory results to be simulated are in Figure 3.
The constants of the model were obtained with the Differential Evolution algorithm described in [28], considering ܵ ୰ୣୱ = 0. Two simulations were done: the first using the model proposed, WRM, in which ܵ sat < 1 at zero suction (assuming the presence of isolated air bubbles); and the second adopting a simplified version, SWRM, where ܵ sat = 1. Tables 2  and 3 represent the values of parameters of the WRM and SWRM, respectively.
In Figure 4, the numerical results obtained with the constants presented in Tables 2 and 3 are shown together with the laboratory data. The numerical results are qualitatively acceptable for both WRM and SWRM models. It should be highlighted the capability of the model to reproduce the wetting and drying cycles of Tests 1 and 2. However, the WRM is able to reproduce much better the complex paths that characterise the laboratory results. The incorporation of the possibility of having ܵ sat < 1 at zero suction (isolated air bubbles) and its evolution with ‫,ݒ‬ improves the fitness of the numerical results.

Conclusions
This paper presents a mathematical model able to reproduce the soil-water retention behaviour. The model was defined in tridimensional space ‫,ݏ(‬ ܵ r , ‫)ݒ‬ and incorporates the possibility of having unsaturated states with zero suction. The numerical simulations were performed imposing d‫ݏ‬ and d‫,ݒ‬ and a Differential Evolution algorithm was used to determine the constants of the model for each group of tests.
The model, WRM, and its simplified version, SWRM, were able to successfully reproduce the soilwater retention curves of the Pearl clay samples prepared with different initial void ratios and subjected to drying/wetting cycles. In particular, it was possible to reproduce hysteresis and the transition between successive cycles observed. The simplified version, SWRM, was not as good as the full version, WRM, as the quality improves when assuming the existence of a degree of saturation lower than 100% when suction tends to zero in a wetting path.   Table 3. Values of the parameters of the SWRM, obtained by a Differential Evolution algorithm [28], to simulate the soil-water retention behaviour Tests 1 and 2 [31].