Adjustment of a numerical model for pore pressure generation during an earthquake

Summary

This article proposes methodology for evaluating the accuracy of the pore pressure generation model devised by Byrne, as implemented in a commercial software program using a Mohr-Coulomb-type failure criterion and a Finn constitutive model. The different empirical formulas of liquefaction developed by Seed and Idriss are reviewed, as well as various constitutive models specified in the literature, emphasizing the selection of the Finn model for the liquefaction study. In the analysis a comparison is carried out using the factors of safety against liquefaction (FSLs) devised by Seed and Idriss and the adapted formula by Boulanger and Idriss. The analysis assumes a hypothesis to verify whether a soil element is liquefied. The results are then compared with those of a numerical model that simulates a soil column, the base of which is subjected to the same seismic inputs of varying magnitudes, M_w, and peak ground accelerations, Pga, to which the empirical model was subjected. Adjusted equations are provided on the based on that comparison to allow for the calibration of the Byrne equation using the (N₁)₆₀ value obtained via a standard penetration test (SPT), for the study of liquefaction problems in situations in which there are earthquakes of varying magnitudes.

Keywords:

Engineering and technology – Plant resistance to abiotic stress – Valleys – Seeds – Deformation – Shear stresses – Statistical distributions – Chemical composition

Introduction

In studies of soil dynamics, situations are presented that totally account for neither the adjusted response of a numerical calculation model with respect to the data input from field tests or site observations, or the relationship between the model and the field data. To this end, we aim to facilitate such studies by improving in the aforementioned data input scheme (in particular, with the (N₁)₆₀ value obtained from the SPT), correcting the values provided by a field study and soil characterization tests, to ensure a better fit to the formula stored in a computer program and provide a calculated solution with greater accuracy.

On this basis, new adjusted equations are provided to enable data input calibration in a behavior model within a numerical calculation model using the Byrne equation [1], by using the (N₁)₆₀ value obtained via the SPT, in situations in which earthquakes of varying magnitude occur.

Two summaries of previous studies, regarding the numerical models of soil behavior and the various empirical models devised by Seed and Idriss since 1971 are provided in detail below. Subsequently, a calculation adjustment is presented with the methodology proposed by the authors to provide the equations that will control the input data used to study of the behavior of granular soil, which is susceptible to liquefaction, against dynamic actions.

Current state of knowledge

Liquefaction is defined as the loss of the shear resistance of soil subjected to monotonic or cyclical loads due to the tendency of less competent soils to reorganize their structure under shear stresses. The term was coined by Mogami and Kubo [2].

In many geotechnical projects, the model of Seed and Idriss is used to determine whether the soil can undergo liquefaction. Generally, if the results of the Seed and Idriss model indicate the possibility of liquefaction, some actions are carried out (soil improvements) to prevent it.

Geotechnical numerical models allow the total stresses to be calculated to determine if liquefaction occurs. The difficulty of using numerical models is in deciding how pore pressures are generated during the earthquake. To analyze, liquefaction using a numerical model, an adjustment to the equation of Byrne [1] is proposed to guarantee an improvement in modeling results.

If a simple case is analyzed (if a soil is liquefiable) by considering a specific earthquake, a certain discrepancy between the results obtained with the model of Seed and Idriss and those with a numerical model will arise. The purpose of this article is to propose a correction of the Byrne model that aligns the results of these two methods. This correction will allow more complex situations to be evaluated through numerical modeling.

Empirical formula for the calculation of FSL

The state-of-practice simplified methods for the evaluation of liquefaction potential were developed using historical cases and field measurements at sites characterized with on-site tests (such as SPT). These simplified methods are generally expressed as a deterministic and semiempirical models.

In this approach, Seed and Idriss [3] concluded that the liquefaction of a soil would occur if FSL (assuming gently sloping ground with free-field conditions) is ≤1. FSL is defined as the relationship between the cyclic resistance ratio (CRR) and the cyclic stress ratio (CSR); therefore, if FSL>1, then the liquefaction of the soil would be unlikely.

Due to the use of field data collection during real earthquakes, researchers have progressively supplemented the initial Seed and Idriss formula [3] through, adaptations such as those shown in Table 1:

