Modelling of Static Liquefaction of Partially Saturated Non-Cohesive Soils

Static soil liquefaction is widely known to be a serious danger to the stability of structures. The phenomena governing pore water generation, which leads to liquefaction in fully saturated soils, are already quite well described. However, much less is known of these phenomena occurring in partially saturated porous media, although this, too, is an important issue in geotechnics. This study presents the application of a semi-empirical model to predict the response of partially saturated soils under undrained conditions. The model proposed is based on an incremental equation describing the pre-failure undrained response of partially saturated non-cohesive soils during monotonic shearing in a standard triaxial test. Improved differential equations taking into account pore fluid compressibility were implemented together with empirical coefficients describing soil skeleton compressibility during the unloading phase. Model coefficients were determined in triaxial compression tests. The influence of the saturation level represented by Skempton’s parameter B on the full spectrum of predicted stress paths was shown. For the analyzed saturation range, the maximum stress deviator normalized by initial mean effective stress varied from 0.38 to 1.67 for B values between 0.93 and 0.29, respectively. Model predictions were confronted with the results of triaxial tests for two types of non-cohesive soils (quartz medium sand and copper ore post-flotation industrial tailings). Good agreement between experimental data and theoretical predictions was achieved.


Introduction
It is well known that monotonically loaded non-cohesive saturated soils under undrained conditions may be prone to liquefaction due to a static [1] or cyclic [2,3] load. It was originally thought that, in order for soil to liquefy, voids should be completely filled with water which corresponds to full saturation of the medium. This was due to the fact that any presence of air would significantly limit or prevent the process of pore pressure build-up. Such an approach was physically justified because the presence of air results in full or partial dissipation of the pore pressure excess and consequently scales it down by the Skempton's parameter B decreasing liquefaction potential [4]. However, observations in areas particularly vulnerable to earthquakes have shown that liquefaction can also occur in soils that are not fully saturated. The first studies proving that partially saturated soils may also undergo liquefaction concerned the response of the soil to cyclic loads [5,6]. There followed many others related to the same issue, e.g., [7][8][9][10][11][12][13].
In order to theoretically predict the behaviour of partially saturated non-cohesive media, the MODSOL model proposed by Bian and Shahrour [20] has been elaborated, though it is dedicated mainly to cyclic loads as well. An interesting approach to modelling the response of soil to monotonic loading is the elasto-viscoplastic model proposed by Kimoto et al. [21], which was derived for media characterized by relatively high suction values and low saturation levels. Another theoretical study on the behaviour of unsaturated soils can be found in [22][23][24], in which the authors also generally focus on suction pressure at low saturation levels. The implementation of constitutive laws governing the behaviour of soil in a partially saturated state in a SaniSand-Z model built by Chen et al. can be found in [25]. There are some other theoretical approaches to constitutive modeling of partially saturated soils, e.g., a deviatoric hardening model which can predict static liquefaction proposed by Lü et al. [26], solutions based on thermodynamic equations by Yang et al. [27,28] or a hypo-plastic model of sand developed by Zapata-Medina et al. [29].
According to [4,30,31], in partially saturated media, the location of the instability line tends to increase with a decrease in the saturation level, which causes a change in the moment of triggering of the liquefaction process.
The problem of liquefaction of partially saturated soils seems to be important in many areas, also, from an engineering point of view, e.g., wave-inducted seabed liquefaction [32,33], blast-induced densification in liquefiable sediments [34], submarine largescale landslides triggered by gas dissociation [4], etc. Another interesting case relates to wet mine tailings storage facilities, where the zone of partial saturation could reach a dozen meters [35]. This is a consequence of the successive deposition of sediments in the form of soil-water mixtures by spigotting, which causes cyclic wetting and drying of stored waste mass [36]. As a result, relatively thick zones of partially saturated sediments appear above the full saturation line. In such a situation, the potential generation of pore pressure within this zone could, in an extreme case, lead to the liquefaction of tailing deposits, significantly reducing the stability of the surrounding dams. This is the reason why post-flotation cooper tailings have been selected as one of the investigated soils in this analysis.
In order to examine the influence of partial saturation on the undrained shear strength of soil, an extensive research programme has been carried out in recent years at the Institute of Hydro-Engineering of the Polish Academy of Sciences. The programme was aimed both at the empirical recognition and analysis of this phenomenon under laboratory conditions [37] as well as at its theoretical description [38].
First tests (Series A, [37]) were conducted on reconstructed specimens of post-flotation tailings (OZM50) from TSF Żelazny Most, the largest industrial copper tailings facility in Europe. Next, a complementary series of tests, the results of which are presented in this paper, was performed on the same soil (Series B), supplemented by tests carried out on medium quartz Skarpa sand (Series C). All tests were conducted under triaxial conditions. The first series of tests provided a basis for theoretical modelling of liquefaction phenomena of monotonically loaded partially saturated media [38]. The initial theoretical approach reproduced test results fairly well in the qualitative sense but quantitative agreement was insufficient, showing that the model needed some corrections. The present paper describes and discusses the results of the development and improvement of the theoretical model, reproducing the response of partially saturated non-cohesive soil subjected to monotonic shearing under undrained conditions in a standard triaxial test. In addition, verification of the model by laboratory tests on two different types of non-cohesive soils is also presented. The experimental results shown in the paper provide knowledge about unsaturated soil liquefaction mechanisms and constitute a basis for the theoretical description and prediction of partially saturated soil responses in triaxial conditions.

