Modeling the soil-water retention curves for highly deforming soils

A porous-solid model based on the grain and pore size distributions of the soil is coupled with a mechanical model to simulate the soil-water retention curves while the material is deforming. During the determination of the main drying curve, the soil is subjected to high suctions which induce important volumetric deformations. These volumetric deformations modify the pore size distribution of the sample affecting both the drying and the wetting retention curves. Although, most deformation occurs at drying, the drying curve is only slightly affected by soil deformation. In contrast, the wetting curve shows important shifting when volume change is considered.


Introduction
Different models have been proposed to account for the density of soils during the determination of the soil-water retention curves (SWRCs). See for example [1][2][3][4][5][6][7][8][9][10][11][12], among others. The main drawback of these methods is that none of them considers full hydro-mechanical coupling. When an increment in suction is applied, the soil deforms and the SWRC modifies. This in turn affects the value of the effective stress parameter ߯ used to calculate the next volumetric strain. As Della Vecchia [7] points out, the resulting SWRC is not the same when the retention curve is simulated using the final pore size distribution (PSD). In addition, none of these models includes simultaneously the case for double structured soils along with the simulation of scanning curves during wetting-drying cycles.
In this paper, a porous-solid model based on the grain size distribution (GSD) and current PSD of the soil, is used to simulate the evolution of retention curves while the soil is deforming. This porous-solid model is coupled with an effective stress mechanical model to determine the volumetric strains with suction. The results of the proposed model are compared with some experimental results and some conclusions are withdrawn.

The porous-solid model
In order to correctly simulate the retention curves of deforming materials, hydro-mechanical coupling is required. A simple way to generate coupled models for unsaturated soils, is to build up a hydraulic model based on the current PSD of the material. A porous-solid model able to simulate the retention curves during wettingdrying cycles has been presented elsewhere [13,14,15] In this section, a porous-solid model able to simulate the retention curves during wetting-drying cycles is briefly presented. This model considers two types of pores: the cavities (C) and the throats or bonds (B). The former are those elements containing most of the volume of water of the soil and are represented by circles or spheres in two or three-dimensional networks, respectively. The bonds are those pores interconnecting the cavities and are represented by rectangles or cylinders in two or three dimensional models, respectively. Fig. 1 shows the porous-solid model including the solids (S) in addition to the cavities and bonds. All these elements have been placed at random in a regular twodimensional network. The random distribution of pores of different sizes guaranties the correct simulation of the phenomenon of hysteresis where the extension of the overlap between the PSDs of cavities and bonds plays a fundamental role. However, during the set-up of the network, a construction principle needs to be verified. This construction principle states that two bonds converging at 90 0 at a cavity should not intersect each other. In the nodes where the construction principle is not verified, a redistribution of cavities and bonds is performed until the principle is fulfilled everywhere in the network. The representation of the network in Fig. 1 is only schematic as the length of bonds in the porous-solid model is only a fraction of their diameter. This means that in fact, the resulting network is highly irregular.
Suppose that the network in Fig. 1 is initially subjected to a large suction, all pores are dry and water is available at the boundaries. For the modelling, it can be considered that all boundaries are connected to the bulk of water although other conditions can be imposed. Then, suction is progressively reduced by steps. When this happens, some pores at the boundaries start saturating. During wetting, the first pores to saturate are the smallest [16]. The sizes of pore that can saturate are obtained from the Young-Laplace equation where ‫ݎ‬ represents the critical radius (i.e., the radius of the largest pores that can be saturated at suction s), ܶ ௦ is the air-water interface tensional force, D is the contact angle between water and soil particles (usually considered as zero). Consider a critical size, meaning that all pores from sizes 1 to critical can saturate at this stage [14]. However, pores only saturate if they are connected to the wetting boundaries by a continuous path of saturated elements. In Figure 1, saturated pores are shown in black while dry pores are shown in white. Additionally, solids are represented in gray tones. Using these rules, it is possible to determine which pores saturate or dry during a wetting or drying processes and thus, obtain the main SWRCs.
Mercury intrusion porosimetry tests (MIPTs) indicate that real soils show PSDs that can be approximated to logarithmic normal distributions. This means that only two parameters, the mean size (ܴ ത ) and standard deviation (߫), are required to describe these distributions. One distribution is used for bonds and the other for cavities. Therefore, four parameters are required to describe a single structured soil and simulate both retention curves. For double structured soils, two different logarithmic distributions are required for both elements (bonds and cavities): one for the macro and the other for the micropores. In addition, a parameter relating the relative volume (߷) between macro and micropores is required. Because there is a one to one relationship between PSD and main SWRCs, when the former is not available, use is made of the last. In such a case, by successively modifying an initially proposed PSD, it is possible to fit the numerical with the experimental SWRCs and, finally obtain a theoretical PSD of the material. Both the wetting and the drying curves are required in this process, in order to determine all parameters of the model. Because the shape of pores in real soils is not spherical, it can be said that the theoretical PSD obtained after the fitting process represents an equivalent PSD of the real soil. Therefore, comparisons between experimental and numerical PSDs may show some differences.
Usually, porous networks with at least 1 million nodes are required to avoid the influence of the size of the network on the SWRCs. However, as the size of the network increases, the computing time and storage requirements represent two main concerns. In order to avoid these problems, probabilistic porous-solid models have been developed. These type of models are based on the same principles of the network models with the advantage of simulating infinite porous-solid networks. Probabilistic models greatly reduce the computing time and storage requirements and simplify the fitting process of the experimental SWRCs [15].
The effect of volumetric strains on the PSD can be easily included in the porous-solid model. According to the experimental results reported by Simms and Yanful [17], only the volume of macropores reduces when a soil shows plastic volumetric strains during compression. Therefore, these plastic volumetric strains represent the reduction in the volume of macrocavities as the volume of macrobonds is negligible, as stated before. This reduction in the volume of macrocavities can be included into the model in three different ways: by a reduction in the number of macrocavities, by a reduction in the mean size of macrocavities or by a combination of both. The experimental results reported by Simms and Yanful [17] indicate that soils subjected to suction increase show mainly a reduction in the mean size of macrocavities. Thus, this same approach has been included in the poroussolid model for the case of suction increase. Also, these experimental results show that the reduction in the mean size of macropores is approximately of one order of magnitude from the saturated to the fully dry condition. Therefore, by considering these assumptions, it is possible to determine the updated size of macrocavities that accounts for the current volume of the sample. For macrobonds, the same relative variation in their size is considered. This procedure neither requires any additional parameter for the model nor previous calibration.
Once the updated PSD of the soil has been established, the next point of the SWRC is obtained by simulating the SWRC from the initial suction to the current value of suction. This is so because the water in the soil, redistributes with each variation in the size of pores. In other words, each point of the SWRC is obtained from the updated PSD of the material generated by the current suction. This procedure allows determining the SWRCs while the material is deforming.