Analyzing the previously published data and observing the recent research detailed in the previous table, the advances developed by Boulanger and Idriss [17] were selected to be used in the empirical formulation for the calculation of liquefaction potential. The authors find that the details of these updates to the CSR formulas are the most appropriate.

Constitutive models of soil subjected to dynamic loads

The study of soil dynamics has been addressed in recent years by a large number of researchers. These works provide constitutive models that integrate both the characteristics of the soils (damping and resistance to cyclic loads) and the properties of the boundary conditions (site).

Over time, numerical models, based on the global constitutive models, have been developed in the study of soil dynamics problems and are incorporated in specialized calculation programs. Some of these models are detailed in Table 2:

Over the past few years, the Research Group of the Department of Engineering and Ground Morphology of the Advanced School of Civil Engineering (ETSICCP, as per its Spanish acronym) at the Technical University of Madrid (UPM, as per its Spanish acronym), has made advances in numerical modeling, especially with the Finn model (with Byrne equation), which is the model most adjusted to suit the work of the Research Group.

The development of the study and calibration of the Finn model is reflected in the analysis of Galindo [27] on the laboratory cyclic shear tests carried out by Patiño [28], conducted under different combinations of stresses representing dynamic loads at different frequencies, and then applied to a real case study by Soriano [30]: the sinking of some concrete caissons in the mouth of the Port of Barcelona due to the liquefaction of the foundation sands.

The Finn model is defined by an equation (of volumetric deformation) that depends on the (N₁)₆₀ value from a SPT test, with allows for the calculation of pore pressure generation, unlike the empirical method that utilizes the formulation of Boulanger and Idriss to determine the occurrence of liquefaction. Therefore, the only parameter that could be varied was the N of the SPT.

Therefore, the Finn model, implemented in the explicit finite difference program FLAC3D, is chosen to carry out the study presented in this article.

Finn Model

Martin et al. [23] described the pore pressure generation mechanism of this model, highlighting that the relationship between the irrecoverable volume deformations and the range of the cyclic shear strains depends on the confining stress.

Based on the work of Martin et al. [23], Byrne [1] established a more simple formula for the Finn model that depends solely on one value of the soil, the (N₁)₆₀ value, which takes into account the usual parameters (angle of internal friction, deformation modulus, etc.). These two formulas are intrinsically related; in fact, if the exponential term of Byrne’s formula was written with reference to the first two terms of the Taylor series expansion, Byrne’s model would be the same as Martin’s formula.

This model, implemented in FLAC3D, the explicit finite difference program, is chosen to carry out the study presented in this article.

Proposed methodology

An improvement in the numerical modeling used to perform dynamic calculations in liquefiable soils is proposed by adjusting the Byrne equation. In particular, the input of the (N₁)₆₀ value of the SPT test is modified here because it is the fundamental parameter of the Byrne model.

To this end, the empirical method of global reference and the Finn model with the formulation of Byrne are used, with the data provided by nine seismic inputs with different M_w and Pga values, to establish the generation of pore pressures within the FLAC3D environment and the FSL of Seed and Idriss (adapted by Boulanger and Idriss).

The following development plan, based on equalizing the FSLs to obtain the adjusted Byrne formula [1], is presented:

• Selection of a group of representative seismic inputs to be applied to the base of the model; each group is characterized by its moment magnitude (M_w) and by its peak ground acceleration (Pga).

• Determination of the (N₁)_60cs value for each soil element with FSL = 1, according to the formulation of the liquefaction model by Seed and Idriss [3], adapted by Boulanger and Idriss [17]. These factors are all dependent on the features of the input acceleration and the depth of the soil element considered.

• Formulation of the numerical model and assignment of the soil properties considered using the (N₁)_60cs value obtained in the second section.

• Solution of the numerical model for each seismic input and the subsequent obtainment of the defined FSL as the coefficient between the effective mean stress before the earthquake σm′(t0) and the excess pore pressure immediately after: Δpp = pp(t_f)-pp(t₀).

• If the achieved results are not satisfactory, i.e., if the obtained FSL, according to the previously described definition, is not equal to 1, an iterative calculation is carried out to determine to the (N₁)_60cs value, which guarantees FSL = 1.

Selection of input accelerograms

