Numerical Simulations of Selected LEAP Centrifuge Experiments with PM4Sand in FLAC

Increasing the confidence in estimates of liquefaction-induced deformations from numerical models requires careful validation of both constitutive models and numerical simulation approaches. This validation should be performed with proper consideration of the intended use of the simulation results and the uncertainty inherent in comparing numerical simulations to laboratory or physical models. This paper describes the comparison of a set of calibrated numerical simulations to the results of centrifuge experiments with a submerged, liquefiable slope. These simulations were performed as part of the 2017 LEAP simulation exercise using the constitutive model PM4Sand and the numerical platform FLAC. The simulation results were compared to those from the centrifuge to identify which aspects of the experimental response the simulations were adequately able to capture. This paper provides a description of the simulation approach including the constitutive model and numerical platform. The process of calibrating PM4Sand to the results of laboratory experiments is described. This is followed by a description of the system-level simulations and a comparison of the simulation results with the experimental results. A sensitivity study was performed to examine trends in the simulation response and a discussion of the findings from this study is presented.


Introduction
Liquefaction-induced ground failure continues to be a major source of damage to civil infrastructure during earthquakes. The ability of engineers to predict the response of infrastructure to liquefaction requires analysis techniques that can provide reliable estimates of the response of soils to seismic loading. For important structures, engineers are increasingly using nonlinear effective stress analyses that combine a constitutive model capable of reproducing important aspects of soil behavior with a finite element or finite difference platform that can capture the system response subjected to the applicable boundary and loading conditions. The quality of these analyses will at a minimum depend on the selection and calibration of the constitutive model, the ability of the numerical simulation approach to capture the important system-level responses, the quality of the soil characterization, and the quality of the documentation and review processes. Validation of the constitutive model and numerical platform is often accomplished by comparing results from numerical simulations to carefully performed and documented centrifuge tests. These comparisons provide a means to estimate the current capabilities and limitations of the numerical tools commonly used in liquefaction analyses, to identify major sources of uncertainty in both the numerical and experimental techniques, and to provide insights for future pathways.
This paper documents a set of numerical simulations performed as part of the 2017 LEAP simulation exercise. The simulations described in this paper used the constitutive model PM4Sand (Version 3, Boulanger and Ziotopoulou 2015) and the numerical platform FLAC (Version 8.0, Itasca 2016). The simulations are compared to selected experimental results with the goal of identifying aspects of the experimental response the numerical simulations were able to reasonably capture and which may require further investigation. A brief description of the simulation approach is presented, including the constitutive model and numerical platform used for both the calibration and system-level simulations. The calibration of the model to the results of laboratory experiments conducted on the designated soil for LEAP (Ottawa F-65 sand) is described. This is followed by a description of the simulations performed for the Type-B (Kutter et al. 2015) predictions of the centrifuge model tests, and a critical comparison between simulation and experimental results. A sensitivity study was performed to identify which of the input parameters had the largest effect on the response. A discussion of some of the findings from this study is presented.

Modeling Approach
The numerical platform used for these simulations was FLAC 8.0 (Itasca 2016). FLAC is a two-dimensional, explicit, finite difference program for engineering mechanics computation. The program simulates the behavior of structures built of soil, rock, or other materials that may undergo plastic flow when their yield limits are reached. Materials are represented by elements, or zones, which form a grid that is adjusted by the user to fit the shape of the object to be modeled. Each element behaves according to a prescribed linear or nonlinear constitutive stress/strain law in response to the applied forces and boundary conditions. FLAC uses an explicit, Lagrangian calculation scheme and a mixed-discretization zoning technique that allows for accurate modeling of plastic collapse and flow. Details of the formulation and validation of FLAC can be found in the software manual (Itasca 2016). The nonlinear constitutive model used for this study is PM4Sand Version 3 (formulation and implementation described in Boulanger andZiotopoulou 2015, compiled for FLAC 8.0 in 2017) which was specifically developed for earthquake engineering applications. PM4Sand is a stress-ratio controlled, critical state compatible, bounding-surface plasticity model, based on the plasticity model initially developed by Dafalias and Manzari (2004) and described in detail by Boulanger and Ziotopoulou (2015). It is cast in terms of the relative state parameter (ξ R ), which measures the difference between the relative density (D R ) and the relative density at critical state (D Rcs ) for the current confining pressure (Boulanger 2003). This formulation allows soil properties to change during the simulation as a function of the change in state (i.e., changes in mean effective stress and/or in void ratio). PM4Sand Version 3 has 22 input parameters, from which only three are considered primary and are required as model inputs. The other 19 (2 flags and 17 secondary parameters) can be left with their preset default values if no other information is available or calibrated to the desired response based on the available lab data. Pertinent information on the parameters of PM4Sand and their calibrated values are given in the following section, while extensive information on the model is provided by Boulanger and Ziotopoulou (2015) and Ziotopoulou and Boulanger (2016).

