Modeling seismic compression of unsaturated soils in the funicular regime

. A semi-empirical elasto-plastic constitutive model with a hyperbolic stress-strain curve was developed with the goal of predicting the seismic compression of unsaturated sands in the funicular regime of the soil-water retention curve (SWRC) during undrained cyclic shearing. Using a flow rule derived from energy considerations, the evolution in plastic volumetric strain (seismic compression) was predicted from the plastic shear strains of the hysteretic hyperbolic stress-strain curve. The plastic volumetric strains are used to predict the changes in degree of saturation from phase relationships and changes in pore air pressure from Boyle’s and Henry’s laws. The degree of saturation was used to estimate cha nges in matric suction from the transient scanning paths of the SWRC. Changes in small-strain shear modulus estimated from changes in mean effective stress computed from the constant total stress and changes in pore air pressure, degree of saturation and matric suction, in turn affect the hyperbolic stress-strain curve’s shape and the evolution in plastic volumetric strain. The developments of the new mechanistic model developed in this study will play a key role in the future development of a holistic model for predicting the seismic compression across all regimes of the SWRC.


Introduction
During cyclic or seismic shearing of unsaturated soils, a reduction in total volume of the soil may occur as the particles rearrange into a tighter configuration and the void space is decreased, a phenomenon referred to as seismic compression. Seismic compression was defined by Stewart et al. [1] as the accrual of permanent contractive volumetric strains in soils during earthquakes and has been recognized as a major cause of seismically-induced damage in earth structures [2,3]. The current state of the practice in prediction of seismic compression during earthquakes is to empirically relate vertical strain profiles with depth to the shear strain induced in a soil layer as part of a site response analysis [4], or to use empirical charts to relate the cyclic stress ratio (CSR) during an earthquake with the standard penetration test results and initial degree of saturation in the field [5]. These are practical approaches as they build upon analyses performed for the earthquake site response but they do not necessarily consider hydromechanical coupling observed in laboratory experiments. However, even when evaluating laboratory experiments, there is also inconsistency in the expected trends in seismic compression with the initial degree of saturation from different experimental techniques having different pore fluid drainage conditions [e.g., 6,7,8]. Accordingly, a mechanistic framework is needed to confirm the trends in seismic compression in unsaturated soils and to understand the role of drainage conditions. * Corresponding author: mccartney@ucsd.edu The objective of this study is to develop a constitutive model that can represent the hydromechanical coupling expected in unsaturated soils undergoing seismic compression during cyclic or seismic shearing, with focus on the funicular regime of the soil-water retention curve (SWRC), where the initial degree of saturation exhibits the greatest variation with matric suction. The greatest amount of seismic compression is expected in the funicular regime due to the presence of both air and water. Also, undrained cyclic simple shear test data for sand specimens initially in the funicular regime of the SWRC are available for calibration and validation of the model [7,8]. This constitutive model is intended to be simple and based on an existing constitutive model (UBCSAND) so that properties can be calibrated with non-specialized laboratory tests and to facilitate implementation into software packages used by practitioners in geotechnical earthquake engineering (e.g., OpenSees, Plaxis, FLAC, etc.). This paper presents the details of the model formulation including the prediction of key variables. The model was calibrated using experimental shear stress-strain backbone curves from drained cyclic simple shear tests and transient SWRC scanning path measurements from undrained cyclic simple shear tests. Then the model predictions were validated using experimental data from undrained cyclic simple shear tests on unsaturated sand specimens with different initial degrees of saturation in the funicular regime reported by Rong and McCartney [8].

