A CPT-based effective stress analysis procedure for liquefaction assessment

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Contents lists available at ScienceDirect

Soil Dynamics and Earthquake Engineering journal homepage: http://www.elsevier.com/locate/soildyn

A CPT-based effective stress analysis procedure for liquefaction assessment Nikolaos Ntritsos *, Misko Cubrinovski University of Canterbury, Civil and Natural Resources Engineering, Private Bag, 4800, Christchurch, New Zealand

A R T I C L E I N F O

A B S T R A C T

Keywords: Cone penetration test Earthquake Effective stress analysis Liquefaction System response

There is a growing awareness of the need for the earthquake engineering practice to incorporate in addition to empirical approaches in evaluation of liquefaction hazards advanced methods which can more realistically represent soil behaviour during earthquakes. Currently, this implementation is hindered by a number of chal­ lenges mainly associated with the amount of data and user-experience required for such advanced methods. In this paper, we present key steps of an advanced seismic effective-stress analysis procedure, which on the one hand can be fully automated and, on the other hand, requires no additional input (at least for preliminary ap­ plications) compared to simplified cone penetration test (CPT)-based liquefaction procedures. In this way, effective-stress analysis can be routinely applied for quick, yet more robust estimations of liquefaction hazards, in a similar fashion to the simplified procedures. Important insights regarding the dynamic interactions in liquefying soils and the actual system response of a deposit can be gained from such analyses, as illustrated with the application to two sites from Christchurch, New Zealand.

1. Introduction Liquefaction assessment including triggering and consequences of liquefaction is routinely carried out by geotechnical engineers using semi-empirical methods which are largely based on observations from case histories (e.g. Refs. [3,20,28,39], among others). Since their first introduction by Seed and Idriss [29], these simplified liquefaction assessment procedures (i.e. simplified procedures in the following) have undergone a number of updates taking advantage of the continuously increasing data from earthquake events, and have gained wide accep­ tance in engineering practice. Despite their widespread use, application of the simplified procedures to numerous case histories from recent earthquakes suggests that simplified procedures may be less effective in the evaluation of liquefiable soils other than uniform clean sands. Such “problematic” soils and ground conditions may involve silts, silty sands with non-plastic, low-plasticity or plastic fines, gravel-sand-silt mix­ tures, and deposits with complex stratigraphy composed of interbedded liquefiable and non-liquefiable soils [11]. In parallel with the development and evolution of simplified pro­ cedures, liquefaction research has provided invaluable insights into the mechanics of soil liquefaction and consequent ground deformation, hence, setting the basis for the development of the seismic effectivestress analysis, which is one of the most advanced numerical methods used in geotechnical engineering today. It allows to simulate many

important aspects of the complex dynamic behaviour of soils during earthquakes including pore pressure development, reduction in the effective stress and resulting deformations, and their effects on foun­ dations and supported structures. Importantly, such analysis considers the response of the deposit as a whole allowing for interactions between layers in the dynamic response (e.g. liquefaction effects on ground motion) and through pore water pressure redistribution and water flow (e.g. seepage-induced liquefaction). On the downside, effective-stress analysis can be challenging in the application due to the specialized data and skills required, commonly including high-quality sampling and laboratory testing of soils, selection of appropriate input ground mo­ tions, complex calibration procedures for sophisticated constitutive models, and also in-depth understanding by the user of the phenomena considered, constitutive relationships used and numerical solutions adopted in the analysis. These drawbacks have limited the application of the effective-stress analysis to critical lifelines and structures where the cost and effort to perform the analysis can be generally justified. In modern seismic design philosophies (i.e. Performance-Based Earthquake Engineering), reasonable estimation of liquefactioninduced damage to land and structures is a key requirement even for less critical projects. Therefore, it is important that engineers have the tools to consider key response mechanisms, and predict ground dis­ placements with reasonable accuracy, but without the need for overly onerous numerical analyses (at least at a preliminary analysis stage).

* Corresponding author. E-mail addresses: [email protected] (N. Ntritsos), [email protected] (M. Cubrinovski). https://doi.org/10.1016/j.soildyn.2020.106063 Received 26 October 2019; Received in revised form 9 January 2020; Accepted 21 January 2020 Available online 29 January 2020 0267-7261/© 2020 Elsevier Ltd. All rights reserved.

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

This can be achieved either by appropriately modifying the simplified methods to account for key aspects in the soil response which are currently either ignored or inadequately addressed in these methods (e. g. material and state-characterization of soils other than clean sands, and system-response effects) or by standardizing advanced effectivestress analysis procedures, allowing for easy calibration and mini­ mizing the required input and modelling decisions from the user. The present paper focuses on the second approach by introducing a step-bystep automated procedure for one-dimensional (1D) seismic effectivestress analysis of free-field level ground based solely on conventional input data from CPTs, similarly to the simplified CPT-based liquefaction assessment procedures. In the following, a brief overview of relevant background of simplified procedures is first given. Then, the main steps of the proposed procedure including determination of a simplified soil profile, estimation of constitutive model parameters, definition of input ground motion and other aspects of numerical modelling and analysis are described in detail. Finally, example applications for two casehistory sites from Christchurch, New Zealand are presented.

The earthquake-induced CSR (for either a case-history or forward site liquefaction assessment) can be estimated at any depth z in the deposit either using site response analysis or alternatively, approxi­ mated as: CSR ¼ 0:65

σ vo amax rd ðMw ; zÞ σ ’vo g

(1)

where σ vo and σ ’vo are the initial total and effective vertical stress at depth z, amax is the maximum ground surface acceleration in g, and rd is the shear stress reduction coefficient which can be expressed as a function of Mw and z [20]. In this relationship, it is assumed that amax and CSR are not affected by excess pore pressures during earthquake shaking. 2.2. Consequences of liquefaction without considering interactions between layers The case histories used in the development of liquefaction triggering correlations essentially reflect the overall response of soil deposits and associated severity of liquefaction manifestation for a specific earth­ quake excitation. However, despite the intent to capture the overall performance of the deposit at a given site, in the simplified procedures each layer is considered in isolation, and a factor of safety against liquefaction triggering (FS), maximum shear (γmax ) and volumetric strains (εv ) are estimated separately, and independently, for each layer. In other words, interactions between different layers in the deposit through the dynamic response, excess pore water pressures, and water flow are ignored. Hence, principal mechanisms of interaction or systemresponse effects of liquefying deposits that potentially contribute to the severity of liquefaction manifestation are not accounted for in the simplified procedures. Liquefaction damage indices, such as LSN [35] and LPI [21] use specific weight functions to quantify the damage po­ tential of liquefying layers depending on their proximity to the ground surface, but still, when calculating the cumulative damage indices a simple superposition of independent effects from each layer is used, while cross-interactions between layers during the development of liquefaction and post-liquefaction triggering are simply ignored. Cubrinovski et al. [14] have identified this limitation as one of the key reasons for mispredictions of the simplified procedures at Christchurch case-history sites.

2. Simplified liquefaction evaluation procedures 2.1. Triggering of liquefaction In the stress-based approach, the most widely used approach for assessing liquefaction triggering, earthquake-induced cyclic shear stresses (CSR) are compared with the cyclic (liquefaction) resistance (CRR) of the soil. Liquefaction is predicted at those depths in the deposit where the induced shear stresses exceed the cyclic resistance or the � factor of safety against liquefaction triggering is less than 1.0 FS ¼ � CRR CSR � 1:0 . Empirical liquefaction triggering curves, such as the one shown in Fig. 1, correlating CRR with some in-situ test index are typically used to estimate the soil resistance against liquefaction. In the deterministic procedure, these curves are in essence boundary lines that separate case histories in which liquefaction was manifested at the ground surface from case histories in which such manifestation was not evident. The representation of all case-history data onto a single graph of CSR versus in-situ penetration resistance is achieved by adjusting the data points for each case-history to a common reference condition of earthquake magnitude Mw ¼ 7.5 and effective overburden stress σ ’vo ¼ 1 atm. The measured in-situ penetration resistance is also adjusted to account for the effects of overburden stress and fines content (FC) on the penetration resistance. The corrected CPT resistance, termed the equivalent clean sand cone penetration resistance qc1Ncs , is the one used in the CPT liquefaction triggering correlation of Fig. 1.

3. Proposed effective-stress analysis procedure 3.1. Introduction The proposed effective-stress analysis procedure involves the following key steps: (1) Determination of a simplified soil profile from CPT data, (2) Determination of characteristic soil behaviour and associated constitutive model parameters for liquefiable soil layers, (3) Determination of constitutive model parameters for nonliquefiable soil layers, (4) Definition of input ground motion(s), and (5) Definition of numerical model and analysis parameters. In the first step, the nearly continuous CPT profile is discretized into a number of distinct layers or depth intervals over which the CPT data can be approximated by constant values. This soil profile “simplification” assists in the identification of characteristic layers in the deposit, and allows for more rigorous engineering scrutiny and interpretation. Steps 2 and 3 involve the definition of a target soil behaviour and the subse­ quent determination of constitutive model parameters that can simulate this target behaviour. For liquefiable soils, the main objective is to accurately simulate the cyclic (liquefaction) resistance for a range of soil densities, confining stresses and number of loading cycles of interest.