Calibration Approach
The goal of the constitutive model calibration performed for this study was to select parameters that can reasonably capture the experimental results, while considering the uncertainty inherent in using laboratory tests to predict centrifuge model test results. The calibration of any constitutive model must focus on the aspects of the model behavior which are important to the problem being analyzed (i.e., the constitutive behaviors that are reasonably expected to be activated by the loading paths of the problem at hand). The most critical aspects of the model response for the current study are expected to be the generation of pore pressure and triggering of liquefaction, the accumulation of strains during shaking, and the subsequent reconsolidation of the soil after shaking has stopped. The current calibration intentionally focused on capturing liquefaction triggering curves provided by the LEAP organizers ( Fig. 24.1) as well as some general aspects of the individual stress-strain and stress path responses for the cyclic triaxial test results (El Ghoraiby et al. 2019). For all other aspects of model behavior, PM4Sand parameters retained their default values that are generally functions of an index property (D R ) chosen to reasonably approximate design correlations.
Calibration of PM4Sand was performed by using single-element simulation of undrained cyclic stress-controlled plane strain compression (PSC) tests and comparing the results to the provided results from undrained cyclic stress-controlled triaxial tests on Ottawa F-65 sand (El Ghoraiby et al. 2017, 2019. The response of cyclic triaxial tests is, as expected, asymmetric with respect to deviatoric stresses. This is not the case with undrained cyclic PSC, so the goal of the calibration was to capture a reasonable overall match within each group of cyclic triaxial tests (each density) without looking to precisely match details from the individual tests.
Triaxial test results were provided for three void ratios (e o ¼ 0.585, 0.542, and 0.515). PM4Sand is cast in terms of D R , so it was necessary to estimate corresponding D R values for each of the tests. The organizers provided a database of index test results for Ottawa F-65 sand, which included minimum and maximum void ratios (e min and e max ) from various authors. There was considerable scatter in some of these properties and so some judgment was necessary in selecting single values. The authors chose to use an e max of 0.739 and an e min of 0.492 based on the results from George Washington University, which also performed the triaxial tests. This selection had a major impact on the comparison to the centrifuge results as will be discussed later. Using these values, the corresponding D R values for the experimental results were 62%, 80%, and 90%, respectively.

Single-Element Simulations
The undrained cyclic PSC simulations were performed using a single element (or zone) in FLAC 8.0. PM4Sand was assigned to the zone using the desired properties at the beginning of the simulation. Stresses were then initialized in the zone to replicate isotropic consolidation at an effective stress of 100 kPa. Cyclic loading was applied by applying a constant vertical velocity to the top boundary of the zone while fixing the bottom boundary against vertical movement. The left boundary of the zone was fixed against horizontal displacement while the right boundary had a constant horizontal force applied to represent the cell pressure. The direction of the velocity of the top boundary was reversed when the internal deviatoric stress reached the desired level as defined by the cyclic stress ratio (CSR ¼ q/2σ 0 vo ). The number of cycles required to reach 2.5% single amplitude axial strain and an excess pore pressure ratio (r u ) of 98% were tracked. Cyclic strength curves for each of the experiments are shown in Fig. 24.1. Two values of CSR are shown for the laboratory data in this figure. The first is the reported CSR provided by the organizers and the second is the average of the peak CSR achieved during each cycle of loading.

Calibrated PM4Sand Parameters
Calibration of PM4Sand requires, at a minimum, the three primary input parameters: (1) the apparent relative density (D R ) which controls the dilatancy and stress-strain response characteristics of the soil; (2) the contraction rate parameter (h po ) which controls the cyclic strength of the soil; and (3) the shear modulus coefficient (G o ) which controls the small strain stiffness of the model (G max ). For the current study, the parameter h po was iteratively adjusted for each D R until the liquefaction strength curves produced from the simulation reasonably matched the experimental data to the extent possible ( Fig. 24.1). The calibration reasonably tracked the pattern of responses across all the reported tests. The simulation results from the individual CSR levels were fit with a power function (CSR ¼ aN Àb ) which is illustrated as a solid line in the figures. The CSR to reach 2.5% single amplitude axial strain in 15 cycles was approximated for all three D R values ( Fig. 24.2a) using the results in Fig. 24.1. This curve is often called a liquefaction triggering curve and can also be used to calibrate h po for untested density levels. This relationship is later used to account for the variations with D R between the different centrifuge model tests.
The most challenging parameter to calibrate was the shear modulus coefficient (G o ), because no information on small strain stiffness or shear wave velocity was provided and G o was not well constrained by calibrating to the results from cyclic triaxial tests. Values for this parameter were obtained by comparing the simulation results to the stress-strain response from the cyclic triaxial tests during the first few cycles of loading. G o was iteratively changed until a reasonable approximation was achieved. Ziotopoulou (2018), in the predictions for LEAP 2015, had selected a G o of 240 for a D R of 65%, which was used as a starting point for this calibration. The selected G o values using this approach are shown in Fig. 24.2b. These values are lower than the default relationship used by PM4Sand would suggest (Fig. 24.2b), but follow the same linear trend with D R as expected. The effect of this difference will be examined in subsequent simulations.
Two of the secondary model parameters were adjusted to improve the comparison of the simulations results to the stress-strain and stress path results from the cyclic triaxial tests. The critical state friction angle φ cv was reduced to 30 from the default value of 33 , to better match the slope of the frictional envelopes (bounding line) in the stress path plots. This value is also more consistent with the values reported by Parra Bastidas (2016) based on a review of available literature. The second default parameter that was modified was n b which controls the bounding ratio and therefore dilatancy and peak effective friction angles. The default value of 0.5 was increased to 0.7 to reduce the rate of strain accumulation in the cyclic mobility regime following the triggering of liquefaction. Table 24.1 lists all model parameters for PM4Sand that were given values other than their default during the calibration and simulations. Default values for the remaining parameters are provided by Boulanger and Ziotopoulou (2015). The index properties that were utilized are: • The as-placed dry density or void ratio of as indicated by the reported results.
• The maximum and minimum void ratios of 0.739 and 0.492.
• The average specific gravity G s ¼ 2.66 as determined through tests (Vasko et al. 2014).

System-Level Numerical Simulations
Simulations were performed for each of the nine centrifuge tests completed as part of the LEAP 2017 exercise. All simulations used the same geometry which was based on the prototype dimensions ( Fig. 24.3). The numerical mesh consists of one uniform layer of sand, 80 zones wide and 16 zones high, yielding 81 gridpoints in the x-and 17 gridpoints in the y-direction (Fig. 24.4). Zones in the simulations presented herein had an average size of 0.25 Â 0.25 m, and a maximum aspect ratio of 2:1 at the bottom right corner of the grid. At the time of the Type-B predictions, no information was provided by the centrifuge facilities on the actual location of the instruments, so the goal of the discretization was to have grid points at the prescribed instrument locations.

Boundary Conditions
The mechanical boundary conditions in the simulations replicated the boundary conditions imposed by the rigid container used in the centrifuge model tests, without explicitly simulating the rigid box that surrounded the soil. All bottom nodes were fixed in the x-and y-directions and all side nodes were fixed in the x-direction. The side walls of the rigid container were assumed to have no friction and thus the soil was able to freely slide vertically during all stages of the simulation. Flow of water was allowed across the top surface of the model and restricted across the container boundaries. The experiments were submerged, so pore pressures and saturation were fixed at the top nodes and a pressure was applied to simulate the weight of the fluid. The values of the surface pressure from the water (both the external pressure and boundary pore pressure) were updated every 0.2 s during the simulation to account for any settlement of the soil surface. The elevation of the free water surface was assumed to be at 5 m. This value was determined using the initial hydrostatic values of the pore pressure histories reported by the centrifuge facilities. The input motion was applied as a horizontal acceleration time history to the base and sides of the model, and as a vertical acceleration time history to the bottom of the model. Time histories were provided for each experiment at two locations and so an average time history was computed. This did not allow for consideration of container rocking which appeared to be significant for some of the experiments. The average time histories were baseline corrected to remove any residual displacement drifts. A quadratic baseline correction was applied to all the motions of all facilities by fitting a quadratic curve with the least squares method and subtracting from the recordings.

Solution Scheme
The numerical platform utilized (FLAC) solves the full dynamic equations of motion for each zone and follows an explicit integration scheme where the equations of motion are used to derive new velocities and displacements from stresses and forces. Strain rates are then derived from velocities, and new stresses from strain rates. This calculation loop is performed at each time increment or time step. The time step is calculated to be small enough that information cannot physically pass from one element to another in that interval, and this way, the computational information is always ahead from the physical information. No iteration process is necessary when computing stresses from strains in a zone, even if the constitutive law is highly nonlinear (as is the case with PM4Sand). The disadvantage of the explicit method is that the required dynamic time step is very small resulting in significant computational time for a simulation. Time steps for the simulations in the current study fluctuated around 7.3eÀ6 s. All simulations were conducted with "large" deformations enabled, which allowed the mesh nodes to update their coordinates during dynamic shaking, and the geometry progressively to change. Rayleigh damping was set to 0.5% at a center frequency of 1 Hz (frequency of input motion) in order to reduce noise in the simulations and to account for small strain damping which is not captured by the constitutive model.

Soil Input Parameters
The achieved densities in the centrifuge tests varied between the different tests and none exactly matched the target value (62%). The model parameters h po and G o are dependent on D R and would thus be expected to vary between the different experiments. To address this variation, four separate sets of simulations were performed with each one taking a slightly different approach in the selection of h po and G o . All other model parameters retained the values from Table 24.1. The first two sets of simulations were considered as the "best estimate" approaches and are discussed in detail below. The third set of simulations examined the effect of using the properties in Table 24.1 without adjusting for any relative density variation and essentially assuming the target density was achieved, while the fourth set of simulations examined the effect of using the CPT data provided by the organizers to estimate the D R . Only the results from the first two sets of simulations (referred to as Calibration 1 and 2) are presented here, as these represented the authors' best estimate of the soil behavior.
Calibration 1: For this calibration, the cyclic resistance ratio (CRR) of the soil was estimated using the relationship shown in Fig. 24.2a. G o was estimated using the calibrated relationship shown in Fig. 24.2b. h po was then estimated using singleelement simulations following the procedure described above.
Calibration 2: This calibration was performed in an identical manner to Calibration 1, except for G o which was calculated using the default PM4Sand relationship shown in Fig. 24.2b. This calibration was performed to examine the sensitivity of the simulation results to G o which, as noted, was poorly constrained by the experimental data.

Simulation Procedure
The simulations were performed using the reported centrifugal acceleration from each facility to establish the pre-shaking stress conditions. The sequence of centrifuge model construction was simulated to establish a reasonable initial stress state for the soil. The soil was sequentially placed (or constructed) in layers under a gravity of 1/N g, which simulates the stresses in the model before centrifugal spinning, where N is the centrifugal acceleration from the experiment. During construction, the sand was assigned a Mohr-Coulomb material model with appropriate stress-dependent stiffness. A K o value close to 0.45 was expected for normally consolidated sand, so the Poisson's ratio for the Mohr-Coulomb model was chosen to produce this value under one-dimensional conditions. The sand was placed layerby-layer to avoid any stress arching at the container walls. Saturation was performed by setting the saturation value to 1 everywhere and allowing pore pressures to reach equilibrium. Next, gravity was increased from 1/N g to 1 g in 10 increments to simulate the spin-up of the centrifuge. Following each increase in gravity, the stressdependent stiffness moduli were updated, and the coefficient of earth pressure at rest (K o ) was evaluated. The slight slope of the model produced K o values that differed only slightly from 0.45. Since the modeling took place for the prototype scale, the shape of the ground surface (which during shaking depends on the direction of shaking relative to the curved g-field) was not curved.
After establishing the final static stresses, the water table was applied at an elevation of 5 m. The fluid boundary conditions were set as described in the Boundary Conditions section of this paper. The hydraulic conductivity of the sand was set at 1.18e-4 m/s based on the results presented by Vasko et al. (2014). The centrifuge experiments all used an appropriate viscous fluid, so the laboratory permeability was used without scaling. Example stress profiles at this stage are shown in Figs. 24.5,24.6,24.7,and 24.8. After initial conditions had been established, the PM4Sand was assigned to all zones using the selected properties. The entire mesh was simulated with the same material model since it was all made of uniform sand with a uniform permeability.
The dynamic shaking was applied to the model using the procedures described in the Boundary Conditions section of this paper. The shaking was continued for the duration of the recorded event for each facility (approximately 26 s). After the end of shaking, the motion of the model base was stopped by applying a vertical and horizontal acceleration opposite in direction to the average velocities along this boundary and with a magnitude to reduce the average velocity to zero in 0.5 s. This is equivalent to the centrifuge container being stopped, but the ramping procedure described above was used to avoid sudden changes in velocity. After the velocity of the model base was brought to zero, the model was allowed to reconsolidate. During reconsolidation, the permeability of the sand was increased by 10 times to account for an increase in permeability of the sand layer which has been observed by others during liquefaction (e.g., Shahir et al. 2012Shahir et al. , 2014. During reconsolidation, the Post_Shake flag of PM4Sand was activated (set equal to 1) with G sed ¼ 0.03 and p sedo ¼ À17,000 (Ziotopoulou and Boulanger 2013).

Comparison to Experimental Results
Simulations were performed for each experiment using the two calibrations described above (Table 24.2). Results were compared in terms of horizontal accelerations, excess pore pressures, and vertical and horizontal displacements at the locations shown in Fig. 24.3. The results generally showed that the two calibrations enveloped the observed responses for many of the tests in terms of excess pore pressures and accelerations. The results also generally overpredicted displacements at the center of the model. Excess pore pressures are compared between the numerical and experimental results for eight of the experiments in Figs. 24.9 and 24.10. Pore pressures are shown for piezometer PP2 which is located near the center of the model (Fig. 24.4). Figure 24.9 shows the time histories for the duration of shaking and shows that both calibrations were able to match the rate and magnitude of pore pressure generation for most tests. Some discrepancies are noted between the maximum excess pore pressure between the simulation and experiments which is attributed to differences between the actual location of instruments and the assumed location in the simulations. Figure 24.10 illustrates an extended version of these time histories to show the dissipation process. Calibration 2 does a better job of matching the dissipation pattern for the majority of the simulations.  Horizontal acceleration response spectra are compared for the numerical and experimental results for eight of the experiments in Fig. 24.11. The response spectra are calculated based on the horizontal accelerations at accelerometer AH3 which is located near the center of the model (Fig. 24.4). The two calibrations produced Fig. 24.9 Comparison between simulation and experimental results for excess pore pressure generation at piezometer PP2 using Calibrations 1 and 2 Fig. 24.10 Comparison between simulation and experimental results for excess pore pressure dissipation at piezometer PP2 using Calibrations 1 and 2 similar results in terms of acceleration spectra. Overall, simulations were able to reproduce the frequency content, but the magnitudes were over predicted in some simulations and under predicted in others, making the extraction of a discernable trend challenging. These discrepancies would be most likely attributed to experimental variation and uncertainty, since the same numerical modeling protocol was used in all simulations and any achieved and reported properties from the centrifuge facilities were honored. It is expected however, that a case-by-case calibration and simulation would resolve these issues.
A markedly lesser degree of agreement was observed in the comparison between simulations and experimental results in terms of deformations (i.e., lateral displacements). The simulations tended to overpredict displacements for almost all of the experiments. Calibration 2 tends to show lower displacements, but still tends to overpredict the observations. An example of the horizontal displacement time histories for the center of the profile (location shown in Fig. 24.3) is shown in Fig. 24.12. Reasons for the discrepancies in the displacements became apparent after the release of the experimental data which included estimates of D R values for each centrifuge test that had previously not been provided. These estimates were based on e max and e min values that were selected by the LEAP organizers and differed from the values selected by the authors during the calibration phase. This difference led to a systematic underprediction of the D R for most of the experiments and an associated overprediction of the displacements, which are very sensitive to small changes in D R (Fig. 24.13). Fig. 24.11 Comparison between simulation and experimental results for acceleration response spectra at accelerometer AH3 using Calibrations 1 and 2

Sensitivity Study
A sensitivity study was performed at the request of the LEAP organizers to identify how the simulation results varied with changes in density and input motion. This sensitivity study included five different input motions, which had different magnitudes (as represented by peak ground acceleration) and different frequency contents, and three dry densities. PM4Sand requires D R values, so these were calculated using the two sets of e max and e min values to account for the uncertainty in these parameters discussed previously. This led to simulations being performed at six different D R values (Fig. 24.13). The sensitivity study was only performed using Calibration 1 (Table 24.2). Lateral displacements for each of the simulations versus D R are shown in Fig. 24.13. The displacements in this figure represent final horizontal surface displacements at the center of the model (Fig. 24.4). The results show that the amount of displacement is very sensitive to changes in D R and peak ground acceleration (PGA). Displacements are especially sensitive to changes in D R when values are between 40 and 60% which is where most of the experiments for this exercise fall. These sensitivity results demonstrate that a systematic underprediction in D R could lead to an overprediction in displacement as was observed in the previous comparisons. It also emphasizes the need for reliable estimates of e max and e min when interpreting void ratios from experimental data.

Summary
Results have been presented from a set of numerical simulations using the constitutive model PM4Sand and the numerical platform FLAC. PM4Sand was calibrated to results from cyclic triaxial tests on Ottawa F-65 sand. The calibrations were performed for three D R values and relationships were developed between the input parameters and D R to allow for simulations to be performed at any D R value. The calibrated model was used to simulate nine centrifuge experiments that examined the response of a gently sloping saturated sand deposit. The simulations were able to envelope most of the recorded responses, but some discrepancies were observed. Recordings of excess pore pressure and accelerations generally matched well with results falling in the range of experimental uncertainty. The simulations tended to overpredict lateral displacements which is likely due to a systematic underestimation of D R for the centrifuge experiments. This underestimation was caused by the choice of e max and e min during the model calibration process. A sensitivity study was performed to investigate the effects of D R and input motion on the displacement results. The sensitivity study showed that displacements are very sensitive to small changes in D R . This result suggests that a more accurate assessment of the index parameters for the soil could have led to better agreement between the simulations and experiments. For future simulations using PM4Sand and FLAC, careful measurement of D R should be performed for the experiments to avoid the need for the numerical modeler to make assumptions regarding the experimental results. Overall, the simulation results suggest that the combination of FLAC with PM4Sand can give reliable predictions of the different response metrics with proper calibration.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.