Nine seismic inputs were selected with a moment magnitude of M_w ≥ 6.5, given that lower magnitudes and the associated peak ground accelerations commonly require (N₁)_60cs values that are extremely low and cause liquefaction in soils. Three moment magnitudes are considered, 6.5, 7.5, and 8.5, each with low, medium and high peak ground accelerations. In Table 3, the records used and their main features are shown.

All the records were originally filtered for low frequencies, generally lower than 0.1 Hz. In addition, their baselines have been corrected. Fig 1 shows one of the nine FLAC3D input accelerograms used, that of Puerto Quellon (PQ).

Application of the seed and Idriss model (1971), adapted by Boulanger and Idriss (2014), and the attainment of (N₁)_60cs

In liquefaction risk assessments, it is customary to use the Seed and Idriss model (1971) or the subsequent updates, as these provide an FSL that describes the possibility of a particular granular soil (characterized by (N₁)_60cs) liquefying during a seismic stimulation with a certain moment magnitude, M_w, and peak ground acceleration, Pga, when FSL<1. The denomination of the (N₁)_60cs value stems from the denomination of Seed and Idriss of the corrected (N₁)₆₀ value of the SPT and with the added “cs” detail, referring to the adjustment for “clean sand”. Therefore, the same name is used in this article.

The Seed and Idriss model is based on the concept of resistance to cyclic loading, using the CRR, the maximum tangential stress of a soil sample before reaching the yield strength in a cyclic simple shear test is, confined under a determined effective vertical stress and a cyclic strain is applied with a given frequency.

It is generally accepted that earthquakes are multifrequency signals. Nevertheless, Seed and Idriss find a certain relationship between the application frequency of cyclic loading and the moment magnitude. Therefore, the CRR parameter can be defined as a function of the features of the sample, its confinement and the predominant frequencies.

With cyclic resistance defined, cyclic loading must be defined such that the relationship between them provides the FSL. In this sense, Seed and Idriss define cyclic loading, or CSR, as the maximum shear stress that a soil element undergoes during an earthquake due to the maximum acceleration. Given that the maximum acceleration is not directly representative of the average effect of the earthquake, Seed and Idriss [3], with the formula adapted by Boulanger and Idriss [17], adjusted this value with a coefficient of 0.65. Furthermore, this model is based on one soil column that behaves as a rigid body; therefore, a parameter that takes into account the deformability of the soil should also be included. The expressions of the CRR and the CSR are shown below.

Formula provided by Boulanger and Idriss in 2014:

The resistance obtained thereby is valid for earthquakes with magnitudes of approximately 7.5 and an effective vertical stress of 1 atmosphere. Therefore, to introduce the characteristics of the earthquake and the depth of the soil element in the calculations, the following modifications are made in the CRR formula:

• Correction for depth, K_σ:

  

  where:

  

• Correction based on the seismic moment, MSF:

  

  where:

  

The formula to apply both corrections is as follows:

Once the resistance of the soil element is calculated, the shear stress exerted on the soil element must be obtained using the next formula [3]:

r_d is the coefficient of the shear stress reduction, the arguments of which appear inside the parentheses (in radians).

The following coefficients, α and β, are the adjustment values and are dependent on the depth, z, below the surface of the soil (in meters):

In conclusion, to obtain the FSL:

By means of the previous expressions, it is possible to determine the (N₁)_60cs values to verify FSL = 1 for each of the earthquakes listed in Table 3 and for different depths (from 0 to 20 m). For this purpose, we use a spreadsheet and a tool (SOLVER) to minimize the sum of residuals, defined as (FSL-1)², by changing the cells of (N₁)₆₀. In Fig 2, an example is shown for the Puerto Quellon (PQ) case.

The results shown in the following Table 4 are obtained by carrying out the calculation of the (N₁)_60cs values for all the seismic inputs and all the depths considered.

Numerical model in FLAC3D: Characteristics and properties

In this stage, the numerical model is defined, and its initial and boundary conditions, constitutive model and input properties are assigned.

The selected model represents a 20 m column of granular soil saturated subjected to a seismic input at its base. For this purpose, a mesh, made of 20 elements with dimensions of 1 x 1 x 1 m³ is constructed in FLAC3D, as shown in Fig 3.

Periodic conditions linking the nodes on one side of the model to those on the opposite side are used as boundary conditions to induce identical behavior. The lateral confinement that cannot be implemented through the anchorage of the lateral displacements in this dynamic model is thus represented.

