A theoretical retention model for unsaturated uniform and graded soils

Modelling the relation between the degree of saturation and the suction (ie retention curve) is an important challenge for geotechnical engineering. It has a huge influence on the behavior of large soil constructions as levees, embankments, road earthworks. We present here a theoretical model of retention curve which considers physical relations of unsaturated soils. With this approach, there is no need to assume particular shapes of the retention curves, which are a consequence of the physical assumptions. The present study is focused on a theoretical model based on elastic spherical particle arrangement. As a first step a uniform model is presented with a single diameter of soil particle. A second step extends the use of the model to graded soils. The model uses only 5 physical parameters. It is compared with the experimental retention curve of two different samples of glass uniform particles and two different graded soils, a graded glass sample and the Livet-Gavet loam. It shows its ability to model the experimental curves and a better agreement than the former theory of Brooks and Corey (1966). This current publication is funded by the French National Project « Terredurable» (ANR 2011),.


Introduction
Development of numerical calculations and increasing power of microcomputers allow determining the behavior of large soil construction as levees, embankments, road earthworks. All these constructions are compacted at the Proctor optimum with a degree of saturation between 80 to 98% and an unsaturated soil. For a precise estimation of the final state of density of the soil, it is needed to model the compaction of the soil along a wetting path. Furthermore when the earth construction is achieve, it supports drying and wetting events linked to the rainy and sunny different days. These events need to develop a model able to simulate drying and wetting paths There are several ways for modelling retention curve. The first method is to use experimental relations. This approach had many successes with the model (1) of Brooks and Corey (1966) which needs the determination of air entry suction s air and the experimental parameter . Later Van Genuchen (1980) used a new experimental relation (2) able to simulate larger retention curve with 3 experimental parameters (, n, m). Gallipoly et al. (2003) follow the same research scheme and proposed a complete experimental relation (3) with 4 experimental parameters (n). The increase of precision is paid by an increase number of parameters.
The second approach is to simulate the retention curves by an incremental elasto-plastic model (Arairo., 2013). This approach needs the knowledge of the experimental curve of saturation and imbibition to predict the behavior of other retention curve. It uses the parameters of Van Genuchen (1980) and two limit curves The third approach is to consider physical modelling. With this approach, there is no need to assume particular shapes of the retention curves, which is a consequence of the physical assumptions. It is the approach used in this paper. The present study is focused on a theoretical model based on elastic spherical particle arrangement. This approach follows a new research way so that the physical phenomenon will be explained and the number of experimental parameters will be limited. (1) 2. The expression of the theoretical retention curve for the unsaturated soil-study

Hypothesis
Many authors have mentioned the existence of four areas of saturation, each with a distinct behaviour. This assumption is resumed in the design of our model, mostly based on the work of Boutonnier (2007).
-Domain D1 : s ≥ s air and S r ≤ S rair The gaseous phase is continuous in the soil. This state gives a suction s higher than the air entry suction and a lower degree of saturation than the degree of saturation of air entry.
-Domain D2 : s ≤ s air and S rair ≤ S r ≤ S re and u w ≤0 In this domain, free air disappeared. The air is occluded in the soil. The air is in contact with the soil particles. The suction has the effect of increasing the strength of interparticle contacts. The water pore pressure is negative -Domain D3 : S re < S r < 1 and u w >0 Air is occluded into the soil sample and is considered as to be independent air bubbles from the skeleton. Capillary strength has no effect on the contact forces between soil particles. We consider here D3 is corresponding with the case of positive pore pressures with a degree of saturation lower than 1.
-Domain D4 : S r = 1 There is no air in gaseous state in the soil. The soil is saturated. The boundary between D3 and D4 can also be expressed through the pore pressure with a degree of saturation equal to 1.