Generation of Pore Pressure in Partially Saturated Soils
The phenomenon of soil liquefaction is closely related to the volume changes which dry soil is subjected to during shearing (see the critical state theory [39][40][41][42][43][44]). In general, the behaviour of non-cohesive soil during shearing depends on its initial state, which determines its contractive or dilative response and/or contractive-dilative intermediate states [39]. Contractive and dilative states are separated by the critical state line (CSL) [39], which is also identified as a steady-state line (SSL) [45] at which sheared soil reaches its critical state [46] (Figure 1a). As a measure of the intensity of soil response (contractive or dilative behaviour), a state parameter [40] is usually used, which is defined as a change in void ratio that must occur for the medium to reach a critical state (steady state) at a constant mean effective stress.
Differences between undrained shear stress paths corresponding to each initial state (schematically shown in Figure 1b) are caused by changes in pore pressure build-up and subsequent redistribution between mean effective and total stresses. For initially contractive soils, this results in a decrease in the mean effective stress and subsequently a reduction of their shear strength. This reduction occurs after crossing an instability line [4,47], where the stress deviator reaches its maximum value. In the extreme case, effective stress vanishes, the soil is unable to transfer the shear stress and liquefies [48], which can be formally written as: The degree to which voids are filled with water may change between full saturation, represented by degree of saturation =1, and a completely dry medium (no water), =0. In between there occur intermediate states identified with unsaturated soils, which also include partial saturation. Partial saturation should be distinguished from an unsaturated state in which the gaseous phase of different levels is present, partial saturation being characterised by the lack of a continuous gas phase and the presence of only isolated air bubbles surrounded by water [18]. The model presented in this paper refers to partially saturated soils in which the saturation degree is usually above approximately 80% ( Figure  2). Furthermore, in non-cohesive media with a degree of saturation corresponding to a partial saturation state, suction pressure is negligibly small or does not occur at all [20,49,50]. In Figure 3, water retention curves for some selected non-cohesive soils, i.e., Hostun sand [20] and Singapore fine sand [51], as well as suction pressures for tested soils (determined in Sand Apparatus for the range of high saturation degree) together with their grain size distribution curves, are presented. It can be seen that for a degree of saturation above 0.8 we are dealing with very small suction pressure values of the order of 1 or 2 kPa and therefore the application of Terzaghi's law given by Equation (1) can be justified [52].  [20] and Singapore fine sand [51].
Soils in a state of partial saturation are characterized by values of Skempton's parameter B [53] in a range from = 1, which is equivalent to full saturation, to ≈ 0 , corresponding to the lower boundary of a partially saturated state.
Isotropically loaded saturated soil has a tendency to decrease in volume, but compaction under undrained conditions is prevented by bonds imposed by hardly compressible pore fluid, which leads to the generation of excess pore pressure, Δ . Skempton's parameter B is defined as a portion of generated pore pressure, , caused by corresponding isotropic stress. In an isotropic compression test in a triaxial apparatus, it is equivalent to the ratio of generated excess pore pressure to the change in cell pressure, Δ : Skempton's parameter B can also be described as a function of the compressibility of the pore fluid, , and the soil skeleton, [53]: where is soil porosity.