Model formulation 2.1 Effective stress-based elastoplastic model
The first step in developing the semi-empirical constitutive model for seismic compression is to define a backbone curve that represents the shear stress-strain curve if the soil were sheared monotonically to failure. Following the UBCSAND model, this study assumes that the backbone curve has the shape of a hyperbola. The equation for a hyperbolic shear stress-strain curve, which provides the shear stress for any applied shear strain , is given as follows [9,10]: where G i is the initial shear modulus (or small-strain shear modulus), τ ult is the ultimate shear stress from the shape of the hyperbola, and R f is a reduction factor that is used to adjust the shape of curve to the actual shear stress at failure f observed in an experiment according to the peak friction angle ' p . G i can be expressed using the following relationship based on the model of Hardin and Black [11]: where and n e are fitting parameters, P a is the atmospheric pressure used for normalization, and σ′ is the mean effective stress, defined using the approach of Lu et al. [12] as follows: where S e is the effective saturation, m is the mean total stress, u a is the pore air pressure, and u w is the pore air pressure. The first term in brackets is the net mean stress while the second term in brackets is the matric suction. For a cyclic simple shear test, the mean total stress is constant and can be calculated from the applied vertical stress v assuming an at-rest stress state as follows: where the coefficient of lateral earth pressure at rest, K 0 , which can be defined using elasticity (K o = / (1-), where = Poisson's ratio). An advantage of using the form of effective stress in Eq. 3 is that a SWRC model relating the effective saturation and matric suction can be incorporated into the definition of the effective stress. Another advantage of Eq. 3 is that independent predictions in u a , u w , and S e during cyclic shearing can be incorporated into the effective stress and the initial shear modulus. In the UBCSAND model, an increment in plastic shear strain can be calculated directly from the equation for the hyperbolic stress strain curve, as follows [9]: where d is the increment in developed stress ratio ( d ' m ) calculated from the derivative of the stress strain curve, as follows: and G p is the plastic shear modulus that can also be calculated from the stress-strain curve as follows: In UBCSAND, a flow rule is used to calculate the plastic volumetric strain as follows: where ' cv is the constant volume friction angle. While initial shearing from the origin will follow the hyperbolic curve in Eq. (1), during cyclic shearing, the shape of the hyperbolic curve is assumed to follow a hysteresis loop described by the following pair of equations [13]: where c and c are the cyclic shear stress and strain amplitudes, respectively. The backbone curve intersects the hysteresis loop at ±( c , c ). An example of a hyperbolic shear stress-strain curve for the first N = 0.25 cycles of shearing along with a hysteretic unloading and reloading path up to N = 1.25 is shown in Fig. 1(a). Plastic shear strains are not generated throughout the entire hysteretic stress-strain curve. During initial loading on the backbone curve, plastic shear strains are generated then during unloading back to zero shear stress it is assumed that the soil is elastic. However, plastic shear strains are generated when shear stresses in the opposite direction are applied. This is important as plastic volumetric strains are only calculated using Eq. 8 for the portions of the hysteretic loop where plastic shear strains are generated as shown in Fig. 1(b). To calculate the cumulative plastic volumetric strain ( ), the increment of plastic volumetric strain during each cycle is added to the value from the previous cycle ( ) as follows: = + (11) Also, it is assumed that the elastic volumetric strain in a cyclic simple shear test on an unsaturated sand in the funicular regime is approximately equal to zero as Rong and McCartney [8] found that changes in mean effective stress were within 10% of the mean total stress.

Degree of saturation model
It is assumed that there will be no change between the initial and final volume of water in an undrained specimen during cyclic simple shearing due to the high bulk modulus of water. Saturation is the ratio between water volume and void volume (i.e., void volume is comprised of the total water volume and air volume), therefore the change in saturation due to plastic volumetric strain (seismic compression) is linked to the volume change of the pore air voids. By incorporating the pore air volume change in terms of the plastic volumetric strain and initial void ratio, the following relationship for the evolution in degree of saturation during cyclic simple shearing: (12) where S o is the initial degree of saturation, e o is the initial void ratio, and V s is the unit volume of solids which can be assumed equal to 1 as the relationship between the different variables in Eq. 12 does not depend on the quantity of material. The derivation of this equation is given in Kinikles and McCartney [14]. As the volume of solids, initial void ratio, and initial degree of saturation are constants, the form of Eq. 12 indicates that the degree of saturation is inversely related to the plastic volumetric strain.

Pore air pressure model
Seid-Karbasi and Byrne [15] adapted UBCSAND to consider the seismic response of unsaturated soils by defining an equivalent bulk modulus of the pore fluids, and assuming that the generation of air and water pressures are equal during cyclic shearing. However, several experimental studies [8,16] found that the pore air and water pressures evolved in different manners during undrained cyclic shearing. Accordingly, separate equations are derived for pore air and water pressure generation during cyclic shearing. The instantaneous bulk modulus of air is assumed to be initially at atmospheric conditions and increase incrementally during undrained cyclic simple shearing. The degrees of saturations considered in this study are in the funicular regime, so the pore air voids are assumed to be continuous and connected. Some voids may be water filled, air filled or have partial air and water menisci. To extend the hyperbolic model to include a poromechanical approach, it is assumed that all volumetric strains correspond to a reduction in air void volume due to the incompressibility of water, the assumption of undrained conditions, and the assumption that no particle compression or crushing occurs. Air is assumed to be an ideal gas and the temperature of the air-water-soil system is at ambient temperature conditions (e.g., T = 294 K) with no fluctuation. Kinikles and McCartney [14] incorporated volumetric definitions of the phase diagram for unsaturated soils into the Ideal Gas Law, Boyle's law and Henry's Law to derive the pore air pressure change corresponding to the accumulation of plastic volumetric strains during undrained cyclic simple shearing. This relationship was found to be a function of the initial degree of saturation and initial void ratio, as follows: (13) where R is the universal gas constant (8.314 Nm/molK), h is Henry's constant, and T is the temperature in K. This pore air pressure equation is an isothermal equation that is comprised of constant values (initial void ratio, initial degree of saturation, Henry's constant, and temperature) and the plastic volumetric strain that will evolve during cyclic shearing according to Eq. 11. Data from Youd [17] indicates that all soils are expected to have compressive plastic volumetric strains occurring during cyclic simple shearing regardless of their initial density, which implies that that air pressure calculated from Eq. 13 will increase during a cyclic simple shear test. Further, the pressurization of the pore air during cyclic simple shearing is expected to affect the seismic compression in unsaturated sands during earthquakes by leading to a decrease in the mean effective stress according to Equation 2.2.2.