Deformation of the spherical particles and volume of meniscus: Domain D1
To complete the original model of non-saturation (Monnet and Boutonnier, 2012) written for domain D2 to D4, it was decided to model the D1 domain. Under the action of the water at the interparticle contact, there is suction on the cross section of the meniscus (Fig.1) and a surface tension as assumed previously (Taibi, 1994). The model of soil is composed of uniform spherical particles of radius R.
According to Laplace's law, suction is the product of the surface tension T c by the sum of two principal curvatures of radius r of the meniscus and the radius b of the torus corresponding to the wetted area (4). Under the action of a normal force F N on the contact plane (5), the two balls deforms elastically. The distance between their centers decreases of the quantity  and the contact is no more punctual, but is a circular plane area of radius (a Fig.1). We assume that the balls remain spherical except for the flat portion of the contact. The values of a and  are given by equations (6) to (8) of the elastic Hertz contact between the two spheres. The value of the normal force F N to the contact plane is given by equation (9) which takes into account le value of the suction and the value of the surface tension T c (4) The contact between the particles and the water has a wetting angle of  c which varies in case of drainage or humidification (Gras, 2011). The calculation of the volume of the meniscus is made with a rotational assumption around the horizontal axis through the centre of the two particles ( Fig. 4) and a relative horizontal symmetry with the EC plane contacting the particles. The volume of the meniscus is the difference between the volume of CDAE frustum, the volume of the soil flattened sphere CDAB, and the volume of air of the torus of AFE area (10). The main variables of the calculation are the angles  1 (11) to  5 (12), l the radius of the contact circle (13), b radius of the meniscus torus (14), r' the radius of the circle which connect the two contact circles (15), R1 the radius of the frustum (18).
The main variables of the problem can be defined, as the volume of the frustum (16), the volume of the air (17) and the volume of the particles (18). (16) The meniscus volume was simulated by a Solidworks® model which allows measuring the numerical volume of the meniscus, and comparing with the analytical value of the meniscus (9). The comparison shows that the analytical relations allow finding the meniscus volume of the Solidworks® model.

Different arrangements of soil particles
The different arrangements of uniform soil particles give a variation of void ratio between 1.315 to 0.343 (Taibi, 1994). The SoildWorks® program was used to a precise modelling of the four possible arrangements, tetrahedral, cubic, octahedral, dodecahedral ( Fig.3 to 6). It shows that the void ratio is independent of the radius of the particles, but depends only of the arrangement between particles. The radius of particles has only an influence on the size of the REV (Representative Volume Element). The variable which describes the different phenomenon along the drying path is the ratio between the air bubble and the particles .
The Solidworks model allows determining: -The number and the angle of contact for each particle -The number of total menisci per arrangement For the drying path at nucleation, the model finds: -The radius of the bubble which appears at nucleation -The number of bubbles able to appear inside VER -The degree of saturation at nucleation -The suction at nucleation For the drying path at percolation, the model finds: -The radius of the bubble which is able to percolate -The suction at percolation The theoretical calculation of the relation between suction and degree of saturation can be performed when the menisci are independent (Fig.7). When the meniscuses merge, it is assumed that the suction remains unchanged at wetting. At the drying path nucleation appears first and percolation of air follows.  volume a single meniscus (20). The volume of the meniscus depends of the radius of the meniscus, the radius of the grain and the suction. This allows calculating the suction associated with a void ratio and a water content (Fig.7).

Compactness and homogeneity of granular mixture for a graded soil
The theoretical model is organized from the larger to the smaller of n different Di diameters for the soil squeleton particles, such as (25). The symbol C is the compactness which is equal to the solid volume of the soil sample (26). The relative compactness to the total sample of the class i is noted Ci (27). As a consequence, the total compactness C of the soil sample is given by (28). In the theory, Ri corresponds to the volume refusal on the sieve of diameter Di , based on the total volume of the sample grains (29). Assuming a single density for all aggregates, so this term is the refusal by sieve reduced to the total mass of the sample