Materials and Methods
To examine the influence of saturation level on the liquefaction potential of soils, tests in a triaxial compression apparatus were carried out.
The first series of tests, described in [37] (Series A), and subsequent additional tests (Series B) were performed on anthropogenic post-flotation tailings, denoted as OZM50, which are sharp-edged fine sand grains with a rough surface ( Figure 4a). The third series (Series C) was carried out on a medium sand called Skarpa, which has rounded quarzitic grains with a smooth surface ( Figure 4b) [54].  The basic index properties and strength parameters of both soils are presented in Table 1; grain size distribution curves are shown in Figure 5. is specific gravity, is mean grain diameter, and are, respectively, minimum and maximum void ratios and is an internal friction angle for soils in loose state.

Figure 5.
Particle size distribution curves of tested soils.
In order to identify the initial state of particular samples of soils tested, critical state lines for both soils were independently determined based on standard triaxial tests (Figure 6): in drained and undrained tests for Skarpa sand [55] and undrained tests for OZM50 tailings [56].
The primary aim of the laboratory investigations was to empirically isolate the effect of saturation, expressed by the measured value of Skempton's parameter B on the response of soils during shearing in triaxial conditions. Therefore, in all tests performed, only the saturation level of samples was being changed in a controlled manner, whereas other parameters representing initial state conditions, for both soil types, were kept constant. Special attention was focused on the careful reconstitution of sandy samples in the triaxial apparatus, particularly by controlling void ratio from the very beginning of the specimens' installation as one of two parameters defining the initial state of non-cohesive soil. In order to fully control the value of the void ratio just prior to shearing, its changes during the first phases of the test (saturation, back pressure ramp, consolidation stages) were measured by local gauges making use of the Hall effect. A typical example of such changes is shown in Figure 7, representing the results obtained for test b1. All tested samples were in an initially contractive state. The samples were first subjected to isotropic consolidation by applying mean effective stress, = 400 kPa and = 200 kPa, for OZM50 tailings and Skarpa sand, respectively. All the tests were carried out in a strain-controlled mode at a constant strain rate, ̇≅ 10 % h ⁄ , and a constant cell pressure, = const. The saturation level of the soil was determined by Skempton's parameter B, measured during the test. In order to obtain assumed values of the parameter, back pressure was gradually increased in a controlled manner, which resulted in a decrease in the volume of air bubbles contained in soil voids, causing a reduction in pore air compressibility and thus a change in the compressibility of pore fluid as a mixture.
In total, thirteen experiments were carried out, covering a broad spectrum of saturation states represented by values of Skempton's parameter B between 0.29 and 0.93. In all cases, B was measured under a mean effective stress of 20 kPa after sample saturation and back pressure ramp but before consolidation. Such a level of mean effective stress was applied in each test for the sake of proper preparation of soil specimens to achieve the assumed initial void ratios as well as to install local gauges to control sample deformation prior to shearing. Table 2 contains the basic initial parameters of each experiment. is a back pressure, is mean effective stress at the start of shearing, is a relative density, is the state parameter, end are void ratios corresponding to an initial mean stress of 20 kPa and , respectively.
The location of the critical state lines for each soil in relation to initial void ratios just after sample formation and before applying vertical load are shown in Figure 8.

Experimental Results
Basic results of the laboratory investigations are shown in Figures 9-14. A detailed description and deeper analysis of the results from the first experimental programme (Series A) can be found in the study [37]. Detailed laboratory data for all the experiments analysed in this paper are enclosed in the online Supplementary Materials (Repository S1: Experimental Data).
Changes in pore pressure, normalized by initial mean effective stress in relation to vertical strain are shown in Figure 9 for Skarpa sand and in Figure 10 for OZM50 tailings, respectively. It can be seen that generated excess pore pressure increased with saturation and led to the liquefaction of samples with the highest saturation level.  In turn, changes of corresponding deviatoric stress in relation to vertical strain, , are presented in Figures 11 and 12 for Skarpa sand and for OZM50 tailings, respectively. It can be seen that lower saturation results in higher maximum and residual stress deviators.    In the presented research results, two stress paths corresponding to tests a6 and c4 differ from the rest. In both cases, failure of the sample took place before the boundary surface characteristic of each soil was reached. In these tests, unlike the others, clear shear surface was observed, typical for the deformation process in a highly compacted non-cohesive soil specimen.