Fig. 1. CPT-based liquefaction triggering correlation of Boulanger and Idriss [3], for Mw ¼ 7:5 and σ ’vo ¼ 100 kPa. 2

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

The target liquefaction resistance is determined using the simplified procedure for liquefaction triggering based on empirical CPT charts (e.g. Fig. 1). For non-liquefiable soils, the target cyclic stress-strain relation­ ship is defined using strain-dependent stiffness degradation and damp­ ing ratio curves, commonly employed in site response analyses. In the fourth step, an automated ground motion selection algorithm is used to select ground motions representative of the seismic hazard and earth­ quake scenarios of interest. In the final step, the numerical model is defined by selecting appropriate model dimensions, mesh (element) size, boundary conditions and initial stress state of the soil. Also, anal­ ysis parameters such as computational time increment, integration scheme and numerical damping are adopted, and the dynamic effectivestress analysis is then executed. Fig. 2 summarizes the key steps in the proposed CPT-based effectivestress analysis procedure. In the following sections, important aspects for the main elements are discussed in detail, with particular emphasis placed on the description of steps 1 and 2, as they are essential to the application of the proposed procedure for liquefaction assessment. A key feature of the proposed procedure is that it has been automated in a way that it requires the same input data and uses the same definition of liquefaction resistance as the simplified procedures. These features of the proposed effective-stress analysis procedure significantly facilitate its application into practice, and also provide a basis for rigorous com­ parisons of the outcomes of effective-stress analyses and simplified procedures. A set of Matlab programs were developed to process and interpret the input CPT data, prepare the required input for the effective-stress analysis code, call the code to perform the numerical simulations, and post-process the analysis output. In the numerical simulations, the finite element code DIANA-J [30] and the Stress-Density constitutive model [12,13] were employed. It is important to note however that the pro­ posed procedure is generally applicable to any finite element or finite difference code and liquefaction-oriented constitutive model provided that these can capture key aspects of the cyclic soil behaviour in a fully coupled (soil-water) dynamic analysis, as discussed later in this paper.

empirical correlations between soil (behaviour) types and the various CPT measurements. These features of the CPT allow its direct use in automated analysis programs for liquefaction evaluations which perform point-by-point calculations (typically at 2 cm intervals) without the need for explicit consideration of soil sample data and site stratig­ raphy. Appealing though it may seem, in reality, such an approach prevents a rigorous soil profile scrutiny and interpretation in various aspects of the liquefaction assessment, including characterization of interface zones of soil layers with highly contrasting penetration re­ sistances [2], and also consideration of different levels of uncertainties in the estimates of the cyclic resistance associated with varying quality of material and behavioural characterization of different soil types [11]. The importance of this latter point can be elucidated if one considers, for example, the fact that in the simplified procedures (e.g. Refs. [3,28]), a significant adjustment of the measured penetration resistance (and consequently, cyclic resistance) is prescribed for soils other than clean sands, based solely on the highly uncertain fines content of these soils. Careful inspection of the CPT results along with soil sample data and thereafter definition of a simplified soil profile with distinct soil layers allow to address these issues (at least), and also assist in the identifica­ tion and characterization of critical layers in the deposit. This process can substantially enhance the engineering interpretation of the overall site characteristics and facilitate the identification of relevant systemresponse mechanisms. Furthermore, the determination of a simplified soil profile is essential for the numerical model where the problem domain, the soil deposit in this case, needs to be discretized into an assemblage of finite elements with specific model properties and constitutive behaviour. To demonstrate the equivalence of the proposed effective-stress analysis procedure with the CPT-based simplified procedures, in this paper we define simplified soil profiles solely based on CPT measure­ ments. Needless to say, however, when additional geologic information, in-situ and laboratory test data are available, one could make best use of all available information in determining a simplified soil profile for more robust liquefaction analysis. Simplified soil profiles can be determined from CPT data, by iden­ tifying depth intervals over which the corrected penetration resistance (qc1Ncs ) and the soil behaviour type index (Ic ) can be approximated by constant values, as illustrated in the two leftmost plots of Fig. 3 for an actual site from Christchurch. In this analysis, the actual qc1Ncs values (irregular traces in Fig. 3b) were calculated using the equivalent clean sand adjustments for fines content specified in Ref. [3] with the fines content inferred from the correlation between Ic and FC using a fitting parameter CFC ¼ 0. Key requirements in the layering definition are that the discretization is fine enough to allow detection of thin seems of liquefiable soils and that the Ic and qc1Ncs profiles of the simplified soil profile (solid straight lines in Fig. 3a and 3b) are as close as possible to the actual CPT traces. On this basis, distinct soil layers with represen­ tative qc1Ncs and Ic values can be defined. Even though the profile simplification inevitably results in some deviations from the actual Ic and qc1Ncs traces, the effects of these differences on the outcomes of the liquefaction assessment are negligible, as illustrated in Fig. 3c, 3d and 3e where the simplified profile yields practically identical factors of safety (FS) and liquefaction damage indices (LPI and LSN) with the analyses using the actual Ic and qc1Ncs profiles.

3.2. Determination of simplified soil profile from CPT data (step 1) 3.2.1. General considerations Two principal advantages of CPT over other in-situ tests are that it provides a practically continuous record of the penetration resistance throughout depth and that it allows to infer the soil profile from

3.2.2. CPT profile discretization algorithm Determining a simplified profile from the original CPT traces can often be quite tedious, particularly when a large number of CPT profiles has to be processed. Several approaches for automating this process have been proposed in the literature (e.g. Refs. [9,36]). For the purposes of this study, a more practical yet efficient algorithm for automating the CPT profile simplification was developed, driven by the key re­ quirements outlined in the previous paragraph. A flow chart describing the main steps and details of this algorithm is presented in Fig. 4. The main goal of the algorithm is to identify a simplified profile with

Fig. 2. Key steps in the proposed CPT-based effective-stress analysis procedure. 3

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Fig. 3. Determination of simplified soil profile from CPT data: (a) Ic ; (b) qc1Ncs ; and calculation of: (c) Factors of Safety (FoS), and damage indices: (d) LPI, and (e) LSN, throughout the top 10 m of the deposit; irregular lines refer to the original CPT traces whereas solid straight lines indicate the adopted values for the simplified profile.

Ic and qc1Ncs step-functions of depth (i.e. piecewise constant), denoted as Ic and qc1Ncs respectively, that minimize the dispersion of the actual Ic and qc1Ncs data around them. Apart from the actual Ic and qc1Ncs traces of the CPT, the algorithm requires four input parameters: the maximum tolerable coefficients of variation for Ic and qc1Ncs within each layer, denoted as CVIc and CVqc1Ncs respectively, and the minimum (Tmin ) and maximum layer thicknesses (Tmax ) for the simplified profile. Input parameter values adopted for the analyses presented in this paper are indicated in the flow chart of Fig. 4. After inspection of the output of the algorithm (i.e. generated simplified profile), these parameters can be adjusted accordingly if a finer or coarser layering is desired for the analysis. To assess its effectiveness in determining simplified soil profiles for actual soil deposits, the algorithm described in Fig. 4 was applied to 50 sites from Christchurch for which high-quality CPT data were available. The soil profiles at these sites varied from relatively uniform clean sand deposits to highly stratified deposits consisting of interbedded layers of sand-silt mixtures, silts and clays. After a positive visual inspection of the obtained simplified profiles for each site, triggering analyses and eval­ uation of liquefaction damage indices (LSN and LPI) were performed using both the simplified (Ic and qc1Ncs ) and the actual (Ic and qc1Ncs ) profiles, in a similar fashion to the verification analysis performed for the example site shown in Fig. 3. In these analyses, two earthquake scenarios (i.e. two different combinations of peak ground surface ac­ celeration amax and earthquake magnitude Mw ) were considered for each site. Results from all analyses for the simplified and actual profiles are comparatively shown in Fig. 5, in terms of LPI (Fig. 5a), and LSN (Fig. 5b). The vast majority of the computed damage indices in these plots fall on or close to the 1:1 line indicating that the estimated lique­ faction performance from the analyses using the simplified profiles is