Pore water pressure model
While the generation in pore air pressure is assumed to be due to the compression of the air-filled voids during a reduction in total volume of the soil, it is assumed that generation of pore water pressure is due to shearinduced compression of the water-filled voids associated with relative particle movement and rearrangement. However, it is difficult to mathematically divide the energy from cyclic shearing into the independent generation of pore air and pore water pressure as this likely depends on the initial degree of saturation (the relative amounts of air and water in the soil) as well as the distribution in the two phases throughout a soil layer. Accordingly, a simplified approach was followed in this study to estimate the pore water pressure from the changes in degree of saturation estimated during cyclic shearing using Eq. 12. Specifically, the approach followed in this study is to estimate the changes in matric suction from the changes in degree of saturation, then to use the definition of the matric suction to calculate the pore water pressure from the pore air pressure calculated in Eq. 13. Specifically, the pore water pressure can be calculated from the pore air pressure and matric suction as follows: Use of this equation implies that the water retention in the soil is due primarily to capillarity. It is assumed that the thermodynamic state of the soil is constant throughout the duration of cyclic simple shearing as surface tension and capillarity are affected by temperature.
During undrained cyclic shearing, Rong and McCartney [8] found that unsaturated sand specimens initially on the primary drying path of the SWRC followed a wetting transient SWRC scanning path, as shown in Fig. 2. The transient scanning path links the initial state of a soil specimen on the primary drying path of the SWRC and the primary wetting path. The data of Rong and McCartney [8] also indicate that the slope of the transient SWRC scanning path during undrained cyclic shearing is dependent on the initial degree of saturation. Using the relationship between degree of saturation and matric suction, Kinikles and McCartney [14] proposed a log-linear relationship between the saturation and matric suction, as follows: where S o and o denote the initial position on the primary drying path of the SWRC, and m is the slope of the transient SWRC scanning path. This equation is valid for initial degrees of saturation on the SWRC in the funicular zone where there is a clear transition between the primary drying and wetting paths. For higher initial degrees of saturation (greater than approximately S o = 0.6), it is likely that the changes in pore water pressure and pore air pressure are approximately equal. In this case, a reconsolidation analysis using pore water pressures predicted using the model of Seid-Karbasi and Byrne [15] or the semiempirical approach of Ghayoomi et al. [5] could be used to estimate post-cyclic shearing volume changes.