The initial stress conditions are due to the existing overburden material load (lithostatic); hence, an at-rest thrust coefficient of k₀ = 1, which is considered to adequately represent the stress state of incompetent grounds, is used. In geological time periods, this coefficient tend to dissipate the deviator stresses and therefore equalize the vertical and horizontal stresses (Heim’s rule).

The Finn constitutive model incorporates the Byrne equation [1] for pore pressure generation due to cyclic strain. This constitutive model, proposed by Martin et al. [23], is a classic model with a Mohr-Coulomb failure criterion that, via an internal algorithm, tests the load-unload-reload cycles for each element in the model and determines the cyclic strain generated by each cycle. Then, the algorithm calculates the increase in volumetric strain associated with the cycle and the shear stress using the following expression:

Lastly, the model generates an increase in pore pressure for each cycle considered; this increase is obtained by multiplying the volume increase associated with the cycle, provided by the equation of the volumetric strain and by the soil bulk modulus.

Given that the Finn model includes a classic Mohr-Coulomb failure criterion, two elastic and three plastic properties are required. For this study, pore pressure generation is analyzed in more detail. Thus, the plastic properties are not very relevant because the model is developed by assuming an elastic case. Therefore, a high value is assigned to avoid breakage of the soil and destabilization in the calculation of the model.

With regard to the elastic properties and given that the relevant variable is the (N₁)_60cs value, the initial shear modulus (G_i) is calculated, in Pa, using the following elastic relationship:

Therefore, drawing on both formulas, the initial shear modulus is:

The properties assigned to the Finn model incorporated in FLAC3D are listed in, Tables 5, 6 and 7.

The ff_latency parameter is employed in FLAC3D to delineate the number of steps in the calculation. This internal parameter remains defined as the minimum number of time steps between the stress reversals, [33]. ff_latency is a control parameter used to stabilize the model and help with convergence, serving as a filter so that there are not too many oscillations during modeling.

The previously identified properties refer to the static model; thus, in addition to these properties, the dynamic behavior of the soil must be considered. This behavior is modeled by adopting one of the hysteretic models implemented in FLAC3D, specifically, the sigmoidal model Sig3. The equation of the aforementioned model is displayed below, and the parameters listed in Table 8.

The reduction curve of the shear modulus generated by these parameters, shown in Fig 4, is selected after various tests, based on those used by Idriss and Boulanger [16], to obtain the reduction factor of the shear stress by depth (r_d), to find which value provides the best results.

Once the properties of the mechanical model are defined, the configuration of the flow model settings must be described. In this case, a model with hydromechanical coupling and impeded flow is used; in other words, the generated volume undergoes variations that the ground cannot dissipate. This hypothesis is justified given that in the brief time lapse of the duration of the seismic stimulation, the produced flow could be considered negligible even if the permeability of the soil is of an average type. Shown in Table 9.

Results

This section presents the results obtained in the verification process of the equation of the pore pressure generation, Byrne [1]. The stimulated numerical model is run with each of the seismic inputs displayed in Table 3 applied to its base, and the results are expressed in terms of the FSL for each soil element, defined as

It is well known that this definition of the FSLs diverges from that recommended by Seed and Idriss in their proposal. However, it is widely held that these two expressions must be equivalent at FSL = 1. Therefore, to verify the usefulness of Byrne’s equation to generate excess pore pressure until reaching the liquefaction of the ground, this assumption is reasonable.

If the Byrne model were to yield the expected results, the obtained FSL should be equal to unity. However, it can be expected that, given that the verification equation is not an analytical expression with boundary conditions and specific initial inputs or with a totally defined dynamic load, differences arise between the empirical and numerical results. Hence, analytical studies such as this work are needed. In this last case, considering the results provided by the Seed and Idriss model (1971), which was adapted by Boulanger and Idriss in 2014, we will proceed to modify the parameters of the Byrne equation to achieve good compatibility between these two formulas for FSL = 1. In addition, a convenient modification for FSL = 1.3 will be researched.

To analyze the suitability of the FSL vs. depth results, an analysis of the data by using basic statistical variables (average, standard deviation, and coefficient of variation) is performed.

Initial FSL (precorrection)