consistent with the results of the analyses with the original profiles. This validates the use of the proposed algorithm for determination of a representative simplified CPT soil profile for liquefaction analysis. That being said, it is important to emphasize that the proposed algorithm is intended to assist users and not replace their engineering judgement in determining a simplified profile. As alluded to previously, prediction of fictitious layers at layer interfaces is not uncommon in these automated procedures, and so the user must always review the output and make adjustments when deemed necessary. The use of automated procedures for correcting CPT data for thin-layer and transition zone effects (e.g. Ref. [2]) prior to the application of the algorithm can significantly help to avoid such issues. 3.3. Modelling of liquefiable soils (step 2) 3.3.1. Stress-Density Model An elastic-plastic constitutive model, called the Stress-Density Model (SDM), is employed in the effective-stress analysis [12,13]. The model utilizes the state-concept approach for modelling the combined effects of density and confining stress on stress-strain behaviour of sand. Conse­ quently, it can simulate the behavior of a given sand at any density and confining stress by using the same set of material parameters. Key as­ sumptions in the elastic-plastic formulation of SDM are: (i) continuous yielding or vanishing elastic region, (ii) combined isotropic and kine­ matic hardening plasticity, (iii) dependence of the plastic strain incre­ ment direction on the stress increment direction (hypoplasticity), (iv) modified hyperbolic stress-strain relationship, and (v) an energy based stress-dilatancy relationship. In terms of soil behaviour, this translates into a capability of the model to accurately simulate highly nonlinear stress-strain behaviour both under monotonic loading (from small 4

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Fig. 4. Flow chart describing the main steps and details of the algorithm for determination of simplified soil profiles from CPT data.

5

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

conditions. Each sample is tested at a different cyclic stress ratio in order to establish the LRC across various CSRs or over a range of approxi­ mately Nc ¼ 1 30 cycles, which is the most relevant number of cycles for earthquake loading. In the absence of experimental data, LRCs can be derived by directly following conventional procedures for liquefaction evaluation based on empirical liquefaction triggering charts. In the present study, the Boulanger and Idriss [3] liquefaction triggering pro­ cedure was used to determine a set of representative LRCs through the following steps: (1) For a given qc1Ncs value, the cyclic resistance ratio CRR corre­ sponding to earthquake magnitude Mw ¼ 7:5 and effective over­ burden stress of σ ’vo ¼ 100 kPa is estimated using Eq. (2): � � qc1Ncs �qc1Ncs �2 �qc1Ncs �3 �qc1Ncs �4 CRRM¼7:5; σ’vo ¼1 ¼ exp þ 2:80 þ 113 1000 140 137 (2)

(2) The maximum value of the magnitude scaling factor (MSFmax ) is estimated using Eq. (3): �q �3 c1Ncs MSFmax ¼ 1:09 þ � 2:2 (3) 180

(3) The b value describing the slope of the LRC in the CSR is computed from MSFmax using Eq. (4): b ¼ c0 þ c1 ðMSFmax Þ þ ​ c2 ðMSFmax Þ2 þ c3 ðMSFmax Þ3 þ c4 ðMSFmax Þ4

Nc space (4)

3:0176, c1 ¼ 7:0217, c2 ¼ 5:7685, c3 ¼ 2:152 and where c0 ¼ c4 ¼ 0:3. Eq. (4) provides an approximate expression of the MSFmax b rela­ tionship presented in Figure A16 of [3]. (4) The number of equivalent cycles corresponding to earthquake magnitude Mw ¼ 7:5 (NM¼7:5 ) is then estimated using Eq. (5) and Eq. (6):

Fig. 5. Comparison between damage index predictions using the actual and simplified profiles for 50 Christchurch sites considering two different earth­ quake excitations: (a) LPI, and (b) LSN values.

(5)

1

NM¼7:5 ¼ Nmin ⋅ðMSFmax Þb where

strains to large strains or steady state of deformation) and irregular cyclic loading. The model was specifically tailored for analysis of liquefaction problems and has been extensively verified through rigorous simulations of down-hole array records at liquefaction sites, seismic centrifuge tests, large-scale shake table tests and numerous case histories [15]. SDM has been implemented in the finite element code DIANA-J [30] and is currently at the final stage of its verification in FLAC2D and OpenSees.

� Nmin ¼

1:0 0:65

�1b � � 3 4

(6)

(5) Next, for a given effective overburden stress σ’vo , the overburden stress correction factor Kσ is estimated using Eqs. (7) and (8): � ’ � σ (7) Kσ ¼ 1 Cσ ln vo � 1:1 Pa

3.3.2. Liquefaction resistance For liquefaction problems, the key requirement from the constitutive model is to accurately simulate the development of excess pore water pressures under irregular cyclic loading (earthquake excitation). This ability needs to be demonstrated through a simulation of target lique­ faction resistance curves in element test simulations. Liquefaction resistance curves (LRCs) represent the combination of shear stress amplitude (CSR) and number of cycles (Nc ) required to cause liquefaction or a certain level of strain in the soil, e.g. 3%, 5% or 7.5% double amplitude strain. They are typically expressed as CSR Nc re­ lationships and are used as a key soil property in the calibration of constitutive models for effective-stress analysis. LRCs are commonly derived from a series of liquefaction tests on soil samples in the labo­ ratory in which samples of “identical” (or similar) density are subjected to cyclic shear stresses of uniform amplitude under undrained

where Pa is the atmospheric pressure, and Cσ ¼

1 37:3

8:27ðqc1Ncs Þ0:264

� 0:3

(8)

(6) For a given number of equivalent loading cycles Nc , the magni­ tude scaling factor MSF is estimated using Eq. (9): � �b NM¼7:5 MSF ¼ (9) Nc

6

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

(7) Finally, the cyclic resistance for the given combination of qc1Ncs and σ’vo at Nc cycles is estimated using Eq. (10): CRR ¼ CRRM¼7:5;σ’vo ¼1 ⋅MSF⋅Kσ

Table 1 Stress-Density Model parameters. Material parameter

(10)

Repeating steps 5 and 6 for different Nc values, over the range of cycles of interest, yields the LRC for the soil with penetration resistance qc1Ncs and effective overburden stress σ ’vo . In the same way, empirical LRCs can be derived for any given qc1Ncs σ’vo combination which is within the range of applicability of the corresponding simplified procedures. The above procedure effectively reduces simplified liquefaction triggering procedures to a set of LRCs. Such LRCs derived from the Boulanger and Idriss [3] liquefaction triggering procedure following the above steps 1–7, for various qc1Ncs values and effective overburden stress of σ’vo ¼ 100 kPa, are shown with solid lines in Fig. 6.

Elastic parameters Shear constant A Poisson’s ratio ν Exponent n Reference lines UR-line (Void ratios and normal stresses in kPa) (eU, pU) QSS-line (Void ratios and normal stresses in kPa) (eQ, pQ)

3.3.3. SDM calibration for liquefiable soils The goal is now to determine SDM parameters that can reproduce the target LRCs (defined as above) in element test simulations. SDM has four group of input parameters: elastic parameters, reference states, stressstrain parameters and dilatancy parameters. Cubrinovski and Ishihara [13] have established a set of SDM parameter values for dry-pluviated Toyoura sand, after a series of laboratory tests, including drained and undrained, monotonic and cyclic (liquefaction) tests. As Toyoura sand is often used as a representative clean sand for liquefaction studies, these values were used as a basis for the determination of a set of new SDM parameters that can simulate the target (empirical) liquefaction resis­ tance curves. In this process, the stress-strain parameters were kept the same as for Toyoura sand, whereas other parameters of the model were varied in a trial-and-error process to identify the best-fit values providing the most accurate simulation for the whole set of the target LRCs. The new SDM parameter values resulted from the above calibra­ tion process are listed in Table 1. For reference, the original SDM pa­ rameters for dry-pluviated Toyoura sand are also shown in this table. It is important to note that changes in elastic, dilatancy and reference line parameters were relatively small and definitely within an acceptable range of variation. Fig. 6 shows with open symbols the LRCs obtained using simulations with the new set of SDM parameters. These SDM simulated LRCs are established through a series of element test simulations in which a soil element, at an initial stress-density state, is subjected to a given

Stress-strain parameters Peak stress ratio coefficients a1, b1 Max. shear modulus coefficients a2, b2 Min. shear modulus coefficients a3, b3 Degradation constant f Dilatancy parameters Dilatancy coef. (small strains) μ0 Dilatancy coef. (cyclic loading) μcyc Critical state stress ratio M Dilatancy strain Sc

Toyoura sand [13]

Generic sand (compatible with [3])

250 0.20 0.60

310 0.25 0.80

(0.895, � 400) (0.877, 1)

(0.895, � 400)

(0.877, 10) (0.873, 30) (0.870, 50) (0.860, 100) (0.850, 200) (0.833, 400)

(0.874, 10) (0.873, 30) (0.872, 50) (0.871, 100) (0.868, 200) (0.860, 400)

0.592, 0.021 291, 55

0.592, 0.021 291, 55

98, 13 4

98, 13 4

0.22 0.00 0.607 0.0055

0.22 -0.02 0.620 0.0040

(0.875, 1)

amplitude of uniform stress cycles. Liquefaction triggering is assumed to occur in the simulation cycle at which the double amplitude (DA) shear strain exceeds 5%. The number of cycles required to cause liquefaction and develop 5% DA strain is then used together with the applied level of cyclic stress to define the soil resistance (i.e. one symbol in Fig. 6) at the given stress-density state. The simulation of the target LRCs shown in Fig. 6 is considered sufficiently accurate for liquefaction analysis. Small discrepancies do occur at higher stress ratios, but overall, simulated and target curves show a reasonably good agreement across Nc ¼ 1 50 cycles. In addi­ tion, the adopted SDM parameter values provide a reasonably accurate simulation of the tail of the target LRCs, at low cyclic stress amplitudes, which defines the threshold CSR separating between liquefaction and no-liquefaction. This detail is particularly important, as underestimation of the liquefaction resistance for low shear stresses can lead to sub­ stantial overprediction of pore water pressures in the seismic effective stress analysis. It is important to note that, a single set of values for SDM parameters was used to simulate the target LRCs across different qc1Ncs values (i.e. soil densities), by only adjusting the initial void ratio e in SDM to achieve a good fit. The use of a single set of parameter values across various densities is a distinctive feature of the SDM model which comes as a result of the incorporation of the state-concept characterization of soil behavior into the model. To facilitate the application of the SDM model over a range of different qc1Ncs values that may result from the profile simplification process, an expression for estimating the void ratio e from qc1Ncs was developed as follows: e¼