Coupling in model
For strain-controlled conditions like those used in cyclic simple shear tests, increments of shear strain are inputted into the model and the shear stress is calculated using Eq. 1 with an initial shear modulus corresponding to the initial mean effective stress. For shear stresses calculated on the path of the hysteresis loop described by Eq. 9 and 10, the applied shear strain is assumed to be a plastic shear strain depending on the segments in Fig. 1(a). New values of volumetric strain, degree of saturation, pore air pressure, matric suction changes are then predicted using Eq. 11, 12, 13, and 15, respectively, a new value of mean effective stress was computed 3. In all cases considered in this study, the mean effective stress decreased during undrained shearing, which led to a softening of the initial shear strain modulus used in subsequent cycles. Each component of the model has a different impact on the evolution in plastic volumetric strains (or seismic compression), and all are dependent on the initial degree of saturation.

Model calibration
The model was calibrated using undrained cyclic simple shear experiments performed on well-graded sand with initial degrees of saturation in the funicular regime from Rong and McCartney [8] and triaxial tests performed by Zheng et al. [18]. The grain size distribution and SWRC of the sand are shown in Figs. 3(a) and 3(b), respectively. The sand classifies as SW according to the Unified Soil Classification System. The parameters of the van Genuchten [19] SWRC relationship are shown in Fig. 3(b) for the wetting and drying paths. The backbone curve was calibrated using the data from drained cyclic simple shear tests at different cyclic shear strain amplitudes reported by Rong and McCartney [8], and the predicted hysteresis loops match those from the experiments relatively well.    [14], and found that the soils with different initial degrees of saturation had a similar initial value of m = 0.021 at the beginning of cyclic shearing as shown in Fig. 4(a) and Fig. 4(b). However, the soils with higher initial degrees of saturation showed an increase in the slope as cyclic shearing continued. Kinikles and McCartney [14] proposed an empirical relationship specific to the SW sand tested by Rong and McCartney [8] describing the increase in m with the initial degree of saturation on the primary drying path as follows: = 0.021 + 0.117( − 0.117) (16) This empirical equation for m is only valid for initial degrees of saturation in the funicular regime between approximately S o of 0.117 and 0.6 for the sand tested by Rong and McCartney (2020b). Although Eq. 16 provides a simple and practical approach to consider the effects of the initial conditions on the evolution in the slope of the transient SWRC scanning path during undrained cyclic shearing, a more general theory is necessary in future work. This is especially the case as the data of Rong and McCartney (2020b) indicate that the transient SWRC scanning path is only log-linear as assumed in Eq. 15 for lower degrees of saturation, but becomes nonlinear for greater initial degrees of saturation as shown in Fig. 4(a).  A summary of the calibrated parameters of the model is shown in Table 1. This table does not include the parameters of the SWRC in Fig. 2(b) or those in the empirical relationship for m in Eq. 16, but they can be considered as additional parameters that need to be calibrated for a given soil. Testing under further conditions may be necessary to confirm if the parameters in Table 1 represent other conditions.