To resolve each of the cases, the process begins by statistically balancing the model and storing the effective mean stress and the pore pressure of each element in one memory position. Then, the corresponding (N₁)_60cs value is assigned according to the element depth (see example Table 4) and the features of the earthquake, namely, M_w and Pga (see Table 3). Finally, the seismic input is applied to the base of the soil column from t = t₀ to t = t_f.

During the development of the earthquake, the pore pressure excess is tested, measured as the difference between the pressure at time t and the previously stored pressure at time t₀; with this detail and the effective mean pressure, also measured and stored, the evolution of the FSL is continuously monitored.

Displayed (Fig 5) below is the output of the model in terms of the evolution of the pore pressure for three depths: 4.5, 9.5 and 14.5 m. This graph supplements the accelerogram of the prescribed earthquake.

In Fig 5, it is possible to confirm when the increase in the pore pressure undergoes its greatest change with the peak ground acceleration; however, after this, the majority of excess pore pressure generated during strong shaking is (generally) low, except when there are acceleration peaks similar to the Pga in question at a later stage.

In Fig 6, the results of the FSL are provided for each range of M_w and each depth. Finally, it is possible to verify the degree of closeness between the Byrne model (incorporated in FLAC3D) and that of Seed and Idriss, whose formula was adapted by Boulanger and Idriss.

It is possible to observe that the obtained result is acceptable, even if cases appear such as that of Imperial Valley (IV) and Coquimbo (CP), which provide an FSL that is sensibly lower than unity. Similarly, the opposite case is shown for Maule (M), which presents a notably higher FSL. For an overview of the obtained results, these cases are grouped together in the histogram of Fig 7.

The distribution of the FSL displayed in the previous histogram clearly reflects a high degree of dispersion and a bias in underestimating the risk of liquefaction with respect to the Seed and Idriss formula; the statistics of the results are displayed below in Table 10.

The table highlights the percentage of the FSL results greater than 1.3, given that FSL>1.3 is generally required for it to be accepted as a reliable prediction of the risk of liquefaction.

Correction of the Byrne equation (1991)

In view of the above results, the authors consider it necessary to revise the parameters of the Byrne equation to better adjust the Seed and Idriss method adapted by Boulanger and Idriss, which is postulated (by consensus) to be correct.

In the first instance, the error with respect to FSL = 1 vs. (N₁)_60cs is plotted to verify some type of relationship, as shown in Fig 8.

The previous figure clearly suggests the possibility of acting on the same (N₁)_60cs value so that the correction would be scarce or null for the high and low values and high for the average values.

To carry out this modification of (N₁)_60cs, a statistical analysis of the data obtained via linear approximation is carried out to improve the results, validating the proposal.

Proceed as follows: the original value of (N₁)_60cs, provided by the model of Seed and Idriss adapted by Boulanger and Idriss, will be tentatively corrected with a linear adjustment, depending on the error, so that a new test value, i.e., (N₁)_{60cs_corr} will be obtained. This correction will be carried out for each seismic input of the nine selected cases.

Once the first attempt is carried out, the nine cases are recalculated in FLAC3D, and their approximation to FSL = 1 is verified; the value provided for the previous equation is a_i = -0.5. Clearly, some cases will be corrected to a greater degree than others; thus, in the following attempt, only those that did not converge to FSL = 1 are corrected.

The results of this approximation process, i.e., (N₁)_{60cs_corr,} which was introduced in the Byrne equation, verify the Seed and Idriss model, as shown in Table 11, and the result of the adjustment in terms of the factor of safety is shown in Fig 9.

With regard to Fig 9, the result is deemed adequate, and the (N₁)_{60cs_corr} value is compared with that of (N₁)_60cs, as shown Fig 10; however, the data for depths under 12 m will not be taken into account to achieve a better adjustment because, at greater depths, there is increased dispersion in the r_d factor of the Seed and Idriss model adapted by Boulanger and Idriss. Beyond that depth, r_d is significantly affected by the frequency content of the seismic signal as Ishihara [34] and Golesorkhi [8] mentioned and collected by Idriss and Boulanger [16].

In the previous figure, the exponential relationship between (N₁)_60cs and (N₁)_{60cs_corr} can be clearly seen. By exploiting this trend, it is possible to determine the values for the earthquake magnitude, M_w, by establishing three correlations, for M_w values of 6.5, 7.5 and 8.5.