0:315½1 þ expð

0:128qc1Ncs þ 18:8Þ�

0:142

þ 0:931

(11)

For a target LRC associated with a specific qc1Ncs value, this expres­ sion was used to derive the respective e value for the SDM in the element test simulations shown in Fig. 6. Using the minimum and maximum void ratios for Toyoura sand, emin ¼ 0:616 and emax ¼ 0:988 respectively, Eq. (11) can be rewritten in terms of relative density Dr as:

Fig. 6. Target liquefaction resistance curves (solid lines) obtained using the Boulanger and Idriss [3] liquefaction triggering procedure and SDM simulated LRCs (open symbols) for different qc1Ncs values (i.e. different void ratios in SDM) and σ ’vo ¼ 100 kPa.

Dr ¼ 0:847½1 þ expð

0:128qc1Ncs þ 18:8Þ�

0:142

þ 0:153

(12)

Fig. 7 shows the relationship of Eq. (12) together with three common 7

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

empirical Dr qc1Ncs expressions for clean sands [20,27,34]. It can be seen that Eq. (12) slightly underestimates Dr for qc1Ncs between 60 and 100, compared to the empirical estimates, whereas it overestimates Dr for qc1Ncs greater than 130. Nevertheless, given the large uncertainty in these empirical expressions, expected differences between laboratory-based and field-based LRCs, and the success of Eqs. (11) and (12) in simulating the target LRCs across a wide range of qc1Ncs values, while using the same set of values for SDM parameters, which was the primary goal in this calibration exercise, these small deviations of Eq. (12) from the range of empirical Dr qc1Ncs correlations are considered acceptable. Another benefit from the embodiment of the state-concept modelling to SDM, is that it accounts for the combined effects of density and confining stress on the liquefaction resistance and cyclic stress-strain behaviour of sand. In the simplified methods, the effect of confining stress on liquefaction resistance is addressed through the use of the Kσ factor (Eqs. (7) and (8)), defined as: Kσ ¼

CRRσ’vo CRRσ’vo ¼1

(13)

Fig. 8. Comparison between Boulanger and Idriss [3] (solid lines) and SDM simulated (dashed lines with open symbols) Kσ σ ’vo relationships for various qc1Ncs values.

where CRRσ’vo is the CRR of a soil under a specific value of σ ’vo , and CRRσ’vo ¼1 is the CRR of the same soil for a reference value of σ’vo ¼ 1 atm. Fig. 8 shows (with solid lines) a set of Kσ σ’vo curves obtained using the Boulanger and Idriss [3] expressions (Eqs. (7) and (8)), and the corresponding SDM simulated relationships (lines with open symbols) using the generic parameters from Table 1, for qc1Ncs values of 70; 100 and 140, and σ’vo from 20 to 400 kPa. The two sets of curves (target and SDM simulated) show a reasonably good agreement, with the SDM simulated curves manifesting a slightly stronger effect of confining stress on liquefaction resistance, particularly near the boundary values of the qc1Ncs principal range of interest (i.e. qc1Ncs � 50 150) where the largest discrepancies occur. However, it is important to note that the trends with respect to the influence of qc1Ncs on Kσ σ ’vo are well captured in these simulations. In summary, the above calibration of the SDM model provides a reasonably accurate modelling of target liquefaction resistance curves over the relevant range of cyclic shear stresses for earthquake engi­ neering (Nc ¼ 1 30 cycles) including the threshold CSR separating between liquefaction and no-liquefaction, and across all densities or penetration resistances (qc1Ncs � 50 150) and confining stresses of interest (σ ’vo � 20 400 kPa). The above are considered key re­ quirements for the calibration of any constitutive model for dynamic analysis targeting liquefaction problems. Finally, note that qc1Ncs and σ ’vo are the only input parameters required by the model.

Fig. 7. Comparison between empirical Dr used in SDM calibration.

3.4. Modelling of non-liquefiable soils (step 3) 3.4.1. Identification of “non-liquefiable” soil layers Natural soil deposits generally have soil layers that are, by compo­ sition, susceptible to liquefaction and others that are not. Ideally, index and cyclic testing of retrieved soil samples in the laboratory should be used to distinguish between liquefiable and non-liquefiable soils (e.g. Ref. [8]). In the absence of soil sample data, the soil behaviour type index Ic or the modified soil behaviour type index IB [26] obtained from CPT measurements may be used for classification of soil behavior type. As a general rule, in the CPT-based effective-stress analysis approach of this paper, soil layers with Ic < 2:6 are considered liquefiable and they are modelled as described in the previous section; all other soil layers are regarded as non-liquefiable and their modelling is discussed in this section. Exceptions to this rule apply for liquefiable soils (Ic < 2:6) that are not expected to respond in a strongly non-linear manner or develop large excess pore water pressures (γmax � e 0:1 0:2%). These might be liquefiable (by composition) soils with high qc1Ncs (e.g. qc1Ncs > 170) or soils at large depths where (and when) the seismic demand is relatively low. Modelling of such soils and conditions must be targeting soil behaviour for strains γ � e 0:1 0:2%, rather than their liquefaction resistance. Taking this into consideration, the modelling strategy described in this section may offer a better alternative also for such soils. 3.4.2. Target stress-strain behaviour The nonlinear shear stress–shear strain (τ γ) response of nonliquefiable soils at a constant effective-stress can be approximated by the small strain shear modulus (Gmax ), the peak shear stress (τmax ) at large strains, and a strain-dependent modulus that describes the tran­ sition (change in stiffness) from small strains to large strains. Under cyclic loading, a measure for the energy dissipated in each loading cycle is also required. It is common in site response analysis to represent the strain-dependent variation in stiffness and energy dissipation via secant stiffness degradation (G=Gmax ), and damping ratio (ξ) curves, such as those depicted in Fig. 9. In this study, the stiffness degradation and damping ratio models proposed by Darendeli [16] were initially adopted as target curves. Key parameters required for the definition of G=Gmax γ and ξ γ curves in Darendeli’s model are the mean effective confining stress (σ’mo ¼ ½ð2Ko þ 1Þ =3�σ ’vo , where Ko is the earth pressure coeffi­ cient at rest), the overconsolidation ratio (OCR) and the soil plasticity index (PI). More often than not, some adjustment of the original Dare­ ndeli G=Gmax γ model is necessary to make the backbone stress-strain curve asymptotically approach the target shear strength (τmax ) of the soil

qc1Ncs relationships and Eq. (12)

8

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

procedure is required in which constitutive model parameters that provide reasonably accurate simulation of the target stress-strain behavior are determined. Any non-linear constitutive model that can concurrently model both the target G=Gmax γ and the target ξ γ curves could be used. In fact, this is often difficult to achieve because the Masing rule [24] that most constitutive models adopt for modelling the hysteretic behaviour of soils during cyclic loading does not necessarily match the corresponding behaviour observed in the lab. Recognizing this limitation, a compromise approach often needs to be taken, where both the target damping and modulus curves are reasonably fitted over the expected range of strains. Fig. 9 shows an example of stiffness degradation and damping ratio curves simulated using the SDM model. It can be seen that the SDM simulation slightly overpredicts shear stiffness in a range of shear strains from 0.001% to 1%, whereas it significantly overpredicts the damping ratio for shear strains greater than about 0.5%. Provided that the maximum shear strains in the nu­ merical simulations do not exceed (about) 0.5%, such compromise modelling is considered acceptable, as the model is expected to be used only over the range of strains where it shows good performance. To verify the above approach for modelling non-liquefiable soils, a series of total stress analyses of Christchurch sites was performed and the results were compared with analogous simulations using the 1-D nonlinear site response analysis program DEEPSOIL [18]. DEEPSOIL was chosen because it provides a more rigorous modelling option (MRDF) for hysteretic stress-strain response, with a nearly perfect fit for both G=Gmax γ and ξ γ target curves even at large strains. A com­ parison between the ground motion and shear strain predictions from total stress analyses with SDM in DIANA-J (DJ-SDM) and GQ/H model with MRDF in DEEPSOIL (DS-GQ/H þ MRDF), at different depths of an

Fig. 9. Example target shear modulus degradation and damping ratio curves (solid lines), together with fitted curves from SDM simulations (dashed lines).

at large strains [37]. In summary, the complete stress-strain (target) model for non-liquefiable soils can be defined if Gmax , τmax , K, OCR and PI are known. The reader is herein referred to the CPT guide of Rob­ ertson and Cabal [27] for a detailed review of available CPT empirical correlations that can be used to infer (most of) the above parameters. As there is no reliable correlation to estimate PI directly from CPT data, PI has to be assumed based on geological or other available information. 3.4.3. Constitutive modelling Once the target stress-strain model has been defined, a calibration

