Physical Modeling of Liquefaction in Various Granular Materials

Abstract

This paper compiles numerous experiences gained from physical models, to highlight the phenomena of triggering, propagation, and mitigation of liquefaction in granular soils. Results of tests at different scales, from the element volume (cyclic triaxial tests) to small-scale models in centrifuge, performed using several granular soils, will be presented to provide behavioral tools for predicting the phenomenon. Furthermore, the efficacy of vertical and horizontal drains as liquefaction mitigation techniques will be discussed.

1. Introduction

Cyclic liquefaction is a sudden phenomenon of loss of shear strength and stiffness, which may occur when a granular and non-plastic soil, whose voids are saturated by an incompressible fluid, is vibrated at a frequency too high for the soil to comply with its tendency to contract discharging pore water.

The tendency to contract produced by cyclic and dynamic shear strain, as those induced by an earthquake, turns into accumulation of excess pore Δu. As Δu rises, hydraulic gradients establish within a deposit, triggering fluid flow from higher toward lower hydraulic load. If the tendency to dissipate excess pore pressure is overcome by the tendency to accumulate pore pressure, the contact pressure between grains may cyclically reset, zeroing shear strength and stiffness, and the soil starts to behave like a viscous fluid.

The effects on existing buildings or structures may range from settlement and tilting up to catastrophic failures.

This chapter deals with earthquake-induced liquefaction and its modeling, from the element volume scale to the small-scale physical models.

The first part of the chapter describes some experimental observations gained testing a natural soil that experienced liquefaction during the 2012 Emilia seismic sequence, in Italy. The most relevant liquefaction manifestations were observed during the May 20 shake in the Ferrara Province. The moment magnitude was M_w = 6.1, with an estimated peak ground acceleration PGA ≈ 0.26 g. The affected sites, located about 15 km SE of the epicenter (estimated PGA ≈ 0.16 g), exhibited various phenomena, including craters, sand boils, surface cracks, and lateral spreading in free field areas. Additionally, there were reports of moderate building settlement and tilting at developed sites. The sandy stratum which underwent liquefaction resulted from the fluvial processes of the Apennine Reno river during the period spanning from 1450 to 1770. The sand within this layer is normally consolidated and loose. This sand was tested monotonically and cyclically to find a relationship between its state, and the cyclic resistance. In addition, it was noticed that sandy deposits of similar origin and age located in regions closer to the epicenter of the earthquake on May 20, 2012, did not experience liquefaction, suggesting the hypothesis of potential partial saturation of these soils. This supposition gains support from the recurring reports of gas emissions from the soil and the presence of gas within the groundwater, particularly within a few kilometers of the earthquake’s epicenter. Therefore, cyclic triaxial tests were conducted on both saturated and unsaturated samples to determine the increase in liquefaction resistance resulting from partial saturation. P-wave velocity was assumed as a measurable variable to estimate the cyclic resistance ratio, CRR of partially saturated sand.