Mixture scattered -general case for a graded soil
With a mixture of 3 parts a loosening effect (Fig. 8) and wall effect may appear. If the class 2 prevails the grains are subjected to a loosening effect exerted by the smaller grains from class 3 and a wall effect exerted by class 1. In This can be generalized for n different classes (30) by finding the virtual C compactness of the mixture, considering that the class i is dominant in equation (31). The relative density corresponding to i coefficient is given by (32). The experimental i coefficient is found by interpolation of the experimental curve ( Fig.9): This is the ratio between the density of the mixture to the particle density versus the particle diameter Di. At a first step of calculation, the total compactness C m of the mixture and the porosity n m are known. Now consider the total size curve (Fig.10). The dominant diameter separates the curve into two portions: -At the right of the dominant diameter, large elements floating into small elements -At the left of dominant diameter, small items fall into the gaps left by the dominant grains Each portion has a density of its own and depends only on the prevailing diameter. Then we can consider two separate samples for which two new dominant diameters appear with their own characteristics of compactness, porosity and void ratio. We can determine from part to part all theoretical porosity of all diameters of the particle size curve.
Let us consider the total mixture of n different particles. For a single rank of aggregate i, the volume of the void V vi into the mixture is given by the equation (37). Under these conditions, the void ratio e théo of the theoretical mixture is known. This void ratio is different from the measured void ratio e m of the mixture, since the former refers to a method of standardized compaction where the grains grow to different relative densities. We assume the proportionality between the void ratio of the different classes. According to equation (39) e mi is the theoretical void ratio of class i; e essai .is the known value of mixture measured during the test; e théo is the void ratio of the mixture compacted at the standard energy (Fig.8), and e i is the void ratio of each diameter into the compacted sample. This void ratio e mi allows knowing the arrangement of the particle class. (33) (35) Figure 10. The grade curve of Livet-Gavet loam separated into two portions, large elements on the right and left small particles We assume that small pores are saturated (S ri = 1), and that large pores are dry (S ri = 0). Under these conditions, the water content can be decomposed as a function of particle size (40), with dry pores (S r = 0) for j <i (large diameter) and saturated pores (S r = 1) for j> = i (small diameters). The rank of the pores involved in the air inlet is i,. In relation (40), the only unknown is the rank i. We can define the diameter of the grains involved in the air inlet, and builds a theoretical retention curve for mixtures. The knowledge of e mi and the water content w NSi allows finding the suction of class i (Fig.7). (36)

Comparison with an experimental uniform sample of glass
The theoretical retention curves shows a large variation of degree of saturation for single value of suction. For the wetting path, there is coalescence of menisci at a precise value of the suction and an unstable situation between the state of separated menisci and complete saturation of void between particles. For the drying path the air bubbles are allowed to percolate between particles at a precise value of suction. when the air percolates, it reaches larger void between particles which allowed sudden decrease of the water content and sudden decrease of the degree of saturation. The vertical parts of the retention curves are linked to the fact that the grade curve is known only by measurements points and not by a continuous curve. The simulation of the experience on two glass particle samples (Indarto, 1991) is made along drainage and wetting paths. The samples had respectively a uniform diameter of 80m and 300m, a unit weight of 18,7kN/m² (e = 0.550; dodecahedral). The unit weight of the glass particle is assumed to be 29kN/m3. The comparison is made for the relation between the suction and the degree of saturation (Fig.11-12). The parameters of the theoretical model are shown (Table 1). It appears that the new model allows finding a correct estimation of the experimental curve on the drying path (cor.0.852-0.695) or on the wetting path (cor.0.843). The difference between the new theory and the experience at 50% of saturation is only 4.2-0.8kPa on drying path and 2.7kPa on wetting path  The simulation of the experience on a grade glass particle sample (Indarto, 1991) and on the Livet-Gavet loam is made along drainage and wetting paths. The sample are respectively graded between 4-140m and 1.6m-50mm of respectively unit weight 18,7kN/m² (e = 0.55) and 15,2kN/m² (e = 0.71)..The comparison is made for the relation between the suction and the degree of saturation ( Fig.13-14). The parameters of the theoretical model are shown (Table 1-2). It appears that the new model allows finding a correct estimation of the experiment curve on drying path (cor.0.954-0.875) and on the wetting path (cor.0.885-0.946). The difference between the new theory and the experience at 50% of saturation varies between 20 to 600kPa on drying path and between 5.6 to 50kPa on wetting path

Conclusion
A new theory is presented for the retention curve of the soil along drying and wetting path. It explain the difference between these two paths by the saturation of the soil particles with water menisci, which are independent at low saturation, then coalescent at high saturation on the saturation path and by the apparition of air nucleation followed by percolation on the drying path This new model uses only 5 physical parameters (mineralogy, void ratio, particle radius and granulometry, wetting angle, surface tension of water) and needs no additional experimental parameters always difficult to determine. It allows finding the retention curve along drying and wetting path.. It has been successfully used to determine the retention curve of uniform or graded glass sample and the natural Livet-Gavet loam.