Fig. 10. Comparison between total-stress site response analysis results using SDM in DIANA-J (DJ-SDM; dashed lines) and GQ/H model with MRDF in DEEPSOIL (DS-GQ/H þ MRDF; solid lines): (a) acceleration time-histories and (b) spectral accelerations at various depths, and (d) maximum shear strains throughout depth; Vs profile of the site is shown in (c). 9

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

example site, is presented in Fig. 10. It can be seen that the two analyses predict almost identical ground motions but the DJ-SDM analysis, as expected due to the overprediction of shear modulus and damping, predicts slightly smaller strains than DS-GQ/H þ MRDF. Nevertheless, the SDM performance in total stress analysis is, clearly, sufficiently good for the purpose of this study. 3.5. Definition of input ground motions (step 4) Several methods for ground motion selection have been proposed in the literature over the last decades (e.g. Ref. [33]). Among them, the generalized conditional intensity measure (GCIM) approach [4,5] allows to rigorously consider multiple intensity measures (IMs) in the selection. Thus, in the selection of ground motions for liquefaction analyses, the GCIM approach allows for concurrent consideration of amplitude-, duration-, and energy-related IMs which may be equally important. The GCIM approach for ground motion selection can be broken down into three main steps. The first step entails dis-aggregation of the seismic hazard curve to obtain the contribution from different sources and earthquake events at a given IM level, referred to as the “conditioning IM” (and denoted as IMj ). The second step involves derivation of the marginal conditional distributions of each (single) considered IM (IMi ) taking into account the contribution of all scenario ruptures to the seismic hazard at the conditioning IM level (IMj ). The obtained marginal IMi distributions are used to generate realizations of the multivariate IM distribution considering the correlation between the considered IMs, which are then used to assess the appropriateness of the candidate ground motions. In the final step, a database of prospective ground motions is searched to find ground motions that fit best the generated realizations of the IM distribution. A weight vector, wi , is used to pre­ scribe the relative importance of the considered IMi and calculate the misfit of each prospective ground motion with respect to the target distribution. Bounds on causal parameters (e.g. magnitude, source-tosite distance, site condition) of prospective ground motions may also be considered prior to conducting IM-based ground motion selection [32]. The main steps described above refer to the case in which the seismic hazard is defined in terms of a seismic hazard curve resulting from a probabilistic seismic hazard analysis (PSHA). An analogous process for ground motion selection can be followed in the case of sce­ nario earthquake ruptures [31]. The GCIM approach (including random realizations of the IM dis­ tribution) has been implemented in the open-source software for seismic hazard analysis OpenSHA [5,17]. Alternatively, a set of Matlab codes that can be used for computing the GCIM distributions and subse­ quently, selecting appropriate ground motions is available from Ref. [7].

Fig. 11. Schematic illustration of the soil-column model used in numeri­ cal analyses.

1 βVs ΔY1 ¼ λmin ¼ 8 8fmax

(15)

ΔY2 ¼ 3ΔYadj

(16)

ΔY3 ¼ 2:0 m

(17)

