**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

*State of the Art and Practice in the Assessment of Earthquake-Induced Soil Liquefaction and Its Consequences*. Washington, DC: The National Academies Press. doi: 10.17226/23474.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

*State of the Art and Practice in the Assessment of Earthquake-Induced Soil Liquefaction and Its Consequences*. Washington, DC: The National Academies Press. doi: 10.17226/23474.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

*State of the Art and Practice in the Assessment of Earthquake-Induced Soil Liquefaction and Its Consequences*. Washington, DC: The National Academies Press. doi: 10.17226/23474.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

**Suggested Citation:**"8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences." National Academies of Sciences, Engineering, and Medicine. 2016.

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

8 Use of Computational Mechanics to Predict Liquefaction and Its Consequences Key Findings and Conclusions ï· Computational methods that use principles of mechanics and incorporate appropriate constitutive relations offer a more detailed and flexible approach to predicting liquefaction and its consequences than do empirical and semiempirical methods. ï· Use of computational methods in engineering practice is hampered by uncertainties in their accuracy and their ability to capture important phenomena, difficulties evaluating necessary constitutive model parameters, and other complexities of their application. ï· Rigorous and innovative experimentation and analyses of the behavior of liquefying soil, along with laboratory and field validation of the results, could raise constitutive relations to a level commensurate with computational abilities. ï· Computational models need to be validated through comparison of calculations with previously vetted, accurate, and precise physical models (e.g., centrifuge test results), well-documented field case histories, or previously validated benchmark analyses that engage the important features of problems under consideration. ï· No current computational approach or computer code can address all aspects of the behavior of liquefiable soil. Computer codes vary in their constitutive models for materials, special features (e.g., the availability of structural elements), and idiosyncrasies during use. ï· Research is needed on ways to implement advanced constitutive models and on incorporating mixed formulations to solve simultaneously for deformations and pore pressures. o Instability and flow problems characteristic of liquefied soil could be addressed with meshless and discrete element methods used in conjunction with multiscale approaches and depth averaged flow models. o The smoothed particle hydrodynamics (SPH) method, the material point method (MPM), and the discrete element method (DEM) hold promise for simulating liquefaction flow and its consequences. The study committeeâs statement of task (see Box 1.4) directs the committee to explore the sufficiency, quality, and uncertainties of various methods used to assess the potential for and consequences of liquefaction triggering. Previous chapters describe the most common PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW 170

COMPUTATIONAL MECHANICS TO PREDICT CONSEQUENCES 171 approaches applied in research and practice today. Those empirical and semiempirical approaches, however, cannot account for all effects that arise from site-specific geologic details, structural configurations, and ground motion characteristics. As a result, the committee explored methods used more commonly for other engineering applications: for example, in the modeling of landslides. Engineering mechanics offers a framework for fundamentally sound approaches to modeling liquefaction and capturing its effects. Computational methods that use the principles of mechanics and incorporate appropriate constitutive relations1 applied to a properly characterized site provide a means to capture such behaviors as pore-pressure generation, redistribution, and dissipation; void (porosity) redistribution; soil and porewater flow; and soil-structure interaction effects. These methods, by necessity, rely on constitutive relations that describe the response of a soil to an external load (either force or displacement). When combined with appropriate generalizations of fluid flow relative to the solid phase, computational models can be used to solve boundary-value problems such as the deformation and pore-pressure response of a soil layer subjected to earthquake loading. Whether for the purpose of establishing boundary conditions or validating modeling results, site characterization and ground motion characterization conducted at a level of detail appropriate to analysis performed through numerical modeling will affect the fidelity of numerical modeling. Numerical models cannot be trustworthy representations of reality without an understanding of the physical site and its ground motion response. The elements of site characterization for case history assessments described in Chapter 3 are also appropriate to inform numerical modeling. Engineering-mechanics-based computational methods are accepted and used widely in structural, mechanical, materials, and manufacturing engineering, in the modern seismological prediction of earthquake ground motions, and in many geotechnical engineering applications not involving liquefaction. Computational methods are also increasingly applied to liquefaction analyses on projects where liquefaction-induced failure could threaten life safety or have substantial economic or environmental impacts. However, the accuracy of those liquefaction analyses is uncertain, important phenomena may not be captured by available models, and a relatively high level of user awareness of the complexities, idiosyncrasies, and limitations of a particular computational model is required for their proper application. Computational mechanics models have been used to investigate consequences of liquefaction in the small deformation or stable regime for a variety of geotechnical problems, including the effect of liquefaction on ground motions (Xu et al., 2003; Peng et al., 2004; Lu et al., 2004; Gyori et al., 2011; Kramer et al., 2011); lateral spreading (Elgamal and Yang, 2000; Elgamal et el., 2002; Yang and Elgamal, 2002; Peng et al., 2004; Lu et al., 2004; Forcellini, et al., 2013); settlement (Elgamal et al., 2005; Shahir and Pak, 2010; Gyori et al., 2011); and the effect of liquefaction on structures such as foundations, lifelines, dams, and retaining structures (Manzari, 1996; Elgamal et al., 2005; Marcuson et al., 2007; Shahir and Pak, 2010). Use of computational models in cases where large deformations and potential instabilities are associated with liquefaction is much less common and reliable than when deformations are small and the overall system remains stable. This may be attributed to the fact that constitutive relations often do not accurately capture the behavior of soil as it approaches liquefaction and enters the post-triggering regime. Many available constitutive relations begin to deviate 1 A constitutive relation is a set of mathematical equations that describe the response of a specific material or element to an external forcing load (or stress) or displacement (or strain). PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