The mechanical model
The elastoplastic framework used to model the volumetric behavior of unsaturated soils is based on the effective stress concept [18]. This framework uses the effective mean stress ‫)′(‬ defined by Bishop´s equation S6 S7 S8 S9 S10 S11 S12 S13 S14 S15 S16 S17 S18 S19 S20 where ‫̅‬ is the net stress and ‫ݏ߯‬ is the suction stress.
Bishop's parameter ߯ is obtained from the relationship [19] where ݂ ௦ and ݂ ௨ represent the saturated and unsaturated fractions of the soil, respectively, while ܵ ௪ ௨ is the degree of saturation of the unsaturated fraction. The saturated fraction is obtained by adding the volume of all solids exclusively surrounded by saturated pores (ܸ ௦ ௦ ) with the volume of these pores (ܸ ௩ ௦ ) and divided by the total volume of the soil (ܸ), in the form ݂ ௦ = (ܸ ௦ ௦ + ܸ ௩ ௦ ) ܸ ⁄ . The unsaturated fraction represents the volume of solids surrounded by a combination of saturated and dry pores (ܸ ௦ ௨ ) added by the volume of these pores (ܸ ௩ ௨ ) and divided by the total volume of the soil, in the form Finally, the degree of saturation of the unsaturated fraction is the volume of saturated pores (ܸ ௩ ௨௦ ) divided by the total volume of pores appertaining to the unsaturated fraction of the sample, in the form ⁄ . For fully coupled models, these three parameters depend on the current value of suction and the whole wetting-drying and loading-unloading paths applied to the sample. Eventually, a dry fraction (݂ ௗ ) may appear at large suctions. The dry fraction is represented by the volume of solids exclusively surrounded by dry pores (ܸ ௦ ௗ ) added by the volume of these pores (ܸ ௩ ௗ ) and divided by the total volume of the soil, in the form With these definitions, the relationship ݂ ௦ + ݂ ௨ + ݂ ௗ = 1 is fulfilled. All these parameters (݂ ௦ , ݂ ௨ , ݂ ௗ , ܵ ௪ ௨ ) can be obtained from the porous-solid model described in the previous section.

The elastoplastic framework
The elastoplastic volumetric behavior of the soil is determined from the following equation where ݁ and ݀݁ represent the current void ratio and its variation, respectively, ‫′݀‬ is the increment of the mean effective stress while ߣ is the compression index in a logarithmic plane formed by the void ratio and the mean effective stress. Because the void ratio of the sample reduces with increments on the effective mean stress, ߣ shows negative values. To account for the elastic behavior, parameter ߣ is substituted by the unloadingreloading index ߢ . These parameters are single valued for any stress path and suction applied to the soil including loading-unloading and wetting-drying cycles. They can be obtained from loading-unloading compression tests on saturated samples Fig. 2 shows the elastoplastic framework for the volumetric behavior of soils in the planes of the logarithm of effective mean stress versus logarithm of void ratio ( Fig. 2(a)) and effective mean stress versus suction ( Fig.  2(b)). Consider a saturated normally consolidated soil sample subjected to a net stress ‫‬ (point A). The soil can follow a drying or loading path, but in any case, the slope of the compression line is given by the same parameter ߣ in Fig. 2(a).
Suppose that in this case, the sample follows a drying path up to suction ‫ݏ‬ reaching the saturated preconsolidation stress ‫‬ * ᇱ (point B) . Fig 2(b) shows the position of both the suction increase yield surface (SIYS) and the loading collapse yield surface (LCYS) at the end of drying. Observe that the LCYS is horizontally displaced from the drying path a quantity equal to the suction stress ߯ ‫ݏ‬ . This displacement is generated by the suction hardening phenomenon and is simulated by coupling the SIYS with the LCYS. This means that when loading is preceded by drying, at the beginning of the loading stage, the sample shows an initial elastic behavior up to the apparent preconsolidation stress ‫‬ ᇱ as shown in Fig. 2(a) (path BC). Once the apparent preconsolidation stress is surpassed, the sample shows elastoplastic strains as indicated in the same figure.

Modeling the SWRC
When the drying SWRCs is being determined, the soil sample is subjected to a continuous increase in suction. When the maximum previous suction of the sample is surpassed, volumetric plastic strains occur and the macropores of the soil reduce in size. Consequently, the SWRC is affected as well as the effective stress according to Eqs. (2) and (3). In order to assess the influence of the volumetric strains on the SWRCs, the porous-solid model can be used advantageously as it is based on the current pore size distribution of the material. The following procedure to simulate the SWRCs while the soil is deforming has been applied in this research. First, the initial conditions of the sample are established: the GSD, the initial PSD and void ratio, as well as the maximum and current net stress and suction applied to the sample. Then the elastic and elastoplastic compression indexes are determined from compression tests. Finally, the porous-solid model coupled with the mechanical model is used to simulate the SWRCs. The volumetric strains produced by suction increases are computed using the elastoplastic framework shown in Fig. 2 along with Eq. (4). The pore size distribution of the soil is continuously updated for each plastic volumetric strain. This means, that each point of the SWRC is obtained from a different PSD when the soil exhibits plastic deformations. Each time that the retention curve is updated, so is the value of parameter ߯. In this way, it is possible to correctly compute the value of the increment of suction stress (Δ(߯‫))ݏ‬ and the associated volumetric strain for the next increment of suction. In this manner, a fully coupled modeling of the SWRCs can be performed.

Tests by Gao et al. [8]
These authors, reported the SWRCs of a clayey silt (called Pearl clay) at three different initial void ratios. The soil is composed of 26% of clay and 74% of silt. This material shows a liquid limit of 43%, a plasticity index of 17.5% and specific gravity of 2.71. The samples were formed by static compaction at a water content of about 26%. The initial void ratios of compacted samples were 1.36, 1.24 and 1.1. Fig. 3(a) shows the experimental (Ex) GSD of this material. Also in this figure, a fitted numerical (N) GSD used to build up the porous-solid model is presented. This numerical GSD was obtained from two logarithmic normal distributions: one for large grains (MS) and the other for small particles (mS). These distributions adopted the following values for the mean size (ܴ ത ) and standard deviation (߫) for the large and small particles, respectively: ܴ ത ‫)ܵܯ(‬ = 2.3, ‫)ܵܯ(ߜ‬ = 2.5, ܴ ത (݉ܵ) = 0.2, ߜ(݉ܵ) = 2.0 while the relative volume of large solids ߷ = 0.02. These parameters take the same value for all soil samples and remain constant during the compression of the sample by suction increase.
For the determination of the SWRC, suction was controlled using both the axis translation and the vapor equilibrium technique. The axis translation technique was used for suctions up to 1.5 MPa and the vapor equilibrium for suctions up to 450 MPa. The experimental retention drying curves for the samples prepared at void ratios of 1.1 and 1.3 are shown in Fig. 3(b).
For the simulation of the retention curves and similar to GSD, two logarithmic normal distributions with their respective mean size, standard deviation and relative volume of macropores were used for both cavities and bonds. Prefixes macro (M) and micro (m) are used for the large and small elements of each type, respectively. By fitting the numerical with the experimental SWRCs for samples with initial void ratios of 1.1, 1.24 and 1.3, the same parameters for the mean size and standard deviation for both the microcavities and microbonds as well as for the relative volume of large pores were obtained, for all three different void ratios. These values are ܴ ത ‫)ܥ݉(‬ = 0.02   This means that only the mean size of macrocavities and macrobonds change for each sample. These values are shown in Table 1. This finding agrees with the experimental results shown in Fig. 3(c) which only shows variations in the volume of macropores for samples compacted at two different void ratios. Even if samples compacted at the same water content show different PSDs, this difference can be accounted for solely with the mean size of macropores. This is why the mean size and standard deviations for both microcavities and microbonds maintain a constant value for all samples.
It is possible that a relationship between the size of macrocavities and the void ratio of compacted samples exists. This possibility can be explored when more experimental results become available.
The relative volume of pores obtained from MIPTs for samples at void ratios of 1.15 and 1.35 are shown in Fig.  3(c). Also in this figure, the numerical PSDs obtained from the parameters indicated in Table 1 are plotted. It can be observed that the experimental and numerical mean sizes for macro and micropores show more or less similar values. In contrast, the numerical maximum relative volume for macro and micropores show smaller and larger values with respect to experimental results, respectively. Moreover, the experimental distributions show sharper curves than the numerical ones. This last characteristic is regulated by the standard deviation of macro and micropores. This means that this parameter should take smaller values than those indicated in Table 1, to ensure a better correlation with the experimental results. However, if the standard deviation is reduced, the SWRC may show a step, clearly indicating the presence of a double structured material. But this is not the case for the experimental results reported in Figure 3(b).
In Fig 3(d), the initial (I) and final (F) PSDs for cavities (C) and bonds (B) are shown. These distributions are for the sample with an initial void ratio of 1.32. The initial PSD was obtained from the parameters indicated in Table 1, while the final PSD was obtained at the end of drying. This figure also includes the GSD of the sample. Observe that, because large cavities and bonds reduce in size, the smaller cavities and bonds show an apparent increase in their relative volume [17]. Each time that the increase in suction produces plastic volumetric strains, the PSD slightly modifies from its initial condition. The next point in the SWRC is then obtained from the updated PSD. Then, parameter ߯ is also updated in order to correctly determine the current suction stress and, finally, the decrease in void ratio for the next increase in suction is computed from Eq. (4). Suction stress evolution during drying for samples with different initial void ratios. Experimental data after [8].
During the experimental determination of the drying curves, the volumetric strains of soil samples were measured for each increment in suction as reported by Gao et al. [8]. The experimental volumetric strains for the three different samples are shown in Fig. 4(a) along with the numerical simulation. Observe that when suction increases beyond 1 MPa, soil samples show a rebound on their volumetric strain. This behavior is generated by the reduction in the suction stress ‫)ݏ߯(‬ which takes place when suction surpasses 1 MPa, as shown in Fig 4(b). This figure shows the variation in suction stress for samples with initial void ratios of 1.1 and 1.32 during drying. Observe that numerical and experimental variations in void ratio with suction adjusts well for all three samples as shown in Fig. 4(a). Notice that experimental results show an erratic behavior when suction surpasses the maximum suction stress. This could be attributed to the vapor exchange technique used for suctions beyond 1.5 MPa.
In order to assess the influence of coupling the hydromechanical behavior of soils during the determination of the SWRCs, Fig. 5 has been prepared. This figure shows the numerical modeling of the main retention curves at wetting (W) and drying (D) using the coupled (Cu) and the uncoupled (Un) model. For the correct interpretation of these curves, it must be considered that bonds control mainly the drying curve while cavities control mainly the wetting curve [13]. In addition, macrobonds and macrocavities control the drying and wetting curves at low ranges of suction, respectively, while microbonds and microcavities control the drying and wetting curves at large ranges of suction, respectively. Also, it must be considered that only macropores change their size during compression while the size of micropores remains constant. Moreover, it must be considered that during drying, the PSD is changing with each plastic volumetric strain. This change in PSD only arrest once the suction stress reaches its maximum value as afterwards, it reduces. Instead, during wetting, the PSD remains constant as only elastic strains occur. Finally, the shift of the retention curves on the suction axis, depends on the plastic volumetric strain experienced by the soil during drying. In this way, it can be observed in Fig. 5 that the coupled and uncoupled drying curves initially superpose. However, when suction increases beyond 0.02 MPa (maximum previous suction of the soil samples), plastic volumetric strains appear and the coupled drying curve displaces to the right hand side on the suction axis. Finally, when the degree of saturation reduces below 0.2, both curves appear close together again. This is so because the last bit of the curve is regulated by the small bonds which do not change in size during compression. With respect to the wetting curve, it can be observed that the reduction in size of macrocavities after drying, produces the displacement of this curve from the beginning although, this displacement reduces as the degree of saturation increases. In this case, the displacement of both the coupled wetting and drying retention curves is small, because the maximum plastic volumetric strain is also small. Larger volumetric strains and therefore, larger displacements on the SWRCs are observed in the next set of tests.

Tests by Salager et al. [20].
For these tests, samples of a clayey sand (CS) statically compacted at densities ranging from 13.5 to 19.5 kN/m 3 and water contents from 18 to 12% were prepared. The soil was constituted by 72% sand, 18% silt and 10% clay. Compacted samples of this material show an optimum water content of 14.5 % and a maximum density of 18.6 kN/m 3 . The liquid and plasticity limits for this soil are 25 and 14.5%, respectively. After compaction, the samples were placed in a pressure plate chamber and saturated by imbibition. At the end of saturation, each sample reached the initial void ratio (݁ ), density (ߩ) and water content ‫)ݓ(‬ indicated in Table 2. For the determination of the SWRCs, the axis translation technique was employed in the range of 0.001 to 1 MPa, while the vapor equilibrium technique was used in the range from 4 to 326 MPa. The numerical and experimental main retention curves at drying for samples indicated in Table 2 are plotted in Fig. 6(a). The fitting parameters for the simulations of the SWRCs for each sample are presented in Table 3. Notice that the mean size of microcavities mC = 0.00025 Pm and microbonds mB = 0.00015 Pm as well as the standard deviation for both cavities and bonds ߜ = 5.0 remains constant for all samples. Therefore, Table 3 only shows the values of the mean size of macrocavities and macrobonds as well as the relative volume between large and small cavities (߷).  Notice that in these tests, both the void ratio and water content change for each sample. The PSDs obtained from the data in this table are depicted in Fig. 6(b). Observe how the mean size and relative volume of macrocavities reduce with the void ratio while the microcavities show an apparent increase in relative volume. The experimental GSD reported by the authors was fitted using the following parameters: ܴ ത ‫)ܲܯ(‬ = 130 Pm, ‫)ܲܯ(ߜ‬ = 3.5, ܴ ത (݉ܲ) = 0.2 Pm, ߜ(݉ܲ) = 3.0 and ߷ = 0.000 0055. These parameters remained constant during the simulations.   For the sample with a void ratio of 1.01, the initial and final PSDs for both cavities and bonds are shown in Fig.  6(c). These variations in the initial and final PSDs are generated by the volumetric strains of the sample as suction increases. In Fig 7(a), both the experimental and the numerical volumetric strains for samples with different initial void ratios are indicated. Figure 7(b) shows the values of suction stress during drying. Observe that a maximum value of suction stress occurs for all samples at around 30 MPa of suction. Therefore, when suction surpasses this value, a small rebound on volumetric strains occurs as observed in Fig 7(a). Fig. 8 shows the comparisons of the main SWRCs for a sample with initial void ratio of 1.01, when simulations are made with coupled and uncoupled models. Observe again that the displacement of the drying curve is only partial while the wetting curve shows a considerable displacement from the beginning although, it reduces as saturation is approached. These results differ from other hydro-mechanical coupled models where the air entry value of the SWRCs is the only parameter affected by volumetric strains resulting in parallel retention curves.

Conclusions
A coupled hydro-mechanical model is used to simulate the SWRCs of soils. The hydraulic model is based on the current PSD of the material which is continuously updated as the soil deforms. Also, parameter ߯ and the mean effective stress are continuously updated to correctly determine the volumetric strain for the next increment in suction.
The numerical and experimental comparisons show the appropriateness of the proposed model. Numerical results show that the drying retention curve suffers a partial shift on the suction axis when compared with nondeforming simulations, even if large plastic volumetric strains occur at this stage. In contrast, the wetting branch shows an important shift from the beginning as it initiates with a modified PSD. However, this shift reduces as the saturated condition is approached. These results differ from other models where the air entry value is the only variable affecting both curves during hydro-mechanical coupling.
The proposed methodology considers that only the number of macropores reduce with plastic volumetric strains, although, it is possible that the size of pores also reduces. This represents a limitation for the procedure proposed herein.