eScience: Nonlinear Soil Effects in Large-Scale 3D Ground Motion Simulations

Ricardo Taborda, Julio López, Jacobo Bielak

We have implemented a rate-dependent plasticity approach to model the plastic behavior of soils in large-scale earthquake simulations using Hercules, the octree-based finite-element earthquake simulator developed by the Quake Group at Carnegie Mellon University. In this first study we have modeled the soil using the von Mises and Drucker-Prager constitutive equations. We test this new implementation in a realistic basin filled with three layers of soft-soil deposits with shear wave velocities (V_s) 200, 350, and 650 m/s, resting on a bedrock of V_s = 2600 m/s. The material within the basin is allowed to deform plastically while the bedrock substrata remains elastic throughout the simulation. Results indicate that there exist three-dimensional (3D) effects that cannot be modeled using approximate two- or one-dimensional approaches. Typical observations from past earthquakes such as permanent deformations, change in the distribution of energy in the frequency domain (reductions in low frequencies and increase of high frequencies), deamplification factors of velocity and accelerations, are all present. In addition, we observe 3D basin, edge, and directivity effects that do not necessarily follow the same patters of the corresponding linear anelastic simulation. We believe that including the nonlinearity effects of soils in large earthquake ground motion simulations is a key ingredient toward the advancement of computational seismology as a tool for improving seismic hazard analysis of earthquake-prone regions. The present implementation and results are a step forward in that direction.

Inelastic Wave Propagation

As in elastic wave propagation problems, we start from the linear momentum equation shown in (1) using indicial notation. σ_ij is the Cauchy tensor of stresses, and f_i and u_i are the body forces and displacements in the i direction. Dots stand for time derivatives and subscripts following a comma represent partial derivatives in space.

Applying finite elements in space using the standard Galerkin method, and keeping the formulation in terms of both displacements and stresses, we obtain the system of ordinary differential equations shown in (2). Here, M is the global mass matrix of the system and f is the assembled vector of nodal forces. The sum over e represents assembling of the finite elements representing the domain. B is the strain matrix. The integral over the element's volume Ω_e is computed numerically using the strains and stresses at the elements' quadrature points. For simplicity, we have omitted terms associated with boundary conditions and intrinsic attenuation.

The stresses in (2) are obtained using Hooke's law of elasticity, by expressing the elastic strain tensor as the difference between the total strain corresponding to the current state of displacements, ε_ij, and the strain due to plastic deformation, , as seen in (3).

Using central differences to express the second derivative of displacements in (2), and decoupling the system with respect to M, the forward solution of displacements for any given node i in the mesh is given by (4), where n represents a given step at time t = nΔt.

Rate-Dependent Plasticity

Since the total strain (ε_ij) at any step n can be computed from the displacements, u_n, it follows from (3) that the solution of the inelastic wave propagation problem depends on the computation of the plastic strains () at each time step n. Onceat n is known, one can proceed to the solution of (4).

In most geotechnical engineering applications where plasticity is properly considered, soils are modeled as continuum rate-independent materials. Although a rate-independent approach is attractive because of its intuitively straight forward formulation, it requires of an iterative procedure to correctly predict at step n from step n-1. This may be computationally intensive, which makes it undesirable for large-scale earthquake ground motion simulations. Therefore, we have opted for a rate-dependent formulation, which is completely step-by-step explicit.

In rate-dependent plasticity, the change in the plastic strain, or plastic strain rate () is known, and can be expressed explicitly as a function of the current state of stresses, as computed based on a particular constitutive model (F, g) and a set of material properties ().

Using a power low of the form proposed by Perzyna (1963, 1966), can be written as shown in (5), where is the material strain rate and m is the strain rate sensitivity.

A great advantage from using rate-dependent plasticity is the fact that as m → 0, the solution of the rate-dependent formulation approaches that of the rate-independent case. In fact, it can be shown that the rate-independent formulation is a special case of rate-dependent theory (Perzyna, 1966).

Soil Constitutive Models

We consider two basic constitutive models for the soil, the von Mises and Drucker-Prager yield criteria (Fig. 1). Although these are among the simplest models in elastoplasticity, they provide a valid first approximation to study the effect of soil nonlinearity at a regional scale. In particular, they allow us to reproduce phenomena that cannot be modeled with other commonly used linear equivalent methods (e.g., permanent deformations). In the future, we are interested in incorporating more complex models that will include other features such as non-associative flow, cyclic mobility, non-constant hardening, and pore-pressure for undrained materials.

Figure 1. Yield surfaces for the von Mises (left) and Drucker-Prager (right) constitutive models.

Case Study: Euroseistest

We build upon the model used for the Euroseistest Verification and Validation Project (Bard et al., 2008; Pitilakis, 2008). The simulation domain covers an area of 16 km x 29 km and has a depth of 41 km that includes a portion of the Mygdonian basin between the Lagada and Volvi lakes near Thessaloniki in northern Greece (Figs. 2 and 3). The model of the basin is composed of three main deposits with V_s = 200, 350, and 650 m/s resting on a bedrock with V_s = 2600 m/s that increases with depth up to 4400 m/s at 41 km. We simulated an earthquake scenario of magnitude Mw 5.2, with hypocenter 5 km right beneath the valley (Fig. 4). The simulation maximum frequency was 4 Hz. Only the soil deposits in the basin were allowed to deform plastically.

Figure 2. Region of interest and horizontal projection of the simulation domain.

Figure 3. Material model in the basin region with respect to the simulation domain surface. Isosurfaces are fictitiously shifted with respect to each other.

Figure 4. Simulation domain, and location and main characteristics of the source.