The area where extensive liquefaction occurred in 2012 was assumed as one of the case studies of the LIQUEFACT project (http://www.liquefact.eu/) and the ground conditions at those sites were taken as a reference for a large series of centrifuge tests, which are partly described in the second part of the chapter. The primary objective of the centrifuge tests was to investigate the seismic response of saturated sandy deposits when subjected to progressively increasing dynamic actions, ultimately leading to liquefaction. Concurrently, the campaign aimed to assess the efficacy of various mitigation strategies for liquefaction.

At the end of the chapter, a simplified method of liquefaction assessment is proposed, based on the analysis of more than 60 centrifuge cone penetration tests CPTs on sandy soils and on the results of cyclic laboratory tests.

2. Liquefaction, a case study in Italy

The onset of liquefaction depends on the cyclic shear loading generated by an earthquake and the cyclic resistance of the soil. The latter is influenced by various factors, including:

• Grain size and mineralogy (plasticity)

• State of the soil (void ratio and confinement stress)

• Degree of saturation

• Aging, bonding, structure

• Ground conditions (level ground, sloping ground, superficial loads)

Manifestations of liquefaction may range from sand boils, craters, and fissures in free field conditions to the floating of buried structures, building settlement and rotations, and lateral spreading of sloping areas. Figure 1 shows, as an example, liquefaction evidences recorded at the site of San Carlo, in Italy, where liquefaction occurred during the 2012 Emilia seismic sequence [1]. Liquefied sand flooded out and submerged large areas of the village, surface ruptures, extensional fissures, sand boils, vents, sinkholes, and craters were observed. Technological networks and sub-services were damaged, while roads exhibited cracks. Tilting of foundations caused damage to the existing buildings.

Figure 1.

Evidence of liquefaction registered in Italy in 2012 at the site of S. Carlo: (a) external inhabited areas, extensional crack due to lateral spreading; (b) ground floor flooded by erupted sand; (c) basement floor lifted by sand under pressure; (d) car lifted up by the ejected sand; (e) building settlement.

The mechanism that takes place during liquefaction can be explained through cyclic undrained simple shear or triaxial tests carried out on saturated samples. Figure 2 shows the grain size distribution of the sand liquefied at San Carlo (S3 in Figure 2a) and its critical state line CLS in the e–p′ plane (Figure 2b). Data of two other sands, discussed in the following sections, are also shown in Figure 2.

Figure 2.

Testing sand’s (a) grain size distribution and (b) critical state lines.

Figure 3 shows undrained cyclic Tx tests on a sample of S3. The test was performed on a reconstituted sample isotropically normally consolidated at a mean effective stress p′_c = 100 kPa. At the end of consolidation, the specimen was characterized by a void ratio e = 0.79 (state parameter ψ = − 0.073 [2]) and was subjected to a stress deviator Δq = Δσ_a = ± 36 kPa (i.e. to a cyclic stress ratio CSR^TX = Δσ_a/2p′_c = 0.18).

Figure 3.

Example of a cyclic triaxial test on S3 (a) the excess pore pressure Δu as a function of the number of cycles N, (b) axial strain ε_a vs. N, (c) deviatoric stress q vs. ε_a, (d) q vs. p′.

Figure 3a and b display the variations in excess pore pressure (Δu) and axial strain (ε_a) plotted against the number of cycles (N). In Figure 3c and d, the deviatoric stress (q) is depicted in relation to both ε_a and the mean effective stress (p′). In Figure 3d, the critical state lines in compression and in extention (CSL) are also plotted, as deduced from monotonic tests.

Throughout the test, the specimen exhibits a characteristic behavior known as “cyclic mobility.” As loading progresses, positive Δu accumulates, while the mean effective stress p′ approaches zero. From the first cycle onward, the specimen undergoes an alternating response with incremental dilation (p′ increasing) and incremental contraction (p′ decreasing). The passage from incremental contraction to incremental dilation is known as phase transformation and the line intersecting transition points, known as the phase transformation line (PTL), is plotted in the q–p′ plane, as shown in Figure 3d.

The onset of liquefaction occurs when a significant accumulation of cyclic axial strain (ε_a) starts after approximately 7–8 cycles. At this stage, the pore pressure ratio Ru = Δu/p′_c approaches 1, and the effective stresses approach zero. Then the stress path reverses direction, oscillating between the critical state line in extension and compression, and the specimen exerts strain hardening, mobilizing sufficient shear strength to withstand the applied deviatoric load. However, as depicted in Figure 3b, ε_a keeps on increasing as the cyclic loading goes on.

A series of samples of S3 were reconstituted at three values of void ratio: high (e_avg = 0.78, which corresponds to ψ_avg = − 0.073), medium (e_avg = 0.73, ψ_avg = − 0.134) and low void ratio (e_avg = 0.64, ψ_avg = − 0.226).

All the specimens exhibited negative ψ after consolidation, as illustrated in Figure 4a, which shows the initial conditions of the tested samples on the ψ – p′ plane. During the cyclic Tx tests, all the samples underwent cyclic mobility.

Figure 4.

Cyclic triaxial tests on S3: (a) end of consolidation state parameter and (b) liquefaction resistance for variable state parameter.

Failure, defined as the states at which double amplitude axial strain ${\mathrm{\epsilon }}_{\mathrm{a}}^{\mathrm{DA}}$ = 5%, is depicted in Figure 4b. In this figure, the cyclic stress ratio applied in the triaxial condition (CSR^TX) is adjusted to the cyclic stress ratio for simple shear conditions (CSR^SS) in accordance with the methodology proposed in Refs. [3, 4].

$${\mathrm{CSR}}^{\mathrm{SS}}={\mathrm{CSR}}^{\mathrm{TX}}\left(1+2{\mathrm{k}}_0\right)/3\left[-\right]$$E1

where k₀ = σ′_r/σ′_a = 0.43 = stress ratio at rest, from the equation of [5].

The pore pressure ratio R_u values at failure were generally larger than 0.9. The data in Figure 4b show that the lower the state parameter, the larger the cyclic resistance, as ψ is an indicator of the direction of volumetric strains, δε_v, (dilation or contraction) during shearing; the stress ratio required to induce liquefaction at a specific number of cycles is inversely proportional to ψ. This is because, at the same stress level, the denser the soil the lower the propensity to contract and develop Δu.

Samples corresponding to a particular average ψ_avg exhibit clear correlations between cyclic stress ratio and the number of cycles. The slope of these relationships in a semi-logarithmic plot is highly influenced by the value of ψ and the interpolating curves can be interpreted through the following function:

$${\mathrm{CSR}}^{\mathrm{SS}}=\frac{\mathrm{a}{\left(1–\mathrm{\psi }\right)}^{\mathrm{b}}}{{\mathrm{N}}^{\mathrm{c}\left(1–\mathrm{\psi }\right)}}\left[-\right]$$E2

where a = 0.115, b = 3, c = 0.145, empirical constants determined by fitting experimental data. For a given number of cycles N, Eq. (2) allows to estimate the cyclic resistance ratio, CRR.

Previous studies have highlighted that the cyclic resistance of soil undergoes a significant increase even with a minor decrease in the degree of saturation [6, 7, 8, 9]. To investigate the impact of saturation on the cyclic resistance of S3, a series of cyclic tests were repeated using partially saturated samples. Achieving specific saturation levels required careful measurement of the amount of deaerated water introduced into the samples during flushing. Water circulation was halted once Sr = 80% or Sr = 90%, was attained. Throughout this process, measurements of the Skempton parameter (B) and compression wave velocity (V_P) were conducted.

Figure 5a and b illustrate the relationship between V_P, the degree of saturation Sr, and B, respectively. V_P falls within the range of 750–800 m/s when Sr = 80% (B = 0.1) and within the range of 900–1200 m/s when Sr = 90% (B = 0.3–0.7). In fully saturated S3 samples (Sr > 97% and B > 0.98), V_P is approximately equal to 1800 m/s. For V_P > 750 m/s, the data in Figure 5a can be approximated using a logarithmic function:

Figure 5.

(a) Compression wave velocity V_P vs. the degree of saturation Sr, (b) V_P vs. Skempton parameter B.

$$\mathrm{Sr}=0.17ln\left({\mathrm{V}}_{\mathrm{P}}\right)-0.29\left[-\right]$$E3

After the partial saturation process, the specimens were isotropically compressed to 100 kPa, assuming negligible influence of partial saturation on the mean effective stress, i.e. p′_c = 100 kPa. The consolidated, partially saturated specimens had e_avg = 0.7 and ψ_avg = −0.16, classifying them as medium-dense samples.

It was observed that, under undrained cyclic conditions, the cyclic stress ratio required to induce double amplitude ${\mathrm{\epsilon }}_{\mathrm{a}}^{\mathrm{DA}}$ increased as the degree of saturation decreased.

The tests revealed that, for the same state parameter and applied cyclic stress ratio, the development of axial strains and excess pore pressure was slower in unsaturated samples. The liquefaction condition was achieved for a larger number of cycles compared to saturated samples.

Unsaturated S3 at Sr = 90% and 80% has, on average, 1.2 and 2.2 times the resistance of fully saturated S3, as can be seen in Figure 6, where the CRR of unsaturated samples is normalized against the relating value at full saturation and plotted versus the degree of saturation, Sr in Figure 6a and versus the compression wave velocity V_P normalized to its full saturation value, V_P,sat, in Figure 6b.

Figure 6.

Cyclic resistance CRR of unsaturated samples normalized to the cyclic resistance of fully saturated samples plotted vs. (a) the degree of saturation Sr and (b) the normalised compression wave velocity V_P/V_P,sat.

The experimental data in Figure 6a can be fitted with an exponential function:

$$\mathrm{R}=\frac{\mathrm{CRR}\mathrm{at}\text{partial saturation}}{\mathrm{CRR}\mathrm{at}\text{full saturation}}=60{\mathrm{e}}^{-4.1\mathrm{Sr}}\left[-\right]$$E4

which, combined with Eq. (3), becomes:

$$\mathrm{R}=197{{\mathrm{V}}_{\mathrm{P}}}^{-0.7}\left[-\right]$$E5

while the data in Figure 6b can be interpreted through Eq. (6):

$$\mathrm{R}={\left({\mathrm{V}}_{\mathrm{P}}/{\mathrm{V}}_{\mathrm{P},\mathrm{sat}}\right)}^{-0.7}\left[-\right]$$E6

and the ratio R can be assumed to correct S3 liquefaction resistance accounting for the effect of partial saturation. Eq. (2) can then be re-written as:

$$\mathrm{CRR}={\left(\frac{{\mathrm{V}}_{\mathrm{P}}}{{\mathrm{V}}_{\mathrm{P},\mathrm{sat}}}\right)}^{-0.7}\frac{\mathrm{a}{\left(1–\mathrm{\psi }\right)}^{\mathrm{b}}}{{\mathrm{N}}^{\mathrm{c}\left(1–\mathrm{\psi }\right)}}\left[-\right]$$E7

3. The ISMGEO seismic centrifuge

The ISMGEO (Istituto Sperimentale Modelli Geotecnici, Italy) geotechnical centrifuge (Figure 7) is a beam centrifuge. Its symmetrical rotating arm has a diameter of 6 m and a nominal radius of 2.2 m to the model base. An external fairing rotates jointly with the arm for aerodynamic purposes. The centrifuge can reach an acceleration of 600 g bearing a payload of 400 kg [10].

Figure 7.

ISMGEO seismic centrifuge: (a) top view and (b) schemes of the centrifuge.

On each side of the symmetric arm, there is a swinging platform that accommodates the model containers, for static tests on one side and dynamic tests on the other side. One end of the arm is instrumented with a single-degree-of-freedom shaking table (as depicted in Figure 7b and 8a).

Figure 8.

(a) shaking table installed on the rotating arm and (b) equivalent shear beam box.

The swinging platform that bears the model for dynamic tests makes contact with the table in flight at 5 g and is released before further accelerating the centrifuge. The shaker is capable of replicating real input motions at the scale of the model [11].

An Equivalent Shear Beam (ESB) box [12, 13], Figure 8b, is used to simulate liquefaction. During the flight, the long side of the container is positioned vertically and aligned parallel to the centrifuge’s rotation axis. This configuration prevents any distortion effects caused by rotation on the central section of the model, where the instruments are placed and minimizes complications associated with Coriolis acceleration, being the direction of shaking also aligned parallel to the centrifuge’s rotation axis.

4. Centrifuge modeling of liquefaction

The LIQUEFACT project involved an extensive series of centrifuge tests conducted at ISMGEO. The objective of the tests was to investigate the seismic response of level ground-saturated sandy deposits, whether homogeneous or stratified, when subjected to progressively intensifying seismic excitations, ultimately leading to liquefaction. The experiments also aimed to assess the effectiveness of various liquefaction mitigation techniques. The tests provided a large database for numerical tools calibration [14, 15]. A more comprehensive description of the experimentation can be found in Refs. [16, 17].

In this section, an overview of the experimental details is provided, along with an analysis of selected results obtained from tests conducted to investigate the triggering mechanism and the efficacy of vertical and horizontal drains to be employed as mitigation techniques.

4.1 Modeling details

The tests reproduced level ground sandy deposits, around 14 m deep. The groundwater table was set at the soil surface. The geometrical scaling factor was N = 50 and the centrifugal acceleration of 50 g was imposed at the model’s base.

The testing sands are Ticino Sand, herein referred to as S1; a natural, liquefiable sand named S3, and S2, which consists of S3 after the removal of particles finer than 0.075 mm. Grains size distribution and critical state line of the test sands in the e-p′ plane are shown in Figure 2. The main physical and mechanical properties of the sands are given in Table 1.

		S1	S2	S3
γmin	(kN/m^3)	13.64	12.55	12.18
γmax	(kN/m^3)	16.67	15.75	15.77
Gs		2.68	2.69	2.69
D50	(mm)	0.53	0.17	0.15
φ′cs	(deg)	34	34.5	34.5
K*	(m/s)	2⋅10^−3	10^−4	8.4⋅10^−5

Table 1.

Soil physical and mechanical properties.

*

End of consolidation values.

The models were reconstructed at low density into the ESB container and saturated under a vacuum pressure of about −60 kPa using a viscous fluid, to avoid the discordance between the scaling ratios for time in dynamic phenomena and in diffusion phenomena [18].

Table 2 reports the average values of the void ratio and relative density of the models here discussed, referring to the pre-shaking condition. The average ψ value is also indicated. Figure 9 depicts the model schemes of the discussed tests, with all measurements provided at the prototype scale.

Test ID	Soil	Density DR* (%)	Void ratio e*	State parameter ψ*
M1_S1_GM17	S1	47	0.761	−0.13
M1_S1_GM34	S1	50	0.748	−0.14
M1_S1_GM31	S1	47.5	0.757	−0.13
M1_S2_GM17	S2	65	0.82	−0.2
M1_S3_GM17	S3	56	0.89	−0.17
M1_S1_VD1 & VD2_GM31	S1	47	0.76	−0.13
M1_S1_HD1 & HD2_GM31	S1	54.5	0.74	−0.15

Table 2.

Test program and model characteristics.

*

Average values.

Figure 9.

Model schemes: (a) homogeneous models (M1) made of S1 sand, tested applying GM17 and GM34 ground motions; (b) M1_S1 model, tested applying GM31; (c) M1 model of S2 sand, tested applying GM17; (d) M1 model made of S3 sand, tested applying GM17; (e) M1_S1 model equipped with vertical drains tested applying GM31; (f) M1_S1 model equipped with horizontal drains tested applying GM31.

The models were instrumented with miniaturized sensors (accelerometers, acc, pore pressure transducers, ppt, displacement transducers, D) placed along the mid-section.

Some models reconstituted using S1 sand were equipped with flexible silicon pipes, meant to simulate prefabricated vertical and horizontal drains. The pipe’s external and internal diameters were 6 mm and 4 mm, respectively (300 mm and 200 mm at the prototype scale,) and their hydraulic conductivity to water was 1.7 × 10⁻² m/s.

The mesh of vertical drains (VDs) was square, with a spacing (S) between drains set at 5 or 10 diameters (30 and 60 mm, equivalent to 1.5 and 3 m at the prototype scale) in models VD1 and VD2, respectively. The number of drains was 30 in VD1 and 12 in VD2 models, giving a treated area of about 45 m² and 54 m² at the prototype scale.

The mesh of horizontal drains was triangular, S was 5 or 10 diameters (30 and 60 mm, 1.5 and 3 m at the prototype scale) in models HD1 and HD2; the number of horizontal drains was 10 and 9, giving a treated area of about 9 m² and 31 m² in HD1 and HD2, respectively.

4.2 Ground motions

An overview of the principal properties of the input motions employed in the tests discussed here is provided in Table 3. Figure 10 provides examples of the time histories generated by the shaking table, and Figure 11 displays the Fourier Amplitude Spectra (FAS) corresponding to all the motions listed in Table 3. The accelerometer signals were converted into Fourier spectra using the FFT (Fast Fourier Transform) algorithm, including the tapering operation; spectra were smoothed using logarithmic smoothing with a triangular window.

Test ID	GMID	PGA (g)	d90 (s)	IA,max (m/s)
M1_S1_GM17	GM17	0.215	15.09	0.348
M1_S1_GM34	GM34	0.222	24.23	0.451
M1_S1_GM31	GM31	0.198	18.63	0.601
M1_S2_GM17	GM17	0.226	13.53	0.32
M1_S3_GM17	GM17	0.211	11.61	0.27
M1_S1_VD1 & VD2_GM31	GM31	0.187	19.83	0.573
M1_S1_HD1 & HD2_GM31	GM31	0.185	19.1	0.467

Table 3.

Input motion characteristics.

GMID = ground motion ID; PGA = peak ground acceleration; d₉₀ = duration calculated on the basis of Arias Intensity; I_A,max = maximum arias intensity.

Figure 10.

Examples of GMs times histories.

Figure 11.

Fourier amplitude spectra of the input GMs.

4.3 Same sand, input motions of increasing intensity

In this section, the results are compared of tests M1_S1_GM17/34/31 (layout in Figure 9a and b). GM17, 34, and 31 had similar PGA but increasing I_A,max. The time history of the motions and their FAS are shown in Figures 10 and 11.

Figure 12 shows the isochrones of excess pore pressure of the three models for instants ranging from 0.3 to 5 times the duration (d₉₀ in Table 3) of the input earthquake.

Figure 12.

Excess pore pressure isochrones of models (a) M1_S1_GM17, (b) M1_S1_GM34, and (c) M1_S1_GM31 at time istants ranging fro 0.3 d₉₀ to 5 d₉₀.

In interpreting the tests, the liquefaction criterion is assumed to be the point at which Δu approaches the pre-shock vertical effective stress σ′_v0 (represented in the Figure by an inclined straight line).

Based on this assumption, complete liquefaction was observed in only one model, M1_S1_GM31, at mid-depth (ppts 3) and near the ground surface (ppts 4), see Figure 12c. In this model, at all monitoring points, Δu initially increased to its maximum value within the first few seconds. Subsequently, at ppts 1 and 2, Δu began to decrease, nearly simultaneously and well before the conclusion of the ground motion. Conversely, at ppts 3 and 4, the excess pore pressure remained at its maximum level throughout the entire dynamic loading and beyond: the excess pore pressure equalized the vertical effective stress up to d₉₀ and up to 1.7d₉₀, at the depth of ppt3 and ppt 4, respectively. In a few cycles of excitation, the model exhibited a clear separation into two distinct halves: the upper part became fluidized, while the bottom part remained in a solid state [19].

The fluid state persisted until the conclusion of the dynamic loading, after which ppt3 began registering a decrease in Δu, signifying that solidification had taken place at that particular depth. It took an additional period equal to 2 times d₉₀ for the solidification front to reach ppt4. Following this, the entire soil column solidified, initiating a reconsolidation process to dissipate Δu.

The occurrence of decay from the model base upward during the earthquake suggests that the response of the sand to the seismic loading is a partially drained process. As the soil undergoes dynamic strain and starts to generate excess pore pressure, it simultaneously begins to dissipate this Δu [20, 21]. During the earthquake, Δu represents a balance between generation and dissipation.

In the lower half of the model, dissipation predominated over generation after a few loading cycles, leading to an upward fluid flow, that is, from higher Δu levels toward lower values at the surface. In the upper half (ppt3 and 4), the generation induced by the shaking and the inflow from greater depths outweighed dissipation. This led to a hydraulic gradient reaching a critical value, causing the soil to liquefy and remain in a fluidized state until the end of the seismic loading (d₉₀ at the depth of ppt3) or even longer (1.7d₉₀ at the depth of ppt4).

The above is confirmed by the evolution of surface settlement, as depicted in Figure 13. At the conclusion of the recording, the total surface settlement S_t amounted to 265 mm. Notably, 30% of this settlement occurred within the initial 2.5 seconds, a period during which all the ppts recorded a rise in Δu; 67% of the total settlement developed during the ground motion, resulting from the combination of drainage, sedimentation, and reconsolidation. Only 33% of the S_t could be attributed to post-seismic settlement, arising from solidification and reconsolidation.

Figure 13.

Ground surface settlement.

Regarding models M1_S1_GM17 and M1_S1_GM34 (depicted in Figure 12a and b), the seismic motions applied had very similar PGA compared to GM31, but they had lower values of the maximum Arias Intensity (I_A,max), as detailed in Table 3. Specifically, I_A,max was 0.348 for GM17, 0.451 for GM34, and 0.601 for GM31.

Both models exhibited a peak Δu at all depths within the first few seconds of loading but they did not liquefy; Δu began decreasing thereafter. The dominant mechanism was excess pore pressure generation only during the initial few seconds, after which dissipation and drainage became predominant.

Furthermore, the surface settlements in both models were predominantly associated with the seismic event and a great part of the final settlement in both models developed during the initial phase of excess pore pressure accumulation: 84% in model M1_S1_GM17 and 50% in model M1_S1_GM34. The settlement attributed to rapid drainage during the initial earthquake phase, characterized by the development of a high hydraulic gradient, was higher than the settlement resulting from post-earthquake reconsolidation.

4.4 Same input motion, different sands

In this section, an analysis is conducted on three tests performed on models reconstructed with various sands and subjected to the same input motions (model M1_S1_GM17, scheme in Figure 9a, results in Figure 12a; models M1_S2_GM17 and M1_S3_GM17, schemes in Figure 9c and d, results in Figure 14a and b). It should be noted that, although the same reconstitution procedure was followed, S2 and S3 exhibited greater compressibility (reflected in the slope of the critical state line, as shown in Figure 2b). During the saturation process and subsequent in-flight consolidation, these models settled more compared to S1.

Figure 14.

Excess pore pressure isochrones of models (a) M1_S2_GM17 and (b) M1_S3_GM17 at time istants ranging from 0.3 d₉₀ to 5 d₉₀.

As a result of these differences, M1_S2 and M1_S3 achieved a higher relative density than M1_S1 at the conclusion of the Ng consolidation phase. Specifically, the relative density of M1_S2 was 18% higher, and that of M1_S3 was 9% higher, as indicated in Table 2.

As model S1 described above, neither model S2 nor S3 experienced liquefaction (Δu < σ′_v0).

Model S2 (Figure 14a) developed very low pore pressure at every depth except at ppt5. At this depth, once the max Δu was attained, dissipation was triggered at a very low rate. The soil beneath ppt5, instead, continued to accumulate pore pressure slightly, even after the dynamic phase had concluded.

The lower values of Δu observed in S2 compared to S1 can be attributed to the lower initial ψ, as indicated in Table 2. This initial condition led to a reduced tendency for the soil to contract during dynamic loading and, consequently, resulted in minor soil surface settlement, as illustrated in Figure 13. In total, 90% of the ground surface settlement was associated with the seismic event, confirming the occurrence of partial drainage during seismic excitation.

In the case of the natural S3 sand model, which experienced a less intense earthquake, the values of excess pore pressure were similar to those observed in model 2 except near the ground surface (ppt5), where the Δu values were lower and dissipation onset during the earthquake but at a very slow rate, while the soil below continued to accumulate pore pressure slightly, even after the dynamic phase had concluded. As a result, the S3 deposit exhibited a lower overall ground settlement, as shown in Figure 13.

Figures 12a, 14a, and b reveal that S1 exhibited a considerably higher dissipation rate compared to both S2 and S3 sands. After the seismic shock had concluded, dissipation in S1 was at 50% at all depths. In contrast, S2 and S3 were still accumulating below the depth of ppt5, which was close to a permeable boundary represented by the surface.

When the data recording was interrupted, the lower half of models 2 and 3 had not yet initiated the dissipation of Δu. This observation aligns with the permeability values of the three sands at the test density, as presented in Table 1. Specifically, S1 exhibited a permeability that was an order of magnitude higher than that of S2 and S3. S3, despite having a fine content of 12%, being slightly looser during the test compared to S2 had similar permeability.

From the Δu measures from tests in Figures 12 and 14, the pore pressure ratio Ru = Δu/σ′_v0 has been computed.

In Figure 15 the max Ru registered at the base ppt (ppt1 in all models), and at the shallower ppts, are plotted as a function of the I_A,max of the input earthquake (Table 3). The results of three additional tests (M1_S1 models equipped with a shallow foundation F) non-discussed herein are also reported (measures from the free field part of the models).

Figure 15.

Ru of shallower and deeper ppts vs. max arias intensity.

The distribution of Ru for both the deepest and the shallowest ppt seems unique and follows a sigmoid function, as indicated by the interpolation curves outlined in the Figure. For I_A,max equal to 0.6, the shallower part of the model reaches the liquefaction condition, meanwhile, at the bottom of the container Ru does not exceed 0.6.

The end of recordings settlements are plotted in Figure 16 vs the max Arias Intensity and describe an S-shape function.

Figure 16.

Superficial settlement vs. max arias intensity.

4.5 Liquefaction mitigation using vertical and horizontal drains

Figure 17 shows the isochrones of excess pore pressure Δu measured in the models equipped with vertical and horizontal drains (test layout in Figure 9e and f, FAS of the applied input options in Figure 11).

Figure 17.

Excess pore pressure isochrones of models (a) M1_S1_VD1_GM31, (b) M1_S1_VD2_GM31, (c) M1_S1_HD1_GM31, and (d) M1_S1_HD2_GM31 at time istants ranging from 0.3 d₉₀ to 5 d₉₀.

As highlighted above, the shallower half of the untreated model M1_S1_GM31 underwent complete liquefaction, while liquefaction was prevented in the treated areas of both configurations of vertical drains.

As expected, the larger the spacing the lower the Δu reduction. The extent of VDs efficacy can be appreciated in terms of Ru in Figure 15, where ppts 1, 3 and 5 are shown and compared with the sigmoid functions obtained from untreated models.

The lower Ru observed at ppt3 and ppt5 is attributed to the influence of the drains. This conclusion is supported by the measurements at ppt1, which are consistent with the Ru function depicted in Figure 15: ppt1 was positioned at the model base, a significant distance (20 diameters) away from the treated area where the drains were installed.

In both VD1 and VD2 configurations, the maximum Δu was reached within a few seconds, followed by Δu decay. This decay process began well before the conclusion of the ground motion. Consequently, by the end of the seismic loading, more than 50% of the measured Δu had dispersed, in stark contrast to the untreated model.

The superficial settlements were exclusively co-seismic. In terms of the balance between Δu generation and dissipation, vertical drains made the latter phenomenon prevalent after just a few cycles of loading. VDs reduced the rate of Δu accumulation, the maximum Δu, and accelerated the rate of Δu decay.

D2 (D1 not working) recorded a final settlement of 244 mm. The difference with the 265 mm measured in M1_S1_GM31, which is slightly less than −10%, is considered to be within the range of experimental variation and is consistent with the settlement expected for the applied intensity based on the Settlement-S-function depicted in Figure 16.

The primary effect of drains is indeed to prevent Δu accumulation. However, if no densification is induced by their installation they do not impact the sand’s tendency to contract and generate excess pore pressure when subjected to seismic vibrations.

Thus, in free field conditions, the effectiveness of drains in limiting settlements is negligible, as drains primarily facilitate the dissipation of excess pore pressure (Δu) without preventing the soil from undergoing strain. However, in the presence of buildings, drains are also effective in mitigating settlements. By preventing soil liquefaction, they help prevent the sinking of buildings [22, 23, 24], provided that they are adequately distributed beneath the entire footprint of the structure [25].

As to the models treated with HDs, account has to be taken for the applied input motion weaker than in models with VDs and in the untreated model. As discussed above, tests on untreated models have shown that I_A,max influences both the max Ru and the superficial settlements.

To assess the effect of horizontal drains (HDs) despite the weaker input motion, the Ru functions of Figure 15 can be used as a reference. For I_A,max = 0.467, at the depth of ppt1 and ppts 3/5 Ru values of approximately 0.3 and larger than 0.9 were expectable. The measured values were 0.26 and 0.5.

Ru = 0.26 is comparable to the expected value and is considered within the range of experimental variation. Ru = 0.5 implies that HDs induced a reduction in Ru of about −45%, confirming their effectiveness.

Regarding surface settlements, the models with HDs developed settlements equal to 95 mm and 130 mm in HD1 and HD2, respectively. These values are consistent with what was expected for I_A,max = 0.467 based on untreated models (as shown in Figure 16). This confirms that HDs are effective in preventing the accumulation of excess pore pressure and liquefaction but do not significantly influence Δu generation and subsequent settlement due to reconsolidation.

As with vertical drains (VDs), in the presence of HDs, the maximum Δu was reached within a few seconds, and Δu decay began prior the end of the dynamic loading, albeit a little slower than in the models equipped with vertical drains.

5. Liquefaction assessment using a simplified procedure

Liquefaction resistance is typically determined through the interpretation of field test results using empirical correlations based on the results of standard penetration tests [26] or, more effectively and widely, cone penetration tests [27, 28, 29, 30, 31].

Both the cone penetration resistance q_c and the undrained cyclic resistance CRR of an uncemented and unaged soil depend on its nature, stress level, and density, that is, on its state parameter ψ, which allows to anticipate the direction of volumetric strains, δε_v, (dilation or contraction) during shearing.

Therefore, ψ can be used to relate CRR to the tip resistance of CPTs.

In this section, the results of centrifuge CPT tests and cyclic undrained triaxial tests carried out on S1 sand are used to link CRR to the cone resistance.

5.1 Cyclic resistance of S1

Undrained cyclic Tx tests on S1 were performed on reconstituted samples, isotropically consolidated at p′_c = 100 kPa.

The consolidated void ratio of the specimens was: e_avg = 0.742, which corresponds to ψ_avg = − 0.132, e_avg = 0.676, ψ_avg = − 0.201 and e_avg = 0.582, ψ_avg = − 0.295.

Figure 18 illustrates the applied cyclic stress ratio CSR as a function of the number of cycles at which the double amplitude axial strain ${\mathrm{\epsilon }}_{\mathrm{a}}^{\mathrm{DA}}$ reached 5% (assumed as failure criterion). The Ru values at failure varied from 0.8 to 0.97. Only the samples TS4_14_01 and TS4_14_03 developed ${\mathrm{\epsilon }}_{\mathrm{a}}^{\mathrm{DA}}$ less than 5%, so their failure condition was considered as the attainment of Ru = 0.95, at N = 60 and N = 900, respectively.

Figure 18.

Cyclic strength of S1.

In Figure 18, CSR^TX is converted into CSR^SSviaEq. (1); for S1, k₀ = 0.44 and CSR^SS = 0.63·CSR^TX. Figure 18 shows also the interpolation curves of the experimental data using Eq. (2), which, calibrated for S1, returns the following calibration coeffcients: a = 0.071, b = 7.8, c = 0.177.

5.2 Calibration chamber centrifuge CPTs in S1

A series of 27 centrifuge CPTs were carried out on dry S1 models, reconstituted at low (e ≈ 0.63), medium (e ≈ 0.73), and high void ratio (e ≈ 0.82). The tests were run at three levels of centrifugal acceleration: 30 g – 50 g – 80 g, corresponding to increasing values of effective stresses.

Calibration Chamber (CC) CPT tests carried out by [32, 33] show comparable penetration resistance for dry and saturated conditions.

A miniaturized electrical piezocone (diameter d_c = 11.3 mm, apex angle of 60°, sleeve friction 11.3 mm in diameter, and 37 mm long) was employed for the tests (Figure 19).

Figure 19.

Model scheme and model picture with a view of the ISMGEO miniaturized piezocone before penetration.

Figure 20 shows the measured tip resistance q_c represented vs. the mean effective stress p′:

Figure 20.

CPTs in S1: tip resistance q_c as a function of mean effective stress p′ (computed adopting k₀ = 0.44).

$${\mathrm{p}}^{\prime }={\mathrm{\sigma }}_{\mathrm{v}}^{\prime }\left(1+2{\mathrm{k}}_0\right)/3\left[\mathrm{kPa}\right]$$E8

σ′_v = vertical effective stress, computed accounting for the acceleration field distortion;

k₀ = σ′_h/σ′_v = stress ratio at rest.

The values of ψ at rest at the beginning and stop of penetration are highlighted in Figure 20. The soil models are rather homogeneous and the tests are repeatable. A less-than-linear rise of q_c with p′, for a given void ratio, can be observed.

As shown by [34, 35, 36, 37], during the penetration the void ratio can decrease (contraction) or, more likely, increase (dilation) due to shearing. Dilation causes a stress increment around the tip, with respect to the stress value at rest, proportional to the dilatancy. Dilatancy can be related to the state parameter ψ at rest, before penetration [37].

Consequently, the two major contributions of q_c can be considered:

• the overburden stresses at the depth of the tip, herein represented by the mean effective stress p′;

• the increment of stresses around the tip caused by the penetration, due to dilatancy, representable by ψ.

In a homogeneous centrifuge model, ψ at rest increases with depth, so the dilative tendency of soil decreases.

So, what happens around the tip during the penetration is an increase in stress level which causes q_c to rise, and a contemporary decrease in the soil dilative tendency, so a rate of increase in q_c reduction.

To separate these two major effects, the following procedures were followed. The effect of stresses at rest at constant ψ has been defined as:

$${\mathrm{q}}_{\mathrm{c}}/{\mathrm{p}}_{\mathrm{a}}=\mathit{f}{\left({\mathrm{p}}^{\prime }/{\mathrm{p}}_{\mathrm{a}}\right)}^{\mathrm{\beta }}\left[-\right]$$E9

The exponential function represents the impact of the overburden stresses at the depth of the cone tip on q_c. The best fit of the experimental constant ψ cone resistance profiles gave β = 0.8, a slightly lower value than 1, as recommended by [38]. In this paper, the measured cone resistance was normalized using this β value as follows:

$$\left({\mathrm{q}}_{\mathrm{c}}/{\mathrm{p}}_{\mathrm{a}}\right)·{\left({\mathrm{p}}_{\mathrm{a}}/{\mathrm{p}}^{\prime }\right)}^{0.8}={{\mathrm{q}}_{\mathrm{c}}}^{\ast }\left[-\right]$$E10

and q_c* is plotted as a function of p′ in Figure 21.

Figure 21.

Normalized tip resistance q* as a function of mean effective stress p′.

Removing the influence of p′ on q_c through the normalization, the effect of ψ on q_c can be appreciated in Figure 21: while crossing less dilative soil (rising of ψ with depth), the reduction of dilatancy around the tip causes a non-linearly q_c* reduction. In a constant ψ soil profile, q_c* would have been nearly constant with depth.

The effect of ψ on q_c* is highlighted by plotting q_c* vs. ψ, as shown in Figure 22. The trend in Figure 22 can be interpreted with the following equation [38]:

Figure 22.

Normalized q* vs. the state parameter ψ.

$${{\mathrm{q}}_{\mathrm{c}}}^{\ast }=\mathrm{k}·{\mathrm{e}}^{–\mathrm{m\psi }}\left[-\right]$$E11

where m and k are dimensionless fitting parameters, equal to:

m = 8.1, k = 28.3 (12,720 data, R² = 0.96).

It’s worth noting that k is the normalized q_c* value when ψ = 0, while m indicates the influence of dilatancy on the cone resistance at a given of ψ.

The normalized cone resistance can be re-written as:

$$\left({\mathrm{q}}_{\mathrm{c}}/{\mathrm{p}}_{\mathrm{a}}\right)·={\left({\mathrm{p}}^{\prime }/{\mathrm{p}}_{\mathrm{a}}\right)}^{\mathrm{\beta }}\mathrm{k}·{\mathrm{e}}^{–\mathrm{m\psi }}\left[-\right]$$E12

which represents the dependency of q_c on the overburden stresses acting at the depth of the cone and on the soil dilative tendency at that depth.

5.3 CRR from CPT trough ψ: a simplified method

Ψ can be assumed as an independent variable governing both the cyclic stress resistance and the normalized cone resistance of the tested soils, and the ciclic resistance ratio, CRR can be derived by combining Eqs. (2) and (12):

$$\mathrm{CRR}=\frac{\mathrm{a}{\left[1+\frac{1}{\mathrm{m}}ln\left(\frac{{\mathrm{q}}_{\mathrm{c}}^{\ast }}{\mathrm{k}}\right)\right]}^{\mathrm{b}}}{{\mathrm{N}}^{\mathrm{c}\left[1+\frac{1}{\mathrm{m}}ln\left(\frac{{\mathrm{q}}_{\mathrm{c}}^{\ast }}{\mathrm{k}}\right)\right]}}\left[-\right]$$E13

where a, b, c, m, and k are the fitting parameters of Eqs. (2) and (12). These parameters have been calibrated for two sands other than S1: S3, described above (see [1]) and Toyoura sand (TOS [39]). These parameters are:

• Toyoura sand: a = 0.037, b = 10.7, c = 0.247, m = 9.8, k = 23.9

• S3: a = 0.115, b = 3, c = 0.145, m = 7.42, k = 27.44

• S1: a = 0.071, b = 7.8, c = 0.177, m = 8.1, k = 28.3

Values of m and k of Eq. (12) are also available for other two silica sands tested in centrifuge [40], Fontainebleau NE34 (FNE34S. 90% quartz, 8% feldspar, 2% mica) and Stava Sand (SS, 55–60% quartz, 20–25% feldspar, 16% fluorite, 7–8% calcite). The latter is a natural sand containing 30% of non-plastic silt, obtained from the mine tailings dam in Stava (Italy), which collapsed the 19th of July 1985 [41]. The fitting parameters are:

• FNE34S: m = 9.4, k = 17

• SS: m = 5.3, k = 57

The q_c*-ψ curves of all the tested sands are shown in Figure 23. A total of 64 CPTs are displayed. For each sand, a trendline is represented in Figure 23. An average relationship for all the sands is drawn in yellow. Its equation is:

Figure 23.

Normalized q_c* vs. ψ for five sands. The yellow curve is the average relationship.

$${{\mathrm{q}}_{\mathrm{c}}}^{\ast }=30.6·{\mathrm{e}}^{–\text{8ψ}}\left[-\right]$$E14

Eq. (13) and the set of average fitting parameters given in Table 4 can be used to derive CRR from q_c measured in situ. The parameter β of Eq. (10) is 0.8 in clean sand and can be assumed equal to 1 in the presence of silty sand.

A	b	c	m	k
0.074	7.17	0.19	−8	30.6

Table 4.

Average fitting parameters of Eq. (13).

Eq. (13) turns into Eq. (15) to account for the degree of saturation:

$$\mathrm{CRR}={\left(\frac{{\mathrm{V}}_{\mathrm{P}}}{{\mathrm{V}}_{\mathrm{P},\mathrm{sat}}}\right)}^{-0.7}\frac{\mathrm{a}{\left[1+\frac{1}{\mathrm{m}}ln\left(\frac{{\mathrm{q}}_{\mathrm{c}}^{\ast }}{\mathrm{k}}\right)\right]}^{\mathrm{b}}}{{\mathrm{N}}^{\mathrm{c}\left[1+\frac{1}{\mathrm{m}}ln\left(\frac{{\mathrm{q}}_{\mathrm{c}}^{\ast }}{\mathrm{k}}\right)\right]}}\left[-\right]$$E15

6. Conclusions

This chapter deals with earthquake-induced liquefaction and its modeling, from the element volume scale to small-scale physical models.

The first part of the chapter describes a case study in Italy presenting some experimental observations gained through cyclic undrained simple shear and triaxial tests carried out on saturated samples of a natural soil that experienced liquefaction during the 2012 Emilia seismic sequence.

The cyclic resistance ratio CRR of a granular soil has been expressed as a function of the applied number of cycles N and on the state parameter ψ, through Eq. (2).

To evaluate the increase in liquefaction resistance due to the partial saturation reference can be made to the velocity of compression waves through Eq. (7).

The second part of the chapter describes the centrifuge modeling of liquefaction; some details of the Italian seismic apparatuses are given. Some results of a large campaign of physical model tests are presented. The tests have shown that in free field conditions, the distribution of max Ru vs. the max Arias intensity follows a sigmoid function.

In the particular test layout discussed here, the function was found to be independent of the specific testing sand. This suggests that it can be adopted to predict Ru for earthquakes with varying maximum Arias Intensity. Also, the superficial free field settlements appear to follow a unique S-shaped settlement function.

Vertical and horizontal drains, applied as a liquefaction mitigation technique, do not affect Δu generation; instead, their action is to prevent the accumulation of excess pore pressure, thereby avoiding liquefaction. However, they do not mitigate the volumetric strains induced by dynamic loadings. As a result, the settlements observed in models equipped with both horizontal drains and vertical drains are comparable to those measured in untreated models for the same earthquake intensity. Therefore, in free field conditions, drains have limited effectiveness in reducing settlement.

The third part of the chapter describes the results of combined centrifuge and laboratory experiments conducted to derive the cyclic resistance of soil from the cone penetration resistance.

The cone penetration resistance may be linked to the state parameter viaEq. (12).

Eqs. (12) and (2), combined, turn into Eq. (13), which can be used to assess the liquefaction resistance of sandy soil from in situ CPTs, and which, with respect to the empirical correlations typically used in the geotechnical practice, has the advantage of a direct experimental calibration.

Acknowledgments

The LIQUEFACT project received funding from the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement No. 700748. This support is gratefully acknowledged by the authors.

All the experimental data available at https://www.zenodo.org/record/1281598#.W-mWVOhKjIU

The authors acknowledge the Istituto Sperimentale Modelli Geotecnici (ISMGEO) for carrying out the experimentation and making the results available for the elaborations and speculations described in this paper.