Evaluation of model predictions
The model was used to simulate the cyclic simple shear tests of Rong and McCartney [8], which were all performed with a constant cyclic shear strain amplitude of 1% applied at a frequency of 1 Hz on soils having different initial degrees of saturation as summarized in Fig. 2(b) but having the same initial void ratio of 0.636, and initial vertical total stress of 50 kPa. A sample of the evolution in the different variables predicted from the model as a function of the number of cycles together with the results from two duplicate experiments having an initial degree of saturation of S o = 0.3 is shown in Fig. 6. A positive aspect of the model is that it has good predictions of the evolution in all different variables until approximately N = 50 cycles. However, the predicted volumetric strains from the model tend to level off for larger numbers of cycles, while the experimental data do not show an asymptotic trend. The lack of an asymptotic trend in seismic volume change of sands was also observed by Youd [17] up to 10's of thousands of cycles. The shape of the relationship between the volumetric strain versus number of cycles from the model is dependent on a number of variables, including the shape of the calibrated hyperbolic stress-strain curve shown in Fig. 4. It is likely that the asymptotic trend in the volumetric change affected the predictions of the model at larger numbers of cycles. Nonetheless, the predictions of the model are good for this initial degree of saturation. Although this comparison is focused on cyclic shearing under a constant cyclic shear strain amplitude, real earthquakes may not have as many cycles of shear strain that are above an amplitude leading to seismic compression, so the asymptotic trend in volumetric strain with number cycles may not be a critical fault when simulating earthquakes. A comparison between the predicted volumetric strains predicted from the model with two duplicate experiments for different initial degrees of saturation is shown in Fig. 7. Similarly, comparisons between the predicted degree of saturation, matric suction, and mean effective stress are shown in Figs. 8(a), 8(b), and 9, respectively. While it is not possible to evaluate the impact of the initial degree of saturation from these figures, they show that the model is able to capture the general trends and magnitudes in the hydromechanical variables with cycles. An asymptotic trend is noted in all the model predictions but not always in the data.     [14]. While the model captured the coupled evolution in hydromechanical variables (pore air pressure, pore water pressure, matric suction, degree of saturation, volumetric strain, effective stress, shear modulus) have a reasonable match after the first 15 cycles of shearing, the predictions do not match the trends in the data after 200 cycles of shearing. Scatter in the experimental data was also observed, partly because preliminary tests not reported by Rong and McCartney [8] were included in the comparison as described by Kinikles and McCartney [14].
After 200 cycles of undrained shearing, a linear decreasing trend between seismic compression and initial degree of saturation was predicted from the model while a nonlinear increasing-decreasing trend was observed in the cyclic simple shear data. A similar linear decreasing trend in seismic compression was noted by Ghayoomi et al. [13] in their seismic compression data, but their experiments were not fully undrained like those of Rong and McCartney [2020b]. Accordingly, this discrepancy may be due to issues with the model in capturing all mechanisms of seismic compression in undrained shearing, calibration of model parameters, or drift in the position of the experimental hysteretic shearstress strain curves. For example, the experimental transient SWRC scanning paths and those predicted by the model are shown in Fig. 11. While the slopes of the curves matched well initially, the range of the predicted changes in degree of saturation decreased with increasing initial degree of saturation. This may indicate that the model used for the transient SWRC scanning paths may need to be refined, or that the assumption of linking the degree of saturation to the shear-induced pore water pressure generation may not be valid for higher degrees of saturation. Another interesting observation was that the initial degree of saturation had a significant effect on the evolution in the pore air pressure with the number of cycles that matched the measured values well [14]. However, after many cycles the pore air pressures approached an asymptote that was nearly the same for all specimens as shown in Fig. 12. This is most likely due to the asymptotic shape of the volumetric strain curve, which drives the evolution in the other variables, but it also indicates that the pore air pressure becomes less influential at large numbers of cycles.

Conclusions and final comments
A new mechanistic model was presented in this paper that may be useful for predicting the evolution in seismic compression for different initial conditions as a function of cycles of shearing. While discrepancies were noted between the model predictions and experimental measurements from undrained cyclic simple shear tests, the model represents a step forward in the ability to predict the coupling between different hydromechanical variables affecting seismic compression. The model components will play a key role in the development of a holistic model for predicting seismic compression across all SWRC regimes, where reconsolidation of excess pore water pressures must also be considered. Insights into future pathways for improving the model were noted, but it is also important that additional experiments be performed with different vertical total stresses, relative densities, and degrees of saturation to fully calibrate and validate the model.