Results and Analysis

Deamplification Factors and Permanent Deformations

We considered two nonlinear cases with moderate and severe levels of plasticity. The yield limits of the soil deposits in the basin were set artificially based on the response of the linear elastic case. We are interested in two effects that have been long regarded as evidence of nonlinear soil behavior during earthquakes: deamplification factors with respect to the linear elastic case (Fig. 5), and the occurrence of permanent deformations away from the source dislocation due to plastic deformations (Fig. 6).

We observe large (> 2) deamplification factors in velocity and acceleration for the moderate and severe nonlinear cases, and amplification factors of up to 1.4 in the surface displacements above the deepest region of the basin for the severe nonlinear case.

In the comparison of the final displacements in the two horizontal components of motion we observe that the severe case exhibits clear directivity effects, whereas the moderate nonlinear case shows a higher spatial variability of permanent displacements that seems to combine both directivity and 3D basin effects. The comparison with respect to the linear elastic case shows a difference of an order of magnitude in the final permanent displacements.

Figure 5. Deamplification factors of the moderate (top) and severe (bottom) nonlinear cases with respect to the linear simulation for the horizontal peak magnitude of displacements (left), velocities (center), and accelerations (right) at the free surface. In parenthesis at the top-left corner of each frame are the minimum and maximum values of each case.

Figure 6. Comparison between the linear and nonlinear simulations for the final permanent displacement in the two horizontal components of motion.

Synthetics in Time and Frequency

When looking at the ground response at particular locations we observe that nonlinearity considerably attenuates the ground response, even for the moderate plastic case. Figure 7 shows a comparison of synthetics for a station on the surface near the epicenter. The severe case shows large permanent deformations in both directions of motion.

In frequency, we observe that the energy content drops for low frequencies (f < 3 Hz) and increases for higher values, for both nonlinear cases as compared to the elastic response. This is shown in Fig. 8 for the case of a station west from the epicenter. This behavior has been observed in past earthquakes and is regarded as evidence of nonlinear soil effects.

Figure 7. Comparison of time synthetics for the linear (blue), moderate nonlinear (red), and severe nonlinear (green) cases at a station S5 (see Fig. 6).

Figure 8. Comparison of normalized Fourier amplitude of velocities for the linear (blue), moderate (red), and severe (green) cases in S2 (see Fig. 6).

Evolution of Stress-Strain Relationships

Finally, we look at the shear stress-strain relationships at particular points in a downhole-like array (Fig. 9). Both nonlinear cases present single and double sided hysteretic curves, some of which exhibit clear 3D effects (coupling between the different components, xy, xz, and yz). This is a behavior that can only be simulated using 3D models and a full 3D nonlinear formulation.

Figure 9. Stress-strain relationships as implicit functions of time for the downhole observation points beneath station S2 (Fig 6). Blue lines correspond to the linear case, red and green lines correspond to the moderate and nonlinear cases, respectively. Solid dots indicate the final state of deformation.

Conclusions

We implement a rate-dependent formulation to incorporate nonlinear soil behavior in large-scale earthquake ground motion simulations. Our first results using simple constitutive models of the soil show that material nonlinearity reduces the final response of the ground, exhibiting basin and directivity effects different from those observed under linear elastic assumptions. Soil nonlinearity also causes permanent deformations, increases the spatial variabilityof the response, and changes the energy distribution in the frequency domain. At the particle level we observe 3D coupling between the different stress/strain tensor components. These are important effects that cannot be simulated using linear equivalent approximate solutions or 1D and 2D models. Incorporation of full 3D nonlinear soil behavior is an important component to be considered in future large-scale ground motion simulations.

REFERENCES

[1] Bard P.-Y., Chaljub E., Hollender F., Manakou M. and Pitilakis K. (2008), "Euroseistest Numerical Benchmark: Phase I", Commissariat à l'énergie atomique (CEC), Report distributed to the participants of the Cashima-Euroseistest Benchmark.

[2] Perzyna, P. (1963). The constitutive equations for rate sensitive plastic materials. Quarterly of Applied Mathematics, 20:321-332.

[3] Perzyna, P. (1966). Fundamental problems in viscoplasticity. In Chernyi, G., Dryden, H., Germain, P., Howarth, L., Olszak, W., Prager, W., Probstein, R., and Ziegler, H., editors, Advances in Applied Mechanics, volume 9, pages 243-377. Elsevier.

[4] Pitilakis K. (2008), "EuroSeisTest report for the Cashima project", Technical report, Department of Civil Engineering, Aristotle University of Thessaloniki, Greece, Report distributed to the participants of the Cashima-Euroseistest Benchmark.

People

LEAD RESEARCHER
Ricardo Taborda (post-doc CEE CMU)

FACULTY
Jacobo Bielak (CEE CMU)
Julio López

Acknowledgements

This research is funded in part by the Southern California Earthquake Center (SCEC) through the NSF Cooperative Agreement EAR-0106924 and USGS Cooperative Agreement 02HQAG0008, in part by NSF-OCI award Towards Petascale Simulation of Urban Earthquake Impacts (NSF OCI-0749227), by the Gordon and Betty Moore Foundation in the eScience project, and an allocation through the TeraGrid Advanced Support Program. Computations were done using Kraken at the National Institute for Computational Sciences.

We thank the members and companies of the PDL Consortium: Amazon, Datadog, Google, Honda, Intel Corporation, IBM, Jane Street, Meta, Microsoft Research, Oracle Corporation, Pure Storage, Salesforce, Samsung Semiconductor Inc., Two Sigma, and Western Digital for their interest, insights, feedback, and support.