where Vs is the shear wave velocity of the considered layer (appropri­ ately reduced to account for stiffness degradation at large strains via a reduction factor β), λmin is the wavelength associated with the highest frequency that is considered in the analysis (fmax Þ, and ΔYadj is the larger thickness of the two adjacent elements. The condition ΔY � ΔY1 (Eq. (15)) ensures accurate transmission of waves with frequency up to fmax [22] whereas the condition ΔY � ΔY2 (Eq. (16)) is adopted to avoid disproportionate changes in the thickness of adjacent elements. For liquefiable soil elements, to avoid impractically small ΔY1 values due to a potential large degradation in Vs , a maximum element thickness of 0.6 m is adopted, irrespective of βVs . The above element size constraints dictate the sizing and number of elements required in the discretization of each layer of the simplified profile.

3.6. Definition of numerical model and analysis parameters (step 5) In the proposed numerical procedure, fully coupled nonlinear effective-stress analyses are performed using 1-D soil column models (i. e. 1-D vertical wave propagation with 2-D quadratic elements con­ strained to deform in simple shear) to simulate the free-field response at level ground sites. The soil-column models are developed based on the simplified soil profile and constitutive model parameters defined in steps 1 to 3. In the following, salient features of the numerical model and analysis are discussed. 3.6.1. Mesh geometry A sketch of the adopted mesh geometry is illustrated in Fig. 11. The sizing of the mesh elements is a trade-off between model accuracy and computational efficiency. The minimum element thickness (ΔY) is taken equal to the minimum layer thickness adopted for the simplified profile (i.e. ΔYmin ¼ Tmin ¼ 0:3 m). The maximum ΔY is taken as the smallest of ΔY1 , ΔY2 and ΔY3 , which are defined as (see also Fig. 11):

3.6.2. Initial stress conditions, boundary conditions and load application Prior to the dynamic analysis, the initial stress state in the model has to be established. This can be done either by simply calculating the static stresses in each element and imposing them “externally” as an initial condition onto the model, or by conducting a self-weight gravity 10

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

analysis for the soil-column. In this analysis, the displacements of the nodes at the base of the column are fixed in both the x- and y-directions whereas the remaining soil nodes are fixed in the x-direction only. In the subsequent dynamic analysis, nodes at the same elevation are free in the x-direction but they are tied to share same displacement, thereby enforcing a simple shear mode of deformation. A Lysmer and Kuhle­ meyer [23] dashpot element with dashpot coefficient c is employed at the base of the soil-column to simulate the compliance of the underlying half-space, as schematically illustrated in Fig. 11. The soil-column is excited at the base by a horizontal force time history f which is pro­ portional to the known velocity time history of the input ground motion _ (u). Theoretically, the soil-column must be extended deep enough to include an elastic bedrock layer at the bottom. Because the exact soil conditions down to bedrock depth may often not be known, it is prac­ tically sufficient to place the base of the soil column at a shallower layer, which is stiff compared to the near-surface soils and is expected to respond with nearly elastic deformations. In this case, the ground mo­ tion selection process described in step 5 should assume soil conditions representative for the underlying half-space.

the free-field response at two Christchurch sites: one in the Hoon Hay suburb and the other in the Avondale suburb. The two sites exhibited vastly different performance during the 2010-2011 Canterbury earth­ quakes. The Hoon Hay site did not show any evidence of liquefaction at the ground surface (neither sand boils nor cracks or fissures) in any of the Canterbury events. In contrast, the Avondale site manifested severe liquefaction in the Mw6.2 22nd February 2011 earthquake, with thick soil ejecta covering the majority of land and roads in the area. This site also exhibited moderate levels of liquefaction in a number of other earthquakes in the series. Moderate-to-severe lateral spreading occurred in this area during the February earthquake. In the following, we examine the response of the two sites due to the 22nd February 2011 earthquake using the proposed procedure. The lo­ cations of the sites with respect to the surface projection of the source faults for the February event is shown in the map of Fig. 12. Both sites are in close proximity to the seismogenic fault (less than 4 km distance). Estimates of the surface peak ground acceleration (PGA) at these sites during the February earthquake were made using the methodology outlined in Ref. [6], and gave a median PGA of 0.46 g for the Hoon Hay site and 0.40 g for the Avondale site. Interestingly though, despite being subjected to a similar shaking intensity, the two sites showed a vastly different performance. This difference in the performance is attempted to be explained through the use of the proposed effective-stress analysis procedure. For comparison purposes, results from simplified liquefac­ tion analyses are first discussed.

3.6.3. Drainage conditions and permeability values In the effective-stress analyses, the soil is treated as a two-phase medium based on Biot’s equations for dynamic behaviour of saturated porous media [1]. The analyses are performed assuming drained con­ ditions allowing for pore water redistribution and vertical water flow through and between layers. Horizontal water flow is restricted in the 1-D soil-column analysis. Permeability values (k) for each distinct soil layer can be approximated from the Ic values of the simplified profile using Eq. (18) [27]: � ð0:952 3:04I Þ c 10 ; 1:0 < Ic � 3:27 k¼ ð ​ m = sÞ (18) 10ð 4:52 1:37Ic Þ ; 3:27 < Ic < 4:0

4.2. Liquefaction evaluation using the simplified method Results from simplified liquefaction analyses of the Hoon Hay and Avondale sites are presented in Figs. 13 and 14, respectively. The shallow soil profiles in the top 10 m of the two sites, including simplified Ic and qc1Ncs profiles, and characteristic soil units throughout depth are shown in the three leftmost columns of these figures. Soil units include clean sands (1:3 < Ic � 1:8, denoted as “C. sand” in the figures); sands with small amounts of fines (1:8 < Ic � 2:1, denoted as “F. sand”); sandy silts and non-plastic silts (2:1 < Ic � 2:6, denoted as “NP silt”); and nonliquefiable plastic silt/clayey/peat soils (Ic > 2:6, denoted as “PL soil”). Earthquake-induced cyclic stresses (CSR), cyclic resistances (CRR) and factors of safety against liquefaction (FS) were calculated for each layer using the deterministic liquefaction triggering procedure of Boulanger and Idriss [3]. Magnitude scaling (MSF) and overburden correction (Kσ ) factors were applied to the earthquake-induced CSR; therefore, both CSR and CRR in Figs. 13 and 14 are adjusted to the reference condition of Mw ¼ 7:5 and σ’vo ¼ 1 atm. Maximum shear strains (γ max ) were esti­ mated from the Yoshimine et al. [38] expressions, as reported in Ref. [20] and using the relative density Dr qc1N correlation proposed therein.

3.6.4. Analysis parameters and output In the final step of the procedure, considering the adopted element sizing and anticipated behaviour, analysis parameters such as compu­ tational time increment, integration scheme and numerical damping are adopted, and the dynamic effective stress analysis is then executed. In the example analyses presented in the following sections, an implicit Newmark method (β ¼ 0:25 and γ ¼ 0:5) was used for time integration, with a computational time step of Δt ¼ 0:005 s and Rayleigh damping with proportionality constants α ¼ 0 and β ¼ 0:005 applied to the mass and stiffness terms, respectively. These parameters were actually shown to provide valid results for a large number of analyzed soil models. Analysis output can be extracted in terms of time-histories of accel­ eration, velocity or displacement at any soil node including at the ground surface, as well as stresses, strains, and pore pressures deter­ mined for soil elements. Note that post-liquefaction reconsolidation volumetric strains (εvc;max ), associated free-field settlements and strainbased liquefaction damage indices can be estimated in a postprocessing phase by utilizing the computed maximum shear strains γ max from the dynamic effective-stress analysis as input in empirical εvc;max γmax relationships [38]. In essence, all key measures of ground response determined from a simplified liquefaction analysis can be easily obtained from the results of the proposed numerical procedure, but the latter would also provide temporal and spatial evolution of the response, while accounting for important dynamic interactions in the response of the soil deposit.

4.2.1. Hoon Hay site The Hoon Hay site (Fig. 13) has three potentially liquefiable zones: a non-plastic silt zone immediately below the groundwater table from 0.8 to 1.4 m depth, a clean sand to fines-containing sand zone from 3.0 to 4.2 m, and a deeper non-plastic silt layer from 4.8 to 5.8 m depth. The three liquefiable zones, which are separated by non-liquefiable layers (Ic > 2:6), are indicated with light shading in the qc1Ncs column of Fig. 13. Simplified liquefaction triggering analysis using the estimated PGA at the site indicates that liquefaction was triggered in all three liquefiable zones during the 22nd February 2011 earthquake. The deepest liquefiable layer at 4.8 m has the lowest cyclic resistance, yields the lowest factor of safety against liquefaction and develops the largest maximum shear strain, however the two other zones in which lique­ faction was predicted by the simplified analysis are perhaps more rele­ vant for liquefaction manifestation at the ground surface due to their shallow depth (0.8 and 3.0 m) and hence, close proximity to the ground surface.

4. Example applications 4.1. Selected sites and their seismic performance The CPT-based effective-stress analysis framework described in the previous sections is herein applied to investigate the characteristics of 11

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Fig. 12. Locations of the two analyzed sites (Hoon Hay and Avondale sites) and the reference site (CACS) used for the deconvolution analysis, with respect to the location of the finite fault planes of the Mw6.2 22nd February 2011 earthquake.

Fig. 13. Soil profile and characteristic results from simplified analysis at the Hoon Hay site, shown for the top 10 m of the deposit.

12

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Fig. 14. Soil profile and characteristic results from simplified analysis at the Avondale site, shown for the top 10 m of the deposit.

4.2.2. Avondale site The Avondale site (Fig. 14) is comprised of vertically continuous liquefiable (by composition) soils in the top 10 m of the soil profile. Silty sands and sandy silts are present above the groundwater table whereas clean sands prevail from 1.9 to 10 m depth. Two sand layers with small amount of fines are also present, one immediately below the water table from 1.9 to 2.6 m depth from the ground surface, and the other from 6.8 to 7.4 m depth. These two layers have the lowest cyclic resistance in the deposit and are again marked with the light shading in the qc1Ncs column of Fig. 14. Liquefaction is predicted to occur in both of these layers as well as at other depths of the clean-sand part of the deposit. More spe­ cifically, liquefaction triggering is predicted in two separate zones, from 1.9 to 3.2 m and from 4.9 to 7.4 m depth. Marginal liquefaction of thin layers is also predicted at greater depths. Clearly, in this case the critical layer for liquefaction manifestation at the ground surface is the finescontaining layer from 1.9 to 2.6 m depth; according to the results from the simplified analysis, this is the shallowest liquefied layer that develops the largest maximum shear strain in the deposit. Overall, the simplified analyses predict that liquefaction should have occurred in the deposits of both sites including at shallow depths. In both cases, the critical layers for liquefaction manifestation were predicted to liquefy and develop large strains. Evaluation of liquefaction damage indices using the results from the simplified analyses yields similar LPI values of 13 and 12 and similar LSN values of 22 and 21 for the Hoon Hay and Avondale sites, respectively, indicating moderate-to-severe liquefaction damage for both sites. Yet, the two sites exhibited highly contrasting performance in the 22nd February 2011 earthquake, i.e. no liquefaction manifestation (Hoon Hay site) and severe liquefaction manifestation (Avondale site).

4.3. Results from effective-stress analyses In the effective stress analyses of the two sites, a ground motion recorded during the 22nd February 2011 earthquake was used in order to closely represent the ground shaking during this event. This is an alternative approach to that discussed in the previous sections for defining the input ground motion that can be used for assessment of a site performance in past events when more appropriate ground motion recordings are available. More specifically, the fault-normal component of the surface ground motion recorded at the Canterbury Aero Club (CACS) strong motion station (a stiff and shallow gravel site that exhibited only minor levels of nonlinear response during this event, see location in Fig. 12) was deconvoluted to the top of the Riccarton Gravel layer (i.e., a firm and dense gravel layer underlying the shallow surface deposits of Christchurch, typically encountered at a depth of about 10 to 40 m below the ground surface, in the western and eastern parts of the city). The deconvoluted ground motion was then scaled to PGA ¼ 0:40 g (scaling factor � 2:1) to approximately represent the intensity of shaking near the fault, at the locations of the two considered sites [25]. This scaled motion was used as a base excitation of the soil-column models of the two sites at a depth of 20 m from the ground surface. 4.3.1. Hoon Hay site Fig. 15 shows characteristic acceleration and excess pore water pressure (ue ) time histories computed at various depths in the analysis of the Hoon Hay site. In fact, the acceleration and ue time series presented in Fig. 15 show the responses of the three key liquefiable layers identi­ fied in the simplified liquefaction analyses discussed in the previous section. Response characteristics throughout the depth of the deposit including the variation of the seismic demand in terms of the Arias in­ tensity of the ground motion (Ia ), excess pore water pressures at specific 13

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Fig. 15. Effective-stress analysis results at selected depths for the Hoon Hay site: (a) accelerations, and (b) excess pore water pressure time histories.

time sections, and computed maximum shear strains are shown in Fig. 16, throughout the top 10 m of the deposit. Initial analyses were conducted assuming fully saturated conditions for soils below the groundwater table. Subsequent examination of compression wave velocity (Vp ) measurements that were available at the site from detailed cross-hole testing [10] indicated that the three liq­ uefiable zones of the soil profile may actually be partially saturated. Partial saturation indicated by measured Vp values substantially below 1500 m/s was often encountered in highly interbedded deposits of liq­ uefiable and non-liquefiable soils. As indicated in the Ic column of Fig. 16, Vp varied from 700 to 1100 m/s in the two deeper liquefiable zones, from 3.0 to 5.8 depth, and was less than 400 m/s in the top 1.4 m of the deposit. To account for effects of partial saturation on the cyclic resistance of liquefiable soils, the effective-stress analyses were repeated using now an increased CRR for partially saturated layers according to the Vp -based model proposed by Hossain et al. [19]. This was achieved by making the dilatancy strain parameter Sc of SDM a function of the measured Vp of each layer: Sc ¼ 0:03⋅1:005

Vp

þ 0:0040

Figs. 15 and 16 refer only to the analysis considering partial saturation effects. The analysis results depicted in Figs. 15 and 16 show that the deepest of the three key layers (zones) in the deposit, at about 5 m depth, rapidly develops excess pore water pressures and liquefies within the first 3 s of strong shaking (from 10.5 to 13.5 s). This early liquefaction results in large maximum shear strains of about 5% in this layer. Notice that the response of this layer is (almost) identical for the full-saturation and the partial-saturation analysis cases (Fig. 15b). As apparent in Fig. 15a and especially in the Ia column of Fig. 16, the early liquefaction of the loose silt layer from 4.8 to 5.8 m depth caused a substantial reduction in the accelerations of this layer as well as in the accelerations transmitted to the overlying soils, above 4.8 m depth. In the full-saturation analysis case, this liquefaction-induced reduction in the seismic demand simply delays the triggering of liquefaction in the shallow layer from 3.0 to 4.2 m (Fig. 15b). However, in the partialsaturation analysis case, the decrease in the seismic demand in conjunction with the increase in CRR due to effects of partial saturation effectively prevented occurrence of liquefaction in the shallow layer, from 3.0 to 4.2 m depth. In both analysis cases, liquefaction did not occur in the shallow non-plastic silt layer from 0.8 to 1.4 m depth. For the sake of completeness, results from simplified analysis at the Hoon Hay site with consideration of the effects of partial saturation on CRR are presented in Fig. 18. In this case, marginal liquefaction is pre­ dicted in the shallow layer from 0.8 to 1.4 m, resulting in reduced maximum shear strains compared with those predicted in Fig. 13, but practically no effect is seen in the deeper layers. This further supports the hypothesis that indeed the combined effects of the two mitigating

(19)

All other SDM parameters remain unchanged, as in Table 1, Eq. (19) results in an increase in the simulated CRR (as a function of Vp ) which is in reasonable agreement with that predicted by the model of Hossain et al. [19], as shown in Fig. 17. Eq. (19) is an extension of the proposed methodology for modelling of liquefiable soils (Step 2) that allows for incorporation of partial saturation effects on LRC. The ue time histories in Fig. 15 are shown for both fully saturated and partially saturated analysis cases whereas all other results presented in 14

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Fig. 16. Soil profile and response characteristics from effective-stress analysis at the Hoon Hay site, shown for the top 10 m of the deposit.

Fig. 17. Modelling of partial saturation effects on CRR based on Vp : SDM simulations using Eq. (19) compared with the regression model of Hossain et al. [19].

mechanisms (i.e., reduction in the seismic demand due to liquefactioninduced “isolation” effect and partial saturation at shallow soils) pre­ vented liquefaction of the shallow layers in the effective-stress analysis. 4.3.2. Avondale site A similar format is used in Figs. 19 and 20 to present the results from the effective-stress analyses of the Avondale site. Fig. 19 shows accel­ eration and excess pore water pressure time histories of key layers in the deposit, whereas Fig. 20 summarizes the analysis results in soil-column

Fig. 18. Characteristic results from simplified analysis at the Hoon Hay site, considering the effects of partial saturation on cyclic resistance. 15

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Fig. 19. Acceleration and excess pore water pressure time histories (ue ) at selected depths for the Avondale site: (a) accelerations; and (b) excess pore water pressures.

plots of key response measures (Ia , ue and γmax ) for the top 10 m of the soil profile. In this analysis, liquefaction is first triggered at about 13.0 s (of the computational time), almost concurrently in the two finescontaining sand layers with the lowest cyclic resistance, at 1.9 and 6.8 m depth respectively. The computed maximum shear strains in these layers are consistent with the early development of liquefaction and are slightly larger in the deeper layer. Note that at the Avondale site, only the shallow critical layer is partially saturated with a Vp of about 1000 m/s. Nevertheless, the slightly increased CRR due to partial saturation effects does not affect the development of excess pore water pressures in this layer, as demonstrated in Fig. 19b. Liquefaction of the deeper layer from 6.8 to 7.4 m depth has two opposing effects on the response of this site. On one hand, as in the case of the Hoon Hay site, it reduces the demand for all overlying soils, as indicated by the abrupt drop in Ia at this depth (Fig. 20). Notice however that, due to the slightly higher penetration resistance (i.e. higher relative density) and smaller thickness of the liquefied layer, this effect is not as pronounced as in the case of the Hoon Hay site. On the other hand, the substantial difference in the excess pore water pressures (and conse­ quently in the total hydraulic head) between the liquefied layer and the surrounding higher resistance soils, at the time of triggering and following the initial liquefaction, creates a high hydraulic gradient outward from the liquefied layer or towards the surrounding soils. Given the high permeability of the surrounding soils, this gradient implies strong flow of water out of the liquefied layer. The effects of this water flow can be seen in the gradual increase of excess pore water pressures within the surrounding soils, i.e., below and above the deep liquefied

layer (see for instance, the gradual increase in ue at z ¼ 3:7 m depth in the time history plot in Fig. 19b). Actually, as shown in Fig. 19b, the excess pore water pressures within the deeper layer (z ¼ 7:3 m) remain at the level of the initial effective overburden stress (σ’vo ) for a long time after the onset of liquefaction, thereby continuing the supply of water to the shallower soils even after the end of strong shaking. This mechanism, together with contributions from shaking-induced pore pressure gener­ ation and additional in-flow of water from the shallow liquefied layer, eventually result in liquefaction of an enlarged thicker critical zone at shallow depth, from 1.9 to about 4.6 m, which includes the critical layer but also extends to deeper layers of much higher resistance. A contrib­ uting factor for the formation of such vertically continuous liquefied zone at shallow depth and the prolonged liquefaction of the critical layer seems to be the presence of the lower permeability silt layer above the initial water table depth, which slows down the dissipation of the developed excess pore water pressures towards the ground surface. In fact, it is highly likely that this surface layer will also liquefy at a later stage of the reconsolidation process, extending further upwards the critical liquefied zone towards the ground surface (Fig. 20). 4.4. Discussion Careful examination of the effective-stress analysis results at the two investigated sites reveals some important aspects of the system-response of liquefiable deposits. The response at the Hoon Hay site was primarily affected by dynamic cross-interaction between layers, where liquefac­ tion first occurred in a loose relatively deep layer resulting in substantial 16

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Fig. 20. Soil profile and characteristic response profiles obtained from effective-stress analysis at the Avondale site, shown for the top 10 m of the deposit.

reduction of the demand for all soils above that depth. This onset of deep liquefaction in conjunction with partial saturation effects at the shallow part of the deposit prevented development of liquefaction in the shallow critical layers and consequent surface liquefaction manifestation, by effectively creating a thick non-liquefied crust from the ground surface to about 5 m depth. Liquefaction of the 1 m – thick layer at 5 m depth would be unlikely to manifest at the ground surface for the seismic de­ mand imposed by the Mw6.2 22nd February 2011 earthquake, at this site. At the Avondale site, the system response was significantly influenced by interactions between layers through pore water pressure re-distribution and water flow. In particular, water flow from early liquefied layers to the surrounding soils, in combination with layer permeability contrasts that impeded fast drainage and dissipation of the excess pore water pressures, contributed to the formation of a thick (partly shakinginduced and partly seepage-induced) liquefaction zone at shallow depth, which could not have been predicted by an undrained analysis without groundwater flow. Indeed, excess pore water pressures are very high (at their maximum or close to it) and increasing with depth throughout the sand deposit indicating that the entire deposit will contribute to the strong water flow. This mechanism involving vertical communication of excess pore water pressures and large volumes of upward flowing water can explain the severe liquefaction effects that were observed at this site after the 22nd February 2011 earthquake. Analyses of the two sites using the simplified procedure were unable to detect the critical cross-layer interaction mechanisms and quantify their significant effects on the response, thereby leading to inaccurate predictions of the severity of surface liquefaction manifestation and associated damage. A more systematic study on the system response of

liquefiable deposits can be found in Ref. [14]. Uncertainties in both the simplified and effective-stress analyses presented herein have to be recognized, particularly those related to the seismic demand which is used as input in these methods. While explicit consideration of these uncertainties is beyond the scope of these example applications, preliminary sensitivity analyses have shown that key observed characteristics and mechanisms in the response are prac­ tically not affected by these uncertainties. 5. Conclusions A CPT-based effective-stress analysis procedure for liquefaction assessment has been presented. Key elements in the proposed procedure can be summarized as follows: (1) As a first step, a simplified, discretized soil profile is determined from the nearly continuous CPT data. A practical algorithm that automates this process has been presented in this paper. Apart from the obvious necessity of a discretized soil profile for the numerical analysis, this profile “simplification” allows for rigorous scrutiny of the overall site characteristics and identification of critical layers in the deposit. Based on their inferred behaviour characteristics, the soil layers resulting from this process are classified as either liquefiable or non-liquefiable. (2) Modelling of the liquefiable soil layers focuses on the simulation of their (cyclic) liquefaction resistance, while using representa­ tive values for elastic and plastic stress-strain parameters for sand. A set of liquefaction resistance curves (LRCs) were derived 17

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

Acknowledgements

over relevant qc1Ncs σ ’v conditions by directly following the simplified liquefaction triggering procedure of Boulanger and Idriss [3]. These target curves were then used to calibrate a constitutive model capable of reproducing the target behavior over all densities and confining stresses of interest, with a single set of values for model parameters. (3) Modelling of the non-liquefiable soil layers targets their cyclic stress-strain response, typically defined in terms of stiffness degradation and damping ratio curves. Reasonably accurate and concurrent modelling of target stiffness degradation and damping ratio curves over the expected range of shear strains is the key requirement in this case. (4) With respect to the seismic input, it is recommended to select an ensemble of ground motions for the earthquake scenario of in­ terest using the generalized conditional intensity measure (GCIM) approach [4,5], considering a suitable set of intensity measures (IMs) relevant to liquefaction problems, including not only amplitude-related IMs, but also duration- and energy-related measures, as they are also critical for liquefaction problems. (5) In the final step, numerical model (i.e. element size, boundary conditions, initial stress state of the soil) and analysis parameters (i.e. computational time increment, integration scheme, numer­ ical damping) are defined. Here, basic rules of a good numerical analysis need to be followed, always taking into consideration the characteristics of the given soil profile and its anticipated behaviour.

The first author would like to acknowledge the EQC-QuakeCoRE scholarship. This is QuakeCoRE publication number 0484. References [1] Biot MA. Theory of propagation of elastic waves in a fluid saturated porous solid, part I–low frequency range; part II–high frequency range. J Acoust Soc Am 1956;28 (1):168–91. [2] Boulanger RW, DeJong JT. Inverse filtering procedure to correct cone penetration data for thin-layer and transition effects. In: Proc 4th international symposium on cone penetration testing, Delft; 2018. [3] Boulanger RW, Idriss IM. CPT and SPT based liquefaction triggering procedures. Report No. UCD/CGM-14/01. Davis: Center for Geotechnical Modeling, Department of Civil and Environmental Engineering, University of California; 2014. [4] Bradley BA. A generalized conditional intensity measure approach and holistic ground-motion selection. Earthq Eng Struct Dynam 2010;39(12):1321–42. [5] Bradley BA. A ground motion selection algorithm based on the generalized conditional intensity measure approach. Soil Dynam Earthq Eng 2012;40(1): 48–61. [6] Bradley BA. Site-specific and spatially-distributed ground-motion intensity estimation in the 2010-2011 Canterbury earthquakes. Soil Dynam Earthq Eng 2014;61-62(1):83–91. [7] Bradley BA. https://sites.google.com/site/brendonabradley/software/ground-mo tion-selection-gcim. [Accessed 20 July 2019]. [8] Bray JD, Sancio RB. Assessment of the liquefaction susceptibility of fine grained soils. J Geotech Geoenviron Eng 2006;132(9):1165–77. [9] Ching J, Wang J-S, Hsein Juang C, Ku C-S. Cone penetration test (CPT)-based stratigraphic profiling using the wavelet transform modulus maxima method. Can Geotech J 2015;52(1):1993–2007. [10] Cox BR, Stolte AC, Stokoe KH, Wotherspoon LM. A direct-push crosshole (DPCH) test method for the in situ evaluation of high-resolution P- and S-wave velocities. Geotech Test J 2018. available online, 10.1520/GTJ20170382. [11] Cubrinovski M. Some important considerations in the engineering assessment of soil liquefaction. In: Proc. 2019 ANZ conference, Perth, Australia: 2019. ANZ Conference; 2019. p. 1–3. Apr 2019. [12] Cubrinovski M, Ishihara K. Modelling of sand behaviour based on state concept. Soils Found 1998;38(3):115–27. [13] Cubrinovski M, Ishihara K. State concept and modified elastoplasticity for sand modelling. Soils Found 1998;38(4):213–25. [14] Cubrinovski M, Rhodes A, Ntritsos N, van Ballegooy S. System response of liquefiable deposits. Soil Dynam Earthq Eng 2019;124(1):212–29. [15] Cubrinovski M, Uzuoka R, Sugita H, Tokimatsu K, Sato M, Ishihara K, Tsukamoto Y, Kamata T. Prediction of pile response to lateral spreading by 3-D soil-water coupled dynamic analysis: shaking in the direction of ground flow. Soil Dynam Earthq Eng 2008;28(1):421–35. [16] Darendeli MB. Development of a new family of normalized modulus reduction curves and material damping curves. PhD thesis. Austin: University of Texas; 2001. [17] Field EH, Jordan TH, Cornell CA. OpenSHA: a developing community-modeling environment for seismic hazard analysis. Seismol Res Lett 2003;74(1):406–19. [18] Hashash YMA, Musgrove MI, Harmon JA, Groholski DR, Phillips CA, Park D. DEEPSOIL 6.1, User manual. 2016. [19] Hossain AM, Andrus RD, Camp III WM. Correcting liquefaction resistance of unsaturated soil using wave velocity. J Geotech Geoenviron Eng 2013;139(2): 277–87. [20] Idriss IM, Boulanger RW. Soil liquefaction during earthquakes. Monograph MNO12. Oakland, CA: Earthquake Engineering Research Institute; 2008. p. 261. [21] Iwasaki T, Tatsuoka F, Tokida K, Yasuda S. A practical method for assessing soil liquefaction potential based on case studies at various sites in Japan. In: Proc 2nd international conference on microzonation. Washington, DC: National Science Foundation; 1978. 1978. [22] Kuhlemeyer RL, Lysmer J. Finite element method accuracy for wave propagation problems. J Soil Mech Found Div 1973;99(SM5):421–7. [23] Lysmer J, Kuhlemeyer RL. Finite dynamic model for infinite media. J Eng Mech 1969;95(EM4):859–77. [24] Masing G. Eigenspannungen and verfertigung beim messing. In: Proc 2nd international congress on applied mechanics, Zurich; 1926. [25] Ntritsos N, Cubrinovski M, Bradley BA. Effective stress analysis of Christchurch strong motion station sites. In: Proc 7th international conference on earthquake geotechnical engineering, Rome; 2019. paper no. 456. [26] Robertson PK. Cone penetration test (CPT)-based soil behaviour type (SBT) classification system – an update. Can Geotech J 2016;53(12):1910–27. [27] Robertson PK, Cabal KL. Guide to cone penetration testing. sixth ed. California: Gregg Drilling & Testing, Inc.; 2015. [28] Robertson PK, Wride CK. Evaluating cyclic liquefaction potential using the cone penetration test. Can Geotech J 1998;35(3):442–59. [29] Seed HB, Idriss IM. Simplified procedure for evaluating soil liquefaction potential. J Soil Mech Found Div 1971;97(SM9):1249–73. [30] Corporation Taisei. DIANA-J3: finite element program for effective stress analysis of two-phase soil medium. Internal report, Software science. 1997 (in Japanese).

Two important features of the proposed procedure are that it can be fully automated in a programming environment, and that it is directly equivalent (in the definition of cyclic resistance and required input data) to the CPT-based simplified procedures for liquefaction analysis. These features allow advanced effective-stress analysis to be routinely applied in practice, in parallel with simplified analysis and inform one another. The benefits from such an approach were illustrated through ana­ lyses of two case history sites from Christchurch that showed vastly different performance in the 2010-2011 Canterbury earthquakes. System-response mechanisms involving cross-interaction amongst different soil layers in the dynamic response and through pore water pressure redistribution and water flow were shown to govern the per­ formance of these sites. In contrast to the simplified analysis, the effective-stress analysis was capable of capturing these mechanisms and could therefore explain the differences in the observed manifestation. Limitations of both constitutive model and numerical analysis have to be recognized, including the fact that these limitations are often model-dependent and problem-dependent. Furthermore, it goes without saying that rigorous effective-stress analysis for accurate quantification of liquefaction effects on critical facilities cannot be based solely on CPT data. Additional data from field investigations and laboratory tests on retrieved soil samples must inform the characterization of soils and liquefaction assessment, in such cases. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Nikolaos Ntritsos: Conceptualization, Methodology, Software, Writing - original draft, Visualization. Misko Cubrinovski: Conceptu­ alization, Methodology, Resources, Writing - review & editing, Supervision.

18

N. Ntritsos and M. Cubrinovski

Soil Dynamics and Earthquake Engineering 131 (2020) 106063

[31] Tarbali K, Bradley BA. Ground motion selection for scenario ruptures using the generalized conditional intensity measure (GCIM) method. Earthq Eng Struct Dynam 2015;44(1):1601–21. [32] Tarbali K, Bradley BA. The effect of causal parameter bounds in PSHA-based ground motion selection algorithm. Earthq Eng Struct Dynam 2016;45(1): 1515–35. [33] Tarbali K, Bradley BA, Baker JW. Consideration and propagation of ground motion selection epistemic uncertainties to seismic performance metrics. Earthq Spectra 2018;34(2):587–610. [34] Tatsuoka F, Zhou S, Sato T, Shibuya S. Method for evaluating liquefaction potential and its application. Rep. on Seismic hazards in the soil deposits in urban areas. Ministry of Education of Japan; 1990. p. 75–109 (in Japanese). [35] Van Ballegooy S, Malan P, Lacrosse V, Jacka ME, Cubrinovski M, Bray JD, Cowan H. Assessment of liquefaction-induced land damage for residential Christchurch. Earthq Spectra 2014;30(1):31–55.

[36] Wang Y, Huang K, Cao Z. Probabilistic identification of underground soil stratification using cone penetration tests. Can Geotech J 2013;50(1):766–76. [37] Yee E, Stewart JP, Tokimatsu K. Elastic and large-strain nonlinear seismic site response from analysis of vertical array recordings. J Geotech Geoenviron Eng 2013;139(10):1789–801. [38] Yoshimine M, Nishizaki H, Amano K, Hosono Y. Flow deformation of liquefied sand under constant shear load and its application to analysis of flow slide in infinite slope. Soil Dynam Earthq Eng 2006;26(1):253–64. [39] Youd TL, Idriss IM, Andrus RD, Arango I, Castro G, Christian JT, Stokoe II KH. Liquefaction resistance of soils: summary report from the 1996 NCEER and 1998 NCEER/NSF workshops on evaluation of liquefaction resistance of soils. J Geotech Geoenviron Eng 2001;127(20):817–33.

19