172 STATE OF THE ART AND PRACTICE IN THE ASSESSMENT OF EARTHQUAKE- INDUCED SOIL LIQUEFACTION AND ITS CONSEQUENCES substantially from some aspects of stress-strain behavior as soil approaches the liquefaction state (i.e., at zero or very small effective stress), and they were never intended to be applied in the post-triggering regime. While many engineering applications are not concerned with large deformations and the post-triggering regime, new models capable of capturing post-triggering features such as phase transitions (i.e., from solid-like to fluid-like behavior), fabric effects, and rate dependence need to be developed and validated to provide a complete picture of liquefaction effects. COMPUTATIONAL LIQUEFACTION MODELING IN ENGINEERING PRACTICE Geotechnical engineers now use a variety of empirical and approximate computational models and codes to analyze liquefaction problems. The simplest are one-dimensional codes that rely on effective stress formulations of soil stress-strain/pore-pressure generation under cyclic loading, and on solving the one-dimensional wave equation to simulate nonlinear site response and the associated pore-pressure evolution due to seismic loading. D-MOD2000 and Deepsoil (Matasovic and Hashash 2012) are typical examples of software available to those in practice. Soil stiffness and strength evolution due to pore-pressure generation and dissipation are accounted for via backbone curves2 that are part of the constitutive model used in the code. The codes typically use pore-pressure models that soften the soil monotonically. Because these models ignore dilation that occurs often when the shear strain increases, they can over-soften the liquefied soil and lead to excessive displacements and unrealistic âbase isolationâ effects: that is, soil layers above the liquefied layer are isolated from the strong ground motions (Anderson et al., 2011). Some one- dimensional codesâfor instance, NOAH (Bonilla, 2001; Bonilla et al., 2005) and PSNL (S. Kramer, unpublished dataâaccount for phase transformation behavior of soils and they are, therefore, more appropriate for capturing the behavior of soil that dilates after liquefaction triggering. Those models are generally not available commercially, however. In addition, two- dimensional codes such as FLACTM (Fast Lagrangian Analysis of Continua), 3 FLIP (Finite Element Analysis of Liquefaction Program), 4 and Open System for Earthquake Engineering Simulation (OpenSees;5 described later in this section) have been used successfully to run one- dimensional analyses if the soils are constrained to behave as shear beams. More advanced codes rely on two- and three-dimensional solutions to the coupled system of equations that describe soil response in boundary-value problems. In engineering practice, one of the most widely used two-dimensional programs for modeling soil liquefaction and its consequences is FLACTM. FLACTM is a finite difference code that allows implementation of user-defined constitutive models, includes structural elements and energy absorbing boundaries, and can employ a Lagrangian large strain formulation to accommodate large displacements. One user-developed constitutive model used widely to model sandy soil liquefaction and its consequences in FLACTM is UBCSAND (Beaty and Byrne, 2011). Figure 8.1 shows a post- liquefaction snapshot of FLACTM-simulated displacements in the Upper San Fernando Dam from 2 A backbone curve connects the tips of the hysteresis loops generated by uniform cyclic loading. In the absence of cyclic degradation, it is analogous to the monotonic loading stress-strain curve for the soil. 3 See http://www.itascacg.com/software/flac. 4 See http://flip.or.jp/flip_english.html. 5 See http://opensees.berkeley.edu. PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

COMPUTATIONAL MECHANICS TO PREDICT CONSEQUENCES 173 the 1971 San Fernando earthquake. Note that the liquefaction-induced deformation of the Upper San Fernando Dam, with a maximum value of approximately 3 meters, was significantly less than that associated with the Lower San Fernando Dam (described in Chapter 1). While the parameters of the UBCSAND constitutive model can be developed from laboratory test data, Beaty and Byrne (2011) provide parameters that allow UBCSAND model results to be consistent with the deterministic standard penetration test (SPT) and cone penetration test (CPT) liquefaction triggering curves for clean sand (described in Chapter 4) for uniform cyclic loading. Boulanger and Ziotoupolou (2012) implemented a modified form of a bounding surface constitutive model (Manzari and Dafalias, 1997) called PM4Sand into FLACTM. Boulanger and Ziotoupolou describe a methodical calibration process that allows their model parameters to be developed from conventional soil property data, facilitating use of the model in engineering practice. Both UBCSAND and PM4Sand were developed to model cyclic mobility and its effect on triggering and post-triggering deformations; to the committeeâs knowledge, their ability to capture post-triggering flow deformations is still under debate, although PM4Sand has a critical state framework and would seem to have no fundamental limitation against simulating flow deformation. FLACTM, a finite difference approach, is sometimes considered on the opposite end of the computational spectrum from finite element approaches such as that employed in PlaxisTM, 6 another commercial computer program widely used among engineering firms that do computational-mechanics-based analyses. PlaxisTM has many of the same features as FLACTM does, including structural elements, energy-absorbing boundaries, and user-defined constitutive models (e.g., including UBCSAND). FLACTM uses both simple but dense finite difference meshes and a solution scheme that treats all problems as dynamic problems. PlaxisTM uses more elaborate elements (e.g., 15-noded triangles) to achieve a high degree of resolution. A significant difference between these two programs is the method used to input an acceleration-time history to the analysis. PlaxisTM uses the conventional approach of inputting an acceleration-time history to the base of the soil layer. FLACTM requires input of a velocity-time history that represents only the incoming seismic wave. Mejia and Dawson (2006) discuss a proper method of inputting an earthquake time history to FLACTM. Discussion of their relative merits is beyond the scope of this report, but these programs are well respected and widely used in the earthquake engineering community. Matasovic and Hashash (2012) provide a brief general discussion of their use for site response analysis considering pore-pressure generation as well as the need for calibration and validation of such analyses. FIGURE 8.1 Displaced shape and displacement vectors at the end of post-earthquake analysis of Upper San Fernando Dam in the 1971 San Fernando Earthquake. This figure shows the ability of a common commercial 7 program to capture relatively large displacements. SOURCE: UBCSAND Constitutive Model on Itasca UDM website. 6 See http://www.plaxis.nl. 7 http://www.itasca-udm.com/media/download/UBCSAND/UBCSAND_UDM_Documentation.pdf. PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

174 STATE OF THE ART AND PRACTICE IN THE ASSESSMENT OF EARTHQUAKE- INDUCED SOIL LIQUEFACTION AND ITS CONSEQUENCES OpenSees is an open source, finite element-based code for calculating the performance of structures and geotechnical systems subjected to earthquakes. The code has an array of sophisticated structural, geotechnical, and interface elements and models that allow it to analyze a variety of soil-foundation-structure interaction (SFSI) problems, including those involving liquefiable soils. A number of constitutive models capable of handling liquefiable soils, including bounding surface (e.g., Manzari and Dafalias, 1997) and multi-yield surface family models (Elgamal et al., 2003; Yang et al., 2003) are available under this numerical platform. There is a steep learning curve associated with the use of OpenSees, however, and because it is an open source program, limited technical support is available. Hence, OpenSees is used much more commonly in research than in consulting engineering practice at this time. ISSUES IN THE COMPUTATIONAL MODELING OF LIQUEFACTION PROBLEMS Computational methods to analyze liquefaction problems solve the equations of motion and deformation for assemblages of solid particles containing fully or partially fluid-filled pore spaces. Mechanical principles require that stresses and displacements of the solid particle framework, and fluid pressures and flows relative to those particles, satisfy the equations of motion and of conservation of mass. Both the deformations of the solid particle framework and the migration of porewater need to be accounted for, as does the transfer of momentum between the solids and liquids. The nonlinear elasto-plastic nature of soil, the added complications of pore fluid pressurization and flow, and the large deformations associated frequently with liquefaction make problems in this domain among the most challenging in modern computational engineering. To simplify the analysis for liquefaction problems, it is generally assumed that the soil is saturated, that the solid constituents (i.e., soil particles) are incompressible8 (although the soil skeleton may be compressible), and that the motion of the interstitial fluid with respect to the soil particle matrix is governed by diffusion equations. These assumptions are similar to the conventional assumptions of Terzaghi consolidation theory, as extended by Biot (1941) to three- dimensional deformation fields. Prevost (1980, 1985b), Prevost and colleagues (1985), Lacy and Prevost (1987), and Zienkiewicz and colleagues (1990a,b, 1999) used continuum mixture theory concepts to generalize Biotâs (1956a,b) dynamic poroelastic formulation into a form suitable for advanced computational modeling. In these approaches, only the equations of overall balance of linear momentum and of mass conservation need to be resolved (Biot, 1941; Prevost, 1980; Zienkiewicz and Shiomi, 1984; Andrade and Borja, 2007). The equations are discussed in Box 8.1. When combined, these sets of equations furnish a coupled system within which the deformation of the solid matrix and the pressure and transport of the pore fluid can be solved. Complexity here stems both from the coupling between porewater pressures and solid matrix deformations and from the intrinsic complexity of soil behavior. Setting the boundary and initial conditions is required for an analytical or a numerical solution. Because of the complexity of the problem, analytical solutions 8 Some formulations use the bulk modulus of the fluid to approximate the pore-pressure evolution. But, in the fully coupled analysis shown in Box 8.1, pore-pressure evolution is computed directly from the balance laws (Equations 1 and 2), without specifying the bulk modulus of the fluid. PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

COMPUTATIONAL MECHANICS TO PREDICT CONSEQUENCES 175 rarely are available. As a result, computational numerical solutions are typically employed in practice. In developing computational solutions for these problems, the analyst must make several choices among alternative approaches, and he or she must consider the trade-offs between analytical accuracy and computational efficiency. Some of these issues include how to discretize the problem in space and time; whether to use small or large strain formulations; whether to solve the coupled equations directly or to decouple them; whether to use a finite element or finite difference model; how to deal with large deformations and discontinuities; how to address questions of stability and flow sliding; how to model the interaction between the soil and any pertinent structures; and, finally, how to validate the computational model. At this time, there is no single computer code that addresses all of these issues satisfactorily, and there is no general agreement within the geotechnical earthquake engineering community on how best to address these issues. The following sections describe briefly the current state of the art and of practice. BOX 8.1 Basic Balance Laws for Fully Saturated Soil-Fluid Analysis Equations 1 and 2 are written using continuum mechanics conventions; a repeated subscript indicates summation over the range of the subscript, and a comma indicates differentiation with respect to the variable whose subscript follows the comma: ð ðð,ð + ðð ð = ðð ð (1) ð ð ð£ ð,ð = âð ð,ð (2) In Equation 1, the balance of momentum involves the derivative (in space) of the total stress tensor (ð ðð ), the density of the soil-fluid mixture (Ï), the acceleration due to gravity (ð ð ), and the acceleration (ð ð ) of the soil skeleton (the formulation is more complicated in the mixture theory context, although some mass- weighted acceleration ð ð can always be defined). In Equation 2, the balance of mass (flow) involves the density of the pore fluid (ð ð ), the derivative (in space) of the velocity (ð£ ð ) of the solid skeleton, and the derivative (in space) of the Darcy mass flux vector (ð ð ). To close this system of equations, the effective stress relation (plus constitutive equations linking the effective stress to the deformation of the solid skeleton) and Darcyâs equation of flow are added. If the relative motion of the two phases is vigorous, a Forchheimer formulation is possible (see, e.g., Mathias et al., 2008), with the pore-pressure gradient balancing the sum of two terms, one being linear in the Darcy mass flux and the other quadratic in that flux. Discretization in Space and Time The two classic approaches used to discretize a domain in space for a numerical solution are finite differences (see, e.g., Moin, 2001) and finite elements (see, e.g., Hughes, 1987; Bathe, 1996; Zienkiewicz and Taylor, 2013). Accuracy of the solution in space increases with the diminishing size of the âmeshâ elements and with the complexity of interpolations within the elements into which the domain is discretized. At the same time, however, computational efficiency decreases. Thus, care must be taken to choose sufficiently small or complex elements to provide an adequate level of accuracy while maintaining computational efficiency. Time discretization (or integration) is typically accomplished using the Newmark family of integrators, which can include, among others, the forward Euler scheme (explicit), backward Euler scheme (implicit), and the Crank-Nicholson mixed scheme (Hughes, 1987), each with its PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

176 STATE OF THE ART AND PRACTICE IN THE ASSESSMENT OF EARTHQUAKE- INDUCED SOIL LIQUEFACTION AND ITS CONSEQUENCES own advantages and disadvantages. Explicit schemes are typically used in practice because of their simplicity, but care must be taken to choose sufficiently small time steps to provide adequate accuracy and stability. Some nonlinear software (e.g., the well-known ABAQUS system)9 includes provisions to address this limitation by varying the time step to maintain a satisfactory energy balance at the end of each time increment. It is also important to note that finer discretization will be required for similar reasons. Implicit schemes, while unconditionally stable, are less efficient computationally and more difficult to implement in a computational model. Hybrid schemes solve the fully coupled system of balance of momentum and flow using an implicit-explicit predictor-corrector scheme. Given their unconditional stability, implicit methods for time domain integration merit further development and increased use in computational modeling of liquefaction. Crank-Nicholson schemes offer a compromise but can produce spurious oscillations under certain conditions. Small and Large Strain Formulations Computational methods used to analyze liquefaction problems can be divided into methods for small strains (typically less than ~1%) and larger strains, although there is no fixed demarcation between what constitutes small and large strain. It must be remembered, however, that large strain and large displacement are not synonymous. A body can undergo large displacements while still experiencing small strains; for example, rigid body motion involves no deformation at all. But, when the body deforms significantly, a large strain formulation may be necessary. Small strain models are computationally simpler and more efficient than are large strain models; therefore, they are generally used when the assumption of small strains is adequate. Determining that adequacy, however, is somewhat subjective. From a computational mechanics perspective, the distinction between large and small strain methods is that small (i.e., infinitesimal) strain methods can be used when expressions that are linear in the first derivatives of the displacements adequately describe the deformations and rotations. In this case, the equations of equilibrium or of motion, which apply rigorously in the deformed configuration, remain adequate when written without accounting for the geometry change. Otherwise, large (i.e., finite) strain methods must be used. Most current formulations presuppose small (i.e., infinitesimal) strains even when, in practice, they are used to simulate finite strains. Coupled and Uncoupled Solutions A fully coupled solution scheme is one in which the changes in pore fluid pressures and deformations in the soil particle framework are obtained simultaneously in each step of the marching procedure. An uncoupled solution scheme is one in which, for example, the pore fluid pressure is computed at some time and then allowed to dissipate without consideration of how this affects the deformation and stress distribution in the solid constituents. Uncoupled solutions for deformation and pore-pressure development tend to be less accurate but simpler to 9 See http://www.3ds.com/products-services/simulia/products/abaqus/latest-release. PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

COMPUTATIONAL MECHANICS TO PREDICT CONSEQUENCES 177 implement, and they are common in commercial computer programs. Uncoupled systems can be combined into staggered systems in which, for example, the dissipation solution is carried forward for one or more steps, the deformations and stresses are updated, the dissipation process is carried forward for a few more steps, the deformations and stresses are updated, and so forth and so on. Staggered solutions and fully coupled solutions are used less frequently than uncoupled solutions, but they yield the most accurate results in computational analysis of liquefaction problems. They require what is referred to as a mixed, or u-p (displacementâ porewater pressure) formulation in which both displacements and porewater pressures are unknowns that must be computed. Finite Element Versus Finite Difference Approaches There is no perfect or ideal computational approach that can address all problems. Each approach has its own strengths and weaknesses. The finite element method is very good for modeling complicated geometries with different material properties in different regions. Indeed, a motivation for developing the method was the difficulty encountered when trying to distort finite difference meshes to represent the shapes of aircraft wings and fuselages (see, e.g., Gupta and Meek, 1996). On the other hand, finite difference methods usually are more efficient at propagating solutions through time. Both finite difference and finite element computer programs are used in research and practice. Both methods can accommodate complex geometries, boundary conditions, and finite deformations. Finite element-based computational methods are by far the most common form of computational solutions in engineering mechanics. Nevertheless, one of the most popular geotechnical programs used by practicing engineers (FLACTM) is based on finite differences. An important consideration in the development of finite element and finite difference programs for seismic analysis is the treatment of boundary conditions to simulate infinite domains. Simple fixed or free boundaries unrealistically reflect stress waves back into the domain of interest and create what is known as a âbathtub effect.â Using large elements to extend the geometry beyond the local domain of interest may also cause problems: the large elements can have difficulty transmitting energy at wavelengths smaller than the element dimension, thereby introducing distortions in the response at those smaller wavelengths. Care should be taken when discretizing a problem with a mesh to ensure that frequencies below the maximum frequency of interest can be propagated through those portions of the mesh that most strongly affect the response. Material softeningâas a result, for example, of pore-pressure generation in liquefaction analysisâis also known to render computations sensitive to mesh size. It has been shown, however, that this sensitivity can be ameliorated by the viscous effect produced by fluid flow (Andrade and Borja, 2007). In selecting a finite element or finite difference code for analysis of dynamic problems, the engineer should be sure that the code selected handles the boundaries correctly. Large Deformations and Discontinuities When large deformations and discontinuities associated with liquefaction occur, both finite element and finite difference formulations have to be enhanced to account for the deformations PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

178 STATE OF THE ART AND PRACTICE IN THE ASSESSMENT OF EARTHQUAKE- INDUCED SOIL LIQUEFACTION AND ITS CONSEQUENCES and discontinuities in a manner consistent with the principles of continuum mechanics (Belytschko et al., 2000; Holzapfel, 2000). Additionally, the spatial discretization techniques (e.g., finite elements or finite difference mesh) of the computational model have to cope with the large deformations resulting from the triggering of liquefaction. Computational modeling of large deformations requires not only large deformation formulations for the governing equations but also computational techniques to accommodate changes in geometry within the soil mass. Examples of such techniques include Arbitrary Lagrangian Eulerian (ALE) techniques (Ponthot and Belytschko, 1998); adaptive remeshing (Berger and Colella, 1989); or meshless methods such as the material point method, smoothed particle hydrodynamics, and others (Sulsky et al., 1994; Randles and Libersky, 1996). The finite elements themselves also have to be updated to avoid the ensuing mesh distortions or discontinuities (Ponthot and Belytschko, 1998; Borja, 2000; Dolbow et al., 2001; Li et al., 2010). Additionally, the validity of constitutive relationsâ developed typically to represent relatively small deformationsâis questionable when the deformations being modeled are large. More experiments (e.g., simple shear, triaxial, true triaxial) are needed to characterize the response of liquefied soil to facilitate development of models that account realistically for the behavior of soil after triggering of liquefaction. The continuum approach has shown promising results in capturing some of the basic features observed in the laboratory and physical models. Nevertheless, more laboratory and physical models are needed to validate computational models and thereby increase the confidence in simulation tools to address liquefaction problems accurately. Stability and Flow Sliding From a mechanics perspective, the unstable regime is one in which relatively small changes in stress conditions yield a large response in strain (Hill, 1958; Nova, 1994; Bazant and Cedolin, 2010). In the liquefaction modeling community, the concept of physical stability is limited generally to the âsmallâ strain domain. Flow slides involve large strains and are considered to be in the unstable regime. Other technical communities (e.g., those that study granular physics) refer to stable regimes as âsolidâ stages and unstable regimes as âfluidâ stages (Jaeger and Nagel, 1992). Modeling the transition from a solid to a fluid stage, or from a stable to an unstable stage, has at least two basic requirements: (1) criteria to flag the onset of the transition, and (2) computational schemes that can accommodate large strains and possible discontinuities in displacement. A circumstance of particular interest with respect to liquefaction consequences is lateral spreading, where instability (i.e., flow failure) may not have occurred but where the soil has undergone large displacements due to significant reduction in stiffness and strength (softening). Lateral spreading in the field may also occur (Kramer, 1996; Elgamal et al., 2002) when the effective stress approaches the zero-pressure region as a result of contractive material response following multiple strain cycles. This results in significant softening but not necessarily in an unstable situation. Soil-Structure Interaction Soil-structure interaction effects, including the impact of lateral spreading on pile foundations and the lateral resistance of piles in liquefiable ground, are particularly important PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

COMPUTATIONAL MECHANICS TO PREDICT CONSEQUENCES 179 consequences of liquefaction (Chapter 7). Interpretation of aerial photographs following the 1964 Niigata, Japan, earthquake (Hamada and OâRourke, 1992; OâRourke and Jeon, 2000) indicates lateral displacements of up to eight meters. While not necessarily indicative of flow-failure instability, these displacements caused significant damage to structures. Computational modeling of these interactions is complicated by the three-dimensional nature of the problems. While several of the previously described computational platforms can accommodate three-dimensional analysis, such analyses are labor and computationally intensive, and there are uncertainties associated with three-dimensional generalizations of constitutive models. Therefore, many soil- pile interaction analyses are conducted two-dimensionally, with the pile represented as a planar element having an approximate lateral stiffness. A similar approach is applied in OpenSees, which uses p-y curves for liquefiable soils (furnished by the use of âpysimpleâ elements), taking into account pore-pressure generation and soil dilation and softening. While these analyses focus on capturing the bending response of the pile, which is critical for most engineering applications, it is sometimes necessary to represent the behavior of liquefied soil flowing around the pile. Furthermore, both two- and three-dimensional analyses of soil-pile interaction require consideration of the behavior of the soil-pile interface, another source of uncertainty. The interface behavior between structural materials and both liquefied and non-liquefied soil is not well understood. Generally, interface behavior is modeled using interface elements that employ simplistic elasticâperfectly plastic constitutive models. For example, there is a need to develop interface elements that can account for void ratio redistribution and the large deformations resulting from the separation between the pile and the soil. Improved computational efficiency for three-dimensional analyses, increased accuracy and validation for three-dimensional constitutive models, and a better understanding of interface behavior and the interface elements that model that behavior are required to obtain more accurate representation of soil-structure interaction in liquefied soil. Computational Model Validation Fundamental to validating computational models is a comparison of calculations with physical models (e.g., centrifuge test results), well-documented field case histories, or benchmark analyses that have been previously vetted for accuracy and precision and that engage the important features of the problem. For example, Peng and colleagues (2004) compare centrifuge and finite element models of lateral spreading of a two-layer soil and its influence on pile foundations. They successfully replicate the behavior of the soil and its interaction with the pile under cycle loading. Boulanger and colleagues (2005) compare analysis of the restraining effect provided by piles (referred to as âpile pinning effectsâ) in an embankment subject to lateral spreading to centrifuge test results. This type of work exemplifies some of the challenges related to predicting lateral spreading and its effect on structures: (a) the ability to capture significant shear deformations as a result of material softening (1-10% shear strain); (b) the ability of constitutive models to reproduce material behavior under numerous loading and unloading cycles and as the effective stress approaches zero; and (c) the ability of the computational scheme (e.g., finite elements) to cope with significant deformations, pore-pressure generation, and even discontinuities in the displacement field in the soil itself or as a result of separation between the soil and the pile. It is important to distinguish the need to capture realistic material behavior and the resulting liquefaction effects (e.g., pore-pressure buildup) from PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

180 STATE OF THE ART AND PRACTICE IN THE ASSESSMENT OF EARTHQUAKE- INDUCED SOIL LIQUEFACTION AND ITS CONSEQUENCES heterogeneities and complexities at the field level. While some continuum models have shown promising results when compared with physical models, more efforts such as the Liquefaction Experiment and Analysis Project (LEAP) (Manzari et al., 2014) and its predecessor, Verification of Liquefaction Analysis by Centrifuge Studies (VELACS) (Arulanandan and Scott, 1993), are still needed. Currently, no computer code or analytical approach can address all aspects of the deformation problems associated with liquefaction. Engineers working on projects for which liquefaction effects can be expected to be important are well advised to use more than one analytical or computational approach and to compare the results. Independent peer review of the computational models and results should be incorporated into important projects. CONSTITUTIVE MODELING OF LIQUEFIABLE SOIL Solution procedures that contain the equations for motion and balance of mass for a coupled system (see Box 8.1) must incorporate constitutive equations that relate effective stress changes to strain increments. In addition, Darcyâs equation for fluid flow in porous media, or possibly an extension of Darcyâs equation to account for deviations from pore-scale laminar flow conditions, is used to complete the system of equations. Because of its central importance in computational modeling, some remarks regarding constitutive modeling are appropriate. Constitutive models need to simulate (a) undrained pore-pressure generation as a function of the amplitude of shear loading and the material state (effective stress and density); (b) the effect of initial shear stress (anisotropic initial loading); (c) transformation of the soil from a solid to a liquefied state and its effects on stiffness; (d) the effects of various soil fabrics and fabric degradation on stiffness; (e) stiffness and dilatancy evolution under cyclic loading; (f) the onset of instabilities (diffuse and localized); and (g) the transition to the steady state (Kramer and Elgamal, 2001). Capturing all of these features is beyond the capability of existing constitutive models. All constitutive models that attempt to describe the complexities of soil behavior need to account for both elastic and inelastic deformations of the soil matrix. It is well documented that the relations between stresses and strains in soils are nonlinear and inelastic and display asymmetric response under compression versus extension (Schofield and Wroth, 1968; Muir Wood, 1990). Perhaps the most widely used framework to account for this complex effective stress behavior is elasto-plasticity, which accounts for both the recoverable (elastic) and the irrecoverable (plastic) deformations of the soil matrix (Desai and Siriwardane, 1984; Prevost, 1985a; Lubliner, 1990; Simo and Hughes, 1998). A number of available constitutive models represent the elasto-plastic stress-strainâpore- pressure behavior of liquefiable soils (see, e.g. Prevost, 1985a,b; Iai et al, 1990; Towhata and Ishihara, 1995; Parra, 1996; Manzari and Dafalias, 1997; Beaty and Byrne, 1998; Cubrinovski and Ishihara, 1998; Kramer and Arduino, 1999; Elgamal et al., 2000; Li et al., 2000; Bonilla et al., 2005; Ellison and Andrade, 2009; Beaty and Byrne, 2011; Boulanger and Ziotoupolou, 2012). These models typically focus on replicating the behavior of granular soils under cyclic loading. Model features generally include state parameter (pressure and density) dependence, dilatancy evolution, fabric degradation, inelasticity, anisotropy, steady-state behavior at large shear strains, and tridimensional loading. By using an elasto-plastic constitutive model and the mixed formulation for the coupling of deformation and porewater pressure described above, it PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

COMPUTATIONAL MECHANICS TO PREDICT CONSEQUENCES 181 ought to be possible to approximate the material and structural responses leading up to liquefaction and the onset of large displacements observed in the laboratory and field (see Figure 8.2). A linear constitutive model is used sometimes for the porewater pressure.10 The bulk modulus may be set equal to that of water for a saturated soil, or it may be set to a lower value if the soil is unsaturated. The interaction between the effective stress and porewater fluid pressure is governed by the coupled equations mentioned in Box 8.1. Based mostly on observations of soil behavior in uniform cyclic laboratory tests, constitutive models have evolved in their level of complexity. Some of the early constitutive models (see, e.g., Martin et al., 1975) attempted to predict the evolution of the effective stress on a cycle-by-cycle basis up to the triggering of liquefaction. (No attempt was made to predict soil behavior after triggering.) Pore-pressure generation was predicted empirically rather than on the basis of the coupling of volumetric deformation and pore pressure. Modern constitutive models use coupled models based on effective stress-strain behavior and the associated pore-pressure development in response to cyclic loads. FIGURE 8.2 Finite element liquefaction modeling and prediction in the laboratory and field. (a) Undrained cyclic triaxial test showing flow related to liquefaction for test conducted by Qadimi and Coop (2007) in loose sands; (b) Simulations of undrained triaxial test using the Manzari and Dafalias (1997) constitutive model. Similar comparisons in different stress paths have been reported in the literature (e.g., Andrade, et al., 2013); (c) field-scale simulation of flow related to liquefaction for a bridge using realistic topography, seismic excitation, and finite elements, and an appropriate constitutive response such as that shown in (b). Red colors imply larger cumulative vertical displacements. SOURCES: (a) and (b) COURTESY: J. Andrade; (c) Arduino 2014 presentation to the committee. 10 Porewater pressure is related to volume change via an appropriate constant bulk modulus assigned to the pore fluid. PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

182 STATE OF THE ART AND PRACTICE IN THE ASSESSMENT OF EARTHQUAKE- INDUCED SOIL LIQUEFACTION AND ITS CONSEQUENCES Most constitutive models based on elasto-plasticity use elastic response descriptors (e.g., elastic shear modulus, Poissonâs ratio), a yield surface, a plastic potential surface, a flow rule, and a hardening rule. This framework requires a number of parameters to fully describe the cyclic behavior of granular materials. It also relies on phenomenological observations in the laboratory to formulate the hardening rules that control the internal (plastic) variables dominating inelastic material behavior and to account for the basic mechanics such as dilatancy, anisotropy, and fabric (Desai and Siriwardane, 1984; Muir Wood, 1990). Elasto-plastic models yield better results when applied under the loading conditions for which they were calibrated (e.g., triaxial compression/extension) than under other loading conditions. This is a direct consequence of phenomenology (i.e., the constitutive equations are construed to represent observed behavior). While the ability of elasto-plastic models to capture the complex behavior of soils is remarkable, there is room to make these computational models more robust. Efforts are needed to develop constitutive relationships that reduce phenomenology in favor of more fundamental characterization of saturated granular flow. New developments in X-ray tomography and micromechanics could shed light on these efforts. These models will more likely be adopted in engineering practice when practical methods are developed to evaluate model parameters from field test data (e.g., from SPT and CPT data) in a manner similar to that presented in Boulanger and Ziotopoulou (2013). Most constitutive models are based on the observed behavior of soil in laboratory tests at relatively small strains (i.e., strains â¤ 10%); thus, they may be valid only up to liquefaction triggering. To model post-triggering behavior of liquefiable soil accurately, constitutive models that extend into the post-triggering large deformation regime are needed, but little is known about post-triggering material response. Given the instability and large deformations that often follow the onset of liquefaction, it is difficult to probe those areas of the stress space in the laboratory. Effects such as phase transition, dilatancy, fabric evolution, pore redistribution, and strain rate effects are important in capturing post-triggering behavior and are, therefore, features that developers of models and numerical approaches have focused on to obtain more accurate liquefaction and post-triggering predictions. More efforts are still needed in this direction, especially at large deformations. RECENT COMPUTATIONAL RESEARCH DEVELOPMENTS APPLICABLE TO LIQUEFACTION ANALYSIS Advanced computational techniques may be able to address such soil liquefaction issues as void ratio redistribution and large displacement flow slides following the triggering of liquefaction, and separations and discontinuities that sometimes accompany ground displacement (e.g., separations between pile and soil and ground cracks in the upstream face of the Lower San Fernando Dam; shown in Figure 8.3). Despite advances in computational capabilities, it is still difficult to use computational techniques to capture pore-pressure generation in general loading patterns beyond simple shear or triaxial conditions, to simulate large deformations, or to predict post-triggering deformations consistent with residual strength and field observations. In particular, when computational models attempt to capture the large deformations associated with liquefaction triggering, they can misrepresent flow associated with liquefaction (and associated displacements) and can fail to represent surface discontinuities associated with lateral spreading. Meshless methods, such as the MPM (Sulsky et al., 1994; Mast, 2013; Ceccato et al., 2014) and PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

COMPUTATIONAL MECHANICS TO PREDICT CONSEQUENCES 183 FIGURE 8.3 (a) Mechanisms of failure for the Lower San Fernando Dam and (b) Finite element simulation of the failure using mixed finite elements. SOURCE: Zienkiewicz, O.C., Y.M. Xie, B.A. Schrefler, A. Ledesma, and N. Bicanic. 1990b. Static and dynamic behaviour of soils: a rational approach to quantitative solutions. II. Semi-saturated problems. Proceedings of Royal Society of London A 429:311â321. By permission of the Royal Society. the SPH method (Randles and Libersky, 1996), offer great promise for accommodating large deformations through computational models based on finite deformation formulations. Other alternatives include depth-averaged flow models (Iverson and George, 2014) and discrete element methods (DEMs) (Cundall and Strack, 1979; Goren et al., 2010; OâSullivan, 2011). Meshless methods have the potential to capture complex material behaviors such as soil skeletonâpore-fluid interactionsâcrucial in modeling liquefactionâand to represent the ensuing large deformations. The continuum approach suffices to capture soil-fluid interaction using the equations of balance of momentum and mass, but the continuum approach is not applicable directly at the granular scale. At that scale, discrete element mechanics can be used to simulate the behavior of the solid soil skeleton. A computational mesh and nodes in the traditional sense are no longer necessary in meshless methods, but implementing advanced elasto-plasticity models and incorporating mixed formulations to solve the coupled dynamic and mass balance equations remain challenging. The MPM and SPH procedures (both of which are meshless methods) still require continuum level constitutive relations, and computational expense severely limits use of the DEMs. Some other difficulties associated with the DEMs will diminish as computational power increases. Despite their limitations, particle morphology and the coupling with fluid flow can be captured more faithfully with discrete element methods. Also, meshless methods have much to offer in capturing large deformations. Figure 8.4 shows simulations of flow problems using the MPM and DEM platforms, displaying the capability of these methods to represent large deformations resulting from flow associated with liquefaction or other phenomena. The aforementioned methods provide alternative numerical solution architectures that may have advantages over the traditional finite element and finite difference methods in modeling liquefaction problems and are further described in this section. PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

184 STATE OF THE ART AND PRACTICE IN THE ASSESSMENT OF EARTHQUAKE- INDUCED SOIL LIQUEFACTION AND ITS CONSEQUENCES FIGURE 8.4 Comparison of simulation of column flow problems using (a) the material point method (MPM) and (b) the discrete element method (DEM). The granular material is initially in the box and then the box is removed suddenly. Similar results have been published in the literature (e.g., Lim, et al.,2014). Both DEM and MPM simulations shown are purely gravity driven, and in the case of the DEM they account for particle features such as shape (discs versus angular) and friction. SOURCE: (a) Mast, 2013. Modeling landslide induced flow interactions with structures using the material point method.; (b) COURTESY: J. Andrade. Depth Averaging The depth-averaged method reduces both the domain size and the discretization needs and, therefore, increases computational efficiency in problems involving instability and flow. The method simulates flow behavior by solving depth-averaged versions of the balance of linear momentum and constant-volume equations. This technique, often used in the geosciences, invokes the Mohr-Coulomb failure criterion to limit the shear stress at the base of the deforming soil mass. Recent formulations also take into account the effective stress principle, observe concepts of critical-state soil mechanics, and account for evolving solid volume fractions and porewater pressures (Iverson and George, 2014). Depth-averaged equations can be solved using conventional finite differences (control volume), finite elements, or meshless techniques such as the MPM and the SPH method. The depth-averaged method has been used to simulate large-scale debris flow problems under controlled conditions and in natural settings (Iverson and George, 2014). Its two-dimensional nature and efficient time integrators make this method capable of simulating flows at the PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

COMPUTATIONAL MECHANICS TO PREDICT CONSEQUENCES 185 kilometer scale with a good level of resolution in topography. It relies fundamentally on the assumption that the depth of the flowing soil mass (the z-axis in Figure 8.5) is negligible compared with its length, making it appealing for long runout flow. When the depth is a significant fraction of the flow length, this approach has limited applicability. FIGURE 8.5 Schematic of depth-averaged flow technique for simulating debris flow by solving the coupled system of balance of linear momentum and mass equations. SOURCE: Iverson and George, 2014. Iverson, R.M., and D.L. George. 2014. A depth-averaged debris-flow model that includes the effects of evolving dilatancy. I. Physical basis. Proceedings of the Royal Society. With permission from the Royal Society. Smoothed Particle Hydrodynamics The SPH method (Randles and Libersky, 1996, and references therein) relies on computational points over which a kernel is defined. A kernel is similar to a shape function in finite elements in that it defines a family of hydrodynamic points that describe the pattern of flow. The accuracy of the method depends crucially on the number of computational points. A formulation of the SPH method can integrate the coupled system of balance of linear momentum and mass equations to obtain simulations such as the one shown in Figure 8.6 of a landslide at different time stations (Pastor et al., 2009). Fluid saturation and pore pressures are simulated, and the formulation is more efficient than are other formulations such as those based on Eulerian finite elements. The formulation is also depth integrated, making it less expensive computationally. Material Point and Discrete Element Methods The MPM (Sulsky et al., 1994) and the DEM, illustrated in Figure 8.4, are promising techniques to replicate material flow. Coupling these methods with pore-fluid flow to solve the coupled system of equations for balance of momentum and mass is a challenge. Nevertheless, these methods can be used to solve flow problems, such as the column drop test shown in Figure 8.4, with great accuracy. The column drop test simulates a granular material-filled container, the walls of which are suddenly removed. The soil particles are allowed to flow as a result of gravity. The boundary conditions for this problem impose large flow and deformations on the soil mass that cannot be captured easily with mesh-based methods. The MPM and the DEM are PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

186 STATE OF THE ART AND PRACTICE IN THE ASSESSMENT OF EARTHQUAKE- INDUCED SOIL LIQUEFACTION AND ITS CONSEQUENCES FIGURE 8.6 Simulated landslides for saturated frictional material considering pore-pressure dissipation. A dark- colored mass of material seen at the top of the slope at time 8 sec fails and flows downhill until it reaches the bottom of the slope by time 60 sec. Calculations are performed on a square mesh of about 1.5 km in each dimension. SOURCE: Pastor, M., B. Haddad, G. Sorbino, S. Cuomo, and V. Drempetic. 2009. A depth-integrated, coupled SPH model for flow-like landslides and related phenomena. International Journal for Numerical and Analytical Methods in Geomechanics 33:143â172. also capable of seamlessly simulating the interaction of soil with structures. Figure 8.7 shows an MPM simulation of landslide interactions with a set of structures. Information about the forces exerted on the structures by the flowing mass is among the computational results of this model. Both the MPM and the DEM can also easily capture discontinuities in the deformation, including complete separation of flowing masses into separate flows. FIGURE 8.7 Evolution of landslide flow into protective barrier system using the MPM. Light (blue) color represents landslide hitting a group of regularly arranged structures at time 6.25 sec. After 25 sec, the figure shows the structures being completely engulfed by the landslide. SOURCE: Mast, 2013. Modeling landslide induced flow interactions with structures using the material point method. The DEM employs three main approaches to capture soilâpore-fluid interaction on a multiscale basis. In order of increasing complexity, these approaches are (a) âdryâ DEM simulations imposing constant volume (undrained) constraints; (b) coarse-grid approaches in which the fluid flow component is simulated using a discretization much larger than the average particle size; and (c) fine-grid approaches in which the fluid flow component is simulated using a discretization finer than the average particle size (OâSullivan, 2011). The first approach is straightforward and suffices for studying elemental tests,11 where undrained conditions hold, 11 Elemental tests are laboratory tests in which uniform stress conditions are imposed on a representative sample. PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

COMPUTATIONAL MECHANICS TO PREDICT CONSEQUENCES 187 roughly, locally (i.e., porewater redistribution effects can be ignored) and the soil and fluid constituents can be assumed incompressible. On the other hand, fluid migration is relevant in liquefaction modeling at scales larger than the element test (e.g., centrifuge, field scale). In these situations, solving for fluid flow by explicitly using either coarse- or fine-grid approaches is necessary to account for hydrostatic and drag forces, thereby accounting for particle-fluid interactions. Coarse-grid approaches tend to rely either on solving Darcyâs equation (Calvetti and Nova, 2004) or the averaged Navier-Stokes equation (Tsuji et al., 1993; Zeghal and El Shamy, 2004), accompanied by the continuity equation. Specifically, Zeghal and El Shamy (2004) have used the coarse-grid approach that solves the averaged Navier-Stokes equation to simulate the effect of fluid flow on a cell of particles, thereby smearing the effect of the coupling over a large set of particles. Figure 8.8 shows a two-dimensional schematic implementation of the averaged Navier- Stokes equation (continuum analysis) coupled with the DEM (discrete analysis). This approach generally works well in two dimensions, but its scale may be limited by its high computational cost. Fine-grid approaches rely on the solution of the full Navier-Stokes equation at the pore or grain level using numerical techniques such as the SPH method (Potapov et al., 2001) or the lattice-Boltzman method, where the Boltzman equation is solved over a fixed grid (lattice) that is finer than the particle size (Cook et al., 2004; Sun et al., 2011a). The lattice-Boltzman approach can be shown to capture the Navier-Stokes equation (Sukop and Thorne, 2006). Fine-grid approaches offer higher fidelity, but they require significantly higher computational effort than do the alternatives, making them best suited for highly detailed investigations rather than large- scale simulations. Further connection to the continuum scale emanating from DEM approaches necessitates bridging techniques usually called multiscale approaches (Liu et al., 2006; Zhodi and Wriggers, 2005). DEM approaches rely typically on homogenizing or averaging (i.e., coarse-graining) the interparticle contact forces using average effective stress tensor expressions such as the one proposed by Christoffersen and colleagues (1981): see, for example, (e.g., Nicot et al., 2005; Wellmann, et al., 2008; Andrade and Tu, 2009; Nitka et al., 2011; Guo and Zhao, 2014). On the other hand, fluid flow, if resolved using fine-grid approaches, necessitates coarse-graining to produce appropriate variables such as pore-fluid pressures (see, e.g., White et al., 2006; Sun et al., 2011b). Multiscale methods connect these disparate representations of material behavior across scales, from the grain scale to the field scale (Zhodi and Wriggers, 2005). Multiscale approaches are finding increased use in geomechanics to inform continuum models based on grain-scale resolution rather than relying on purely empirical relationships for material property evolution (e.g., strength, permeability). Multiscale models are much less computationally expensive than are full-resolution discrete mechanics approaches since they resolve the fine scale at selected locations and at specific instances in time when such level of resolution is needed (e.g., in areas of intense deformation). It is expected that with the increasing level of computational power, use of multiscale approaches will increase in research and practice to provide insight into the physics and mechanics of granular materials and, in particular, of liquefaction phenomena. PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW

188 STATE OF THE ART AND PRACTICE IN THE ASSESSMENT OF EARTHQUAKE- INDUCED SOIL LIQUEFACTION AND ITS CONSEQUENCES FIGURE 8.8 Schematic of the average Navier-Stokes (continuum) coupled with DEM (discrete) approaches used to simulate liquefaction behavior in granular soil deposits. SOURCE: Zeghal, M., and U. El Shamy. 2004. A continuum- discrete hydromechanical analysis of granular deposit liquefaction. International Journal for Numerical and Analytical Methods in Geomechanics 28(14):1361â1383. Hybrid Approaches It may be possible to develop hybrid techniques that use different constitutive models and numerical techniques before and after the onset of liquefaction or large displacements. To develop such methods, however, criteria indicating when the alternative constitutive models and numerical techniques should be applied to cope with large deformations and different material response after liquefaction triggering are needed (Nova, 1994; Borja, 2006; Andrade, 2009; Buscarnera and Whittle, 2013). These criteria are typically based on material (internal) stability analyses such as those proposed by Hill (1958). Some hybrid approaches offer much promise and need to be developed further so their fullest potential can be applied to liquefaction problems. In addition, a new generation of computational models needs to be developed that can reproduce physical tests across scales ranging from the laboratory to the field. This can be achieved by basing computational methods on physics and mechanics to establish a more intimate coupling between computation, modeling, and experimentation. PREPUBLICATION VERSION â SUBJECT TO FURTHER EDITORIAL REVIEW