Governing Equations
Based on the experimental results (Series A), a theoretical model was proposed to describe the behaviour of partially saturated non-cohesive soil in an initially contractive state subjected to axisymmetric monotonic loading under undrained conditions [38]. The core element of the model is a differential equation the integration of which allows the reproduction of any stress path within the stress space, including the generation of pore pressure in relation to the value of Skempton's parameter B characterizing the saturation level. A significant role in the proposed equation is played by two compressibility coefficients, the first corresponding to the compressibility of fluid treated as a mixture of incompressible water and compressible gas, the second to the compressibility of the soil skeleton. In general, these coefficients, which also depend on the type of loading (isotropic, deviatoric), are functions of soil state (void ratio and stress state expressed by mean effective stress and the stress deviator). The functions should be determined experimentally for each type of non-cohesive soil on the basis of isotropic consolidation and pure shearing tests under triaxial conditions.
The presented solutions are based on an idea of a semi-empirical incremental model for non-cohesive soils proposed originally by [57], which makes it possible to predict the liquefaction of fully saturated sandy soils, whereas the idea of developing a model for partially saturated soils was presented originally by Świdziński in an internal report [58] issued in 2015 and then published in the form of a research paper [38]. In general, the model can predict the full spectrum of soil behaviors shown in Figure 1 (from contractive to dilative ones); however, in this paper we focused on analyzing and theoretically reproducing the liquefaction phenomenon, therefore only the contractive state is analyzed. The basic assumptions of the model proposed are detailed below.
Volumetric strains of soil skeleton, , under triaxial compression conditions are dependent on mean effective stress and the stress deviator, .
or alternatively where is the stress ratio expressed by the following equation: = .
The stress ratio, , is a useful and convenient variable that makes it possible to define precisely the loading and unloading process in the stress space.
Functions describing strain-stress dependence (Equations (4) and (5)) take different forms for loading and unloading. Consequently, with a change in load direction, respective parameters have to be adopted.
Applying Equation (6), the change in the volumetric strain of the soil skeleton, , can be formally presented in the following form: Equation (9) is similar in its form to two incremental equations proposed by Sawicki [59] describing the development of volumetric and deviatoric deformations in dry sand through the superposition of the effects of stress tensor invariants: deviator and isotropic stress. Partial derivatives in Equation (9) can be replaced by corresponding functions (recall Equations (7) and (8)) and formally rewritten as: In the present case of partially saturated soils, it is assumed that the condition regarding the equality of potential changes in the volume of the soil skeleton structure and voids remains valid: where is the volume of pores and is the total volume of the soil skeleton structure. Using a definition of soil porosity in the form: the corresponding increments of pore fluid and soil skeleton volumetric strains d , d , respectively, are Combining Equations (11)-(13), we obtain: According to the definition of pore fluid compressibility , we can write: A general model equation describing changes in pore pressure caused by triaxial monotonic loading under undrained conditions is obtained by substituting Formulas (15) and (16) into Equation (10): Equation (17) contains three variables and cannot be solved by numerical methods. However, in the case of a standard triaxial compression stress path, when a constant cell pressure is maintained, we can easily eliminate one of the unknowns [38].
Stress deviator and mean total stress are defined by total principal stresses and as follows: and in the present case: Thus, the corresponding increments of stress invariants can be written as: and finally, we obtain: After transforming and differentiating Equation (6), we obtain: On the other hand, the transformation and differentiation of Terzaghi's fundamental equation linking effective and total stresses (Equation (1)) leads to the relationship: Next, after substituting Equations (23) and (24) into Equation (25), we obtain the dependence of the pore pressure increment on just two unknowns from Equation (17): Finally, substituting Equation (26) to Equation (17) and rearranging it, we obtain the following expression: Integration of Equation (27) makes it possible to reproduce a stress path in the ( ', ), stress space from which one can easily go to the ( ', ) stress space through Relationship (6). In the same manner, we can also solve Equation (17) for an arbitrary stress path defined by any -to-stress ratio.
The originally proposed theoretical model predictions [38], confronted with experimental data, yielded a very good qualitative agreement of the corresponding stress paths for the same B values, but a much worse quantitative fit, especially for low values of Skempton's parameter B, for which the model strongly underestimated the results of the experimental studies.
In order to improve the model and thus obtain better quantitative predictions for the whole range of parameter B, some modifications to the compressibility functions were made.

Pore Fluid Compressibility Function,
The volumetric compressibility of pore fluid, which is a mixture of water and air contained in voids, is defined as [60]: where and are, respectively, the compressibilities of air and water. Pore air can be treated as a perfect gas, so, assuming no temperature changes, its compressibility can be described by the formula: where = + is total pressure. Water is very poorly deformable; its compressibility under atmospheric pressure is = 0.45 × 10 Pa and does not show significant changes with its increase, so it can be treated as a material constant.
In laboratory tests, Skempton's parameter B was used to determine the saturation level of the soil because of its easy measurement during triaxial compression tests, in which a direct determination of saturation is extremely difficult. Since parameter B is not a physical measure nor is it useful in numerical calculations, Relationship (30) was used to determine the initial degree of saturation : Formula (30) was obtained by substituting Equation (3) into (28), whereas many authors use the inverse relationship in their studies, e.g., [61][62][63]. In the above form, the formula makes it possible to determine the saturation level on the basis of Skempton's parameter B, porosity and the isotropic compressibilities of soil skeleton, pore water and pore air. Please note that the effect of solubility has a significant influence on the compressibility of the water-air mixture [49] described by Equation (28). The approach presented to determine the degree of saturation based on an experimentally determined Skempton's parameter B also takes into account the change in pore fluid compressibility due to dissolved air.
In the course of a change in pore pressure, the volume of pores, , also changes, which consequently changes the degree of saturation. Based on the definition of the degree of saturation: = / and using Equations (13) and (16) for the pore voids and the analogical set for pore water, current soil saturation can be expressed in terms of the compressibility of pore fluid components according to the following equation: In a similar manner, by combining the definitions of soil porosity and compressibility of the pore fluid and keeping in mind the constant volume of the soil skeleton we can also easily determine changes in porosity caused by pore fluid compression during the test, which can be described as:

Soil Skeleton Compressibility Functions
In addition to soil porosity and pore fluid compressibility, the basic parameters of the model, there are two functions of soil skeleton compressibility, i.e., the isotropic function, , and the deviatoric function, , whose coefficients should be determined in triaxial tests for soils characterised by similar state parameters.
The function of deviatoric compressibility is determined in pure shearing tests under triaxial conditions during which mean effective stress is kept constant. For dry post-flotation copper tailings OZM50 in contractive state, volumetric strains as a function of stress ratio can be approximated by the following function ( Figure 15): where corresponds to the Coulomb-Mohr yield surface. Figure 15. Volumetric strains versus stress ratio during pure shearing under triaxial conditions for contractive OZM50 tailings.
The coefficients in Function (33), estimated by the least squares method for OZM50 tailings, are = 2.98 × 10 and = 3.11 × 10 . According to Equations (8) and (33), the deviatoric compressibility function can take the following form: In turn, the isotropic compressibility function can be determined in isotropic consolidation tests. Figure 16 shows volumetric strains that developed during isotropic compression applied to a specimen of the same post-flotation tailings OZM50. In this case, they can be nicely approximated by the logarithmic function: whose coefficients estimated by the method of least squares are = 2.97 × 10 and = 6.7 × 10 . The isotropic compressibility function is a derivative of the volumetric strains function, this time calculated with respect to mean effective stress: Under undrained conditions, the contractive response of sheared soil manifests itself in a reduction of mean effective stress as a result of pore pressure build-up. In the primary model, the reduction of mean effective stress which corresponds to unloading was not formally taken into account but it has quite a significant influence on its predictions. Therefore, it was necessary to differentiate between isotropic compressibility for loading and unloading and to introduce into the model the isotropic compressibility under unloading, , , which should also be determined by laboratory tests.
The results of an isotropic unloading test on contractive OZM50 tailings are shown in Figure 17. In order to approximate volumetric strains, which decrease during isotropic unloading, an exponential function was adopted:

Predicted Stress Paths
Before stress path integration, the following initial conditions have been defined: the degree of saturation calculated from the value of Skempton's parameter B (Equation (30)), the value of the stress ratio, (in case of isotropic consolidation = 0), porosity, mean effective stress and back pressure, with corresponding compressibilites of soil skeleton and pore fluid. Equation (27) can be solved numerically, assuming an incremental step of the stress ratio, Δ . For all cases presented in the paper the increment of Δ = 0.001 has been used, ensuring the stability of the solution. At the beginning of a single step increment, the mean effective stress corresponding to the current stress ratio is calculated by using Equation (27). Afterwards, the stress deviator and excess pore pressure are calculated by employing Equations (24) and (26). At the end of a given step, compressibilities, porosity and degree of saturation are updated.
Integration of Equation (27) for different values of the B parameter enables the coverage of stress space ( , ) by corresponding stress paths. Under partial saturation conditions (see Figure 2), stress paths are located between responses characteristic of a fully saturated soil and a completely dry soil. Figures 18 and 19 show the distribution of stress paths predicted by the model for different degrees of saturation and the corresponding Skempton's B parameters for OZM50 tailings and Skarpa sand, respectively.
Predicted changes of degree of saturation (Equation (31)) and porosity (Equation (32)) corresponding to the model simulations shown in Figures 18 and 19 were presented in Figures 20 and 21, respectively.

Verification of the Model
In order to verify the correctness of the model, its predictions were confronted with the results of laboratory tests (presented in Section 2.2).
The functions of soil skeleton compressibility were determined in the laboratory (for OZM50 tailings, see Section 3.3; for Skarpa sand, see studies [57,64]). All functions used together with their coefficients are collected in Table 3. The computations were made using the standard relationship between porosity, , and void ratio, , and Equations (28), (30) and (31), which relate pore fluid compressibility to Skempton's parameter B and pore pressure. Datasets with initial conditions for each test are given in Table 2.
Equation (27) was numerically integrated with the stress ratio increment Δ = 0.001. Predictions of stress paths obtained in this way could be easily transferred from the ( , ) coordinate system to the typical stress space ( , ) using Equation (6). The results of computations and laboratory tests are summarized in Figures 22 and 23.  For both types of soil, the model predicts stress paths quite well, both qualitatively and quantitatively. To quantify the prediction accuracy of the model, the maximum deviators obtained in the laboratory tests were compared with those obtained theoretically. For this purpose, the percentage measure of error, , was applied.
The results presented in Figure 24 show that for all cases the error was less than 16.1% and the average error was 7.7%, which in the domain of soil mechanics is quite satisfactory, although in the majority of tests analysed the model somewhat overestimated the test results.

Discussion and Conclusions
The results of laboratory tests presented in this paper have confirmed that the liquefaction phenomenon may also occur in partially saturated non-cohesive soils in a contractive state subjected to monotonic loading and that the potential for liquefaction decreases as saturation decreases. Furthermore, our investigations have shown that, under conditions of partial saturation characterized by positive values of Skempton's parameter B corresponding to the degree of saturation over 80% or higher, the behaviour of non-cohesive soil changes smoothly from that typical of a dry medium (or a medium under drained conditions) to that characteristic of a fully saturated medium.
Both experimental data and model predictions ( Figure 25) show the decrease of undrained shear strength with the increasing saturation level expressed by Skempton's parameter B. For the analyzed saturation range, the maximum stress deviator normalized by initial mean effective stress varies from 0.38 to 1.67 for B values between 0.93 and 0.29, respectively. Furthermore, linear dependency can be observed across the whole spectrum of partially saturated states. In order to reproduce the test results obtained, a modified theoretical model describing the process of excess pore pressure generation in partially saturated soils has been developed. The modification of the original model consisted in extending the formula for calculating the compressibility of pore fluid using soil skeleton compressibility functions based on dedicated laboratory tests and incorporating soil response under isotropic unloading to reflect the behaviour of soil during mean effective stress reduction due to an increase of pore water pressure. The proposed model reproduces the results of triaxial tests much better than the original one.
The model has been verified for two non-cohesive soils with different granulation characteristics, showing acceptable qualitative and quantitative agreement between predicted values and experimental results.
The research reported here offers a semi-empirical incremental approach to modelling the phenomenon of liquefaction or a significant weakening of partially saturated noncohesive soils under undrained conditions.