Therefore four relationships are proposed, as listed below, to allow a shift from (N₁)_60cs to (N₁)_{60cs_corr}:

This allows us to obtain the (N₁)_{60cs_corr} value and serves as a new entry in the Byrne equation. In this way, the specific coefficients of each moment magnitude range can be used.

Final FSL (postcorrection)

The next step is to verify the results with the proposed modification for the (N₁)_60cs value. The first step consists of correcting values of Table 4 according to the previous equation, with the adequate coefficients, and assigning them to the elements of the numerical model. Then, the FLAC3D model is developed with the nine seismic hypotheses, and the results are plotted for each earthquake in terms of the FSL and depth at the end of the dynamic action in Fig 11.

Visual inspection of the previous figure suggests that the correction carried out via the proposed correlation for four magnitudes M_w improves the result the Byrne equation produces in the numerical model. However, for the precorrection hypotheses, the histogram of values is presented in Fig 12.

The factor of safety distribution is tightly concentrated around FSL = 1, although a percentage of cases persist that undervalue the resistance to liquefaction. Table 12 shows the statistical values of the previous distribution.

Considering the numerical parameters that describe the distribution of the results, the improvement in the degree of dispersion is evident, as it decreases from 28.8% to 13.6%. Furthermore, only 5% of the cases present a factor of safety in excess of 1.3; consequently, overestimation of the resistance to liquefaction is limited to a few cases.

In addition to the verification of the case of FSL = 1, the calculations are repeated for FSL = 1.3 to verify that the results provided by the proposed correction are achievable for other ranges of the factor of safety. Therefore, first, in the spreadsheet prepared for this purpose, we obtain the (N₁)_60cs value that verifies a factor of safety of 1.3 for all earthquakes and depths, which are then introduced into the FLAC3D numerical model for the nine tested cases.

The results of this last hypothesis are shown below in Figs. 13 and 14 and Table 13:

It is important to remember that the definition of the FSL proposed in this report and that used by Seed and Idriss [3] adapted by Boulanger and Idriss [17] are not the same and that they are only comparable in the case of FSL = 1. This discrepancy could cause the results obtained in this last hypothesis to diverge slightly from FSL = 1.3; nevertheless, as the distribution in Fig 14 and Table 13 show, only 10% of the cases are under 1.

Conclusions

This study sought to verify the results of the Byrne equation [1], implemented in the Finn model by using FLAC3D, with respect to the generation of pore pressure in soils subjected to cyclic strain. The output results of the (N₁)₆₀ value of the SPT are compared with the FSL provided by the method proposed by Seed and Idriss [3] and adapted by Boulanger and Idriss [17]. The results are then introduced in the numerical model, with a purpose-built geometry, to which diverse seismic inputs are applied to its base. The FSL is defined as the ratio between effective mean stress and the pore pressure increase. Therefore, this factor of safety is obtained using the numerical model output values.

The factors of safety obtained by the numerical model using the selected (N₁)_60cs value, which, according to the Seed and Idriss method adapted by Boulanger and Idriss, should be unitary, do approach this value. Nevertheless, the results display a certain dispersion and many cases present greater factors of safety (FSL>1.3 ≈ 50%), underestimating the risk of liquefaction.

A correction of the (N₁)_60cs value is proposed according to equation (22) and its coefficients (23) based on M_w, which allows us to obtain the (N₁)_{60cs_corr} value and serves as a new entry in the Byrne equation.

After repeating the modeling by using the parameter of the corrected Byrne equation, a clear improvement is obtained in the output; the distribution of the factors of safety has a significantly reduced dispersion and fewer cases overestimate the risk of liquefaction, FSL>1.3 ≈ 5%.

Additionally, the calculation is repeated with the numerical model for a factor of safety of 1.3, verifying that the correction can obtain consistent FSLs, with a distribution whose average is 1.29 and that it provides only 10% false positives, namely, soil elements for which the numerical model determines them to be liquefied but the analysis does not.

Therefore, an equation is provided, serving as a correction of the Byrne equation for (N₁)₆₀ values of the normalized SPT test between 5 and 35 and soil depths to approximately 12 m in areas affected by earthquakes with magnitudes between 6.5 and 8.5.

Supporting information

S1 File [pdf]
Supporting information.