SciDAC Review - Computational Predictions of EARTHQUAKE Ground Motions and the Response of MAJOR STRUCTURES

GEOPHYSICAL SIMULATIONS

Computational Predictions of EARTHQUAKE Ground Motions and the Response of MAJOR STRUCTURES

In the approximately 100 years since the Great San Francisco Earthquake of 1906, cutting-edge computational capabilities have progressed to the point where complex physical processes associated with earthquakes can be numerically simulated. Realistically simulating ground motions, and the response of complex structures and systems, using high-performance computing can help investigators improve seismic safety and the integrity of critical structures, such as bridges and nuclear power plants.

Earthquakes cause tremendous devastation and suffering everywhere they occur. The Great San Francisco Earthquake of 1906 was a defining event that provided new scientific insight into the mechanics of earthquake processes. The 1906 earthquake, and the fire that followed, occurred in an area that was only beginning to be broadly urbanized and contained few large and complex structures and systems. Today's situation is very different: major cities in seismically active locations, including the San Francisco Bay Area of California, are densely populated and jammed with skyscrapers, elevated highways, and complex subterranean infrastructure systems. As human populations continue to increase dramatically worldwide, and urban centers expand in seismically active regions, the potential for earthquakes to stress complex structures and systems increases just as dramatically. In seismically active regions, infrastructure and building designers are challenged to accurately estimate the likely earthquake motions over the lifetime of a given structure. Once those ground motions are properly identified, the designers then must apply appropriate seismic design features so that the structure can withstand the expected ground motions without catastrophic failure or unacceptable damage.

Today's high-capacity, high-capability computing platforms (HPC machines) are sophisticated tools for investigating the effects of earthquakes anywhere on the planet. Engineers can test components of structures (for example, a single column or a modest-size substructure assembly) under simulated earthquake motions using mechanical shaking tables to gain important insights into substructure and component performance. Prior to 1980, however, computational limitations constrained early simulations to very low-order simulation models that invoked many idealizations in order to render the model computationally tractable, and to approximately simulate the seismic response of large, complex structural systems. In that technological environment, the final assessment of seismic performance was often discovered in actual, infrequent earthquake events, when the structures either performed adequately or failed dramatically (figure 1).

As human populations continue to increase dramatically worldwide, and urban centers expand in seismically active regions, the potential for earthquakes to stress complex structures and systems increases just as dramatically.

Figure 1. The final assessment of seismic performance is often discovered in an actual earthquake, when structures either perform adequately or fail dramatically. These major transportation structures, (upper) an elevated expressway, following the 1995 Kobe Earthquake in Japan and (lower) a viaduct highway in Oakland, CA (the "Cypress Structure"), following the 1989 Loma Prieta Earthquake, failed dramatically.

Over the past three decades, an improving ability to construct accurate and detailed nonlinear finite-element models of complex structures and systems has increased tremendously (see a 3D model of the most heavily traveled bridge in the United States in figure 2, p44). These 3D models can predict the dynamic response of structures to imposed earthquake ground motions. The models now include appropriate interactions among all parts of the structure, and nonlinear damage models, so that structural nonlinearities and damage states can be realistically represented. The latest computational models, enabled by HPC machines, provide enhanced insight into the complex linear and nonlinear dynamics of structural systems. They offer an extremely powerful tool for subjecting structures to earthquakes in virtual space, where it is far better to discover structural shortcomings and even to collapse a structure in computer memory than in Downtown Anywhere. Advanced simulations can reduce the uncertainties in predicting structural performance during real earthquakes. Thus, it is now accurate to state that computational modeling has truly transformed the fields of structural mechanics and earthquake engineering, and has become fundamental to modern seismic structural design.

The latest computational models, enabled by HPC machines, provide enhanced insight into the complex linear and nonlinear dynamics of structural systems.

The prediction of earthquake ground motions, in contrast, persists largely as an empirically dominated process, with estimates of future earthquake motions at a specific site derived from the ground motion measurements of past earthquakes at many vastly distributed sites. This approach conceals significant limitations. From extensive observation, scientists now understand that the complexities of geophysical processes result in localized, very site-specific dependence of earthquake ground motions. To enhance understanding and assist in developing realistic site-specific motions, researchers are compelled to seek geophysics-based simulation models that can accurately represent the complex geophysical conditions and processes at the site of interest. Because of the extraordinary computational demands associated with regional 3D modeling of the Earth, such simulations had previously exceeded the capabilities of even high-end computing resources. Top-end, massively parallel HPC machines now enable a tremendously exciting computational regime that permits, for the first time, simulation-based predictions of future ground motions.

Figure 2. The five elements of the bridge system model are: (a) finite rotation fiber flexure element for the bridge towers; (b) penalty node-to-node contact for bridge impact; (c) tension-only, two-force member with initial stress for cables; (d) composite membrane, truss, and sway stiffness deck model; and (e) caisson block with uplift.

This article examines the challenges of and opportunities for creating accurate HPC seismic simulations by examining two case studies. The first looks at a virtual geophysics reconstruction of the Great San Francisco Earthquake of 1906, an effort focused on evaluating the ability of cutting-edge, high-performance geophysics models to inform us about the scientific aspects and hazard implications of that historically important event. The second example demonstrates an end-to-end simulation that highlights the ability to simulate from fault rupture through to detailed structural response. Such end-to-end simulations, which employ the merging of geophysics and engineering, provide an important new tool to understand, and appropriately design, complex and expensive structural systems.

The Great San Francisco Earthquake of 1906 Reconstructed
At 5:12 on the morning of April 18, 1906, a magnitude 7.8 (M7.8) earthquake struck northern California and propagated seismic waves across the entire western United States. This event occurred along the San Andreas Fault, with an epicenter (point of earthquake origin) located a few kilometers off the coast of San Francisco. When it was over, a 480 km segment of the San Andreas Fault had failed. It has been estimated that more than 3,000 people were killed in what was then a relatively sparsely populated region of California. Known as the Great San Francisco Earthquake of 1906, this event marked the beginning of modern earthquake science. For the first time, a large earthquake and its effects were carefully observed and systematically documented, as a significant basic understanding of earthquake phenomenology was being developed.

To commemorate the centennial of the 1906 earthquake, and to investigate the effects of that earthquake with the best emerging computational tools, a collaborative modeling and simulation effort was engaged to recreate and study the strong ground motions generated by this event. The United States Geological Survey (USGS) organized this multi-institutional effort to include simulation teams and other participants from several institutions and organizations, including the U.S. Department of Energy National Labs and the University of California. In addition to the 1906 earthquake modeling effort, simulations were performed for scenario earthquakes along the San Andreas and other faults in the San Francisco Bay Area of California. Simulations of these hypothetical events predicted ground motions from future earthquakes. The complete simulation suite consisted of three simulations of the 1906 earthquake and seven simulations of other scenario events along the San Andreas Fault. The simulations differed by the seismic source-term used to describe the rupture mechanism along the fault, that is, the scenarios each considered different points of fault rupture initiation and fault rupture propagation.

Top-end, massively parallel HPC machines now enable a tremendously exciting computational regime that permits, for the first time, simulation-based predictions of future ground motions.

In an independent effort, simulation code validation studies were performed by comparing observed and simulated ground motions for the 1989 Loma Prieta Earthquake, which occurred in a mountainous region near Santa Cruz, approximately 60 km south of San Francisco. See the sidebar on verification and validation (V&V; sidebar "Verification and Validation Efforts,"p47) for a discussion of how the regional-scale geophysics codes are verified and validated.

A 3D geologic model of northern California was used for the 1906 simulations. This model, which was recently developed by the USGS, utilizes all available sources of sub-surface data to describe the constitutive properties of the Earth's crust and upper mantle throughout the domain of interest. The computational model domain is 630 km long and 320 km wide. It is oriented in a northwest-southeast direction, roughly parallel to the California coast and the major tectonic features in the region. The model extends from the surface to about 50 km depth, and includes surface topography. The constitutive properties in the model are represented by the compressional and shear velocities of seismic wave propagation, the material densities, and damping factors that define the material-dependent, non-geometric dissipation of seismic energy. The model includes many of the geologic features important for earthquake ground motion prediction, such as low-velocity near-surface sediments and regional basins.

Figure 3. Comparison of finite-element (CMUN) and finite-difference (UCBL) velocity time histories (N-S component) for one of the Sierra Madre fault scenarios. Results are shown at 16 sites.

Each of the teams discretized the geologic model onto either a finite-difference or finite-element grid. The computational domain used by each group spanned all or part of the geologic model. The size of the domain was largely constrained by the computing power available to each team. Some groups included the effects of surface topography and seismic attenuation in their simulations, while others did not. Spatial discretizations between 50 m and 200 m were used, which are appropriate for resolving seismic simulation frequencies of up to 1 Hz. Each research group performed between two and ten simulations.

As an example, one modeling team used the E3D seismic wave propagation code (as well as a public domain counterpart) to simulate the 1906 earthquake and associated events. Selected results from these simulations are shown in figures 5 (p48), 6 (p49), and 7 (p50). In particular, E3D uses a staggered-grid, finite-difference scheme to approximate the elastodynamic formulation of the full wave equation. The numerical accuracy is fourth-order in space and second-order in time. The effects of seismic attenuation and surface topography are included in the model. E3D is used at over 100 government, academic, and commercial institutions, and has been validated in about three dozen studies. It is used in a number of application areas, including earthquake research, solid earth geophysics, petroleum exploration, nonproliferation, carbon sequestration, hydroacoustics, and intelligence-related efforts.

E3D is used in a number of application areas, including earthquake research, solid earth geophysics, petroleum exploration, nonproliferation, carbon sequestration, hydroacoustics, and intelligence-related efforts.

Figure 4. Differences between the simulated ground motions computed with E3D and the ground-shaking intensities observed for the 1906 San Francisco Earthquake. The simulated ground motions match the observed data over much of the model (shown by white and yellow), although the agreement is less pronounced in some locations (shown by orange, red, and blue). The fit is particularly degraded in the areas north of Sacramento. Differences between the simulated and observed ground motions are attributed to inaccuracies in the geologic model, sparse data coverage over parts of the domain (blue dots), and the empirical nature of the relationship that translates qualitative measures of felt shaking into quantitative values. Overall, however, the simulated ground motions fit the observed data very well.

For the earthquake simulations presented here, the spatial discretization of the finite-difference grid was 100 m, with grid dimensions spanning the entire 630 x 320 x 55 km domain of the USGS 3D geologic model. The discretized model contained approximately 11 billion grid nodes and required about 1 TB of computer memory. Parallelization was achieved through domain decomposition. MPI was used to communicate information between the boundary nodes of each subgrid. Three HPC systems at Lawrence Livermore National Laboratory (LLNL) and the Earth Simulator in Japan were used to complete the simulations. Several simulations were performed using 128 nodes (1,024 CPUs) of an AMD Opteron Linux cluster (with 1,152 total nodes). A sustained performance of about 22% peak was observed.

In addition, high-resolution test simulations using grid discretizations of 50 m and 75 m were used to evaluate the numerical accuracy of the results. These simulations utilized up to 40 billion grid nodes and about 4 TB of memory. It was found that the finer discretizations did not significantly alter the computed ground motions and the numerical results were judged to converge for frequencies up to 1 Hz. The resolution of 1 Hz frequency motions is supported by the geologic data available. In the quest to resolve higher frequencies, which may be important for some complex structures and systems, special approaches will be needed to define the geology at shorter spatial resolutions as discussed below.

Figure 5 (p48) presents snapshots (at 30 second intervals) of seismic energy propagating throughout northern California for the preferred simulation of the M7.8 1906 earthquake. This scenario, in which the San Andreas Fault rupture begins just off San Francisco with the rupture propagating bilaterally both north and south, has been accepted by the scientific community as the most likely manner in which the 1906 earthquake actually progressed. In addition, this rupture simulation best fits the limited seismic and geophysical data obtained from the 1906 event. Shown is the seismic wave front as it expands out from both the epicenter and propagating rupture along the San Andreas Fault. The color scheme indicates the Modified Mercalli Intensity (MMI), which is an intensity measure of the maximum ground shaking during an earthquake. Orange and red indicate areas of strong to extreme motion, where heavy damage to structures can be expected. The strongest computed shaking occurs in a band surrounding the San Andreas Fault, which corresponds with the historical observation that severe damage was restricted to a very narrow zone (in terms of east-west damage) for such a massive earthquake. This is consistent with the high attenuation rates observed in the fractured, dissipative geology of the western United States. In addition, heavy shaking is observed in the simulations for several of the sedimentary (soft soil) basins in the region. The low-velocity sediments in these basins cause seismic ground motions to be amplified. Additionally, basins often trap seismic energy thereby creating long-lasting reverberations in the ground, which can place increased demands on a complex structure or system by exposing it to many more oscillatory cycles. The most significant damage often occurs in these regions, as was the case in 1906 when some of the most intense shaking occurred inland in the Santa Rosa Basin (just north of San Francisco Bay). In the last snapshot in figure 5, for example, a red zone indicative of extreme ground shaking extends inland into the Santa Rosa area, which is consistent with the extreme damage observed in that city in 1906.

High-resolution test simulations using grid discretizations of 50 m and 75 m were used to evaluate the numerical accuracy of the results. These simulations utilized up to 40 billion grid nodes and about 4 TB of memory.

Figure 5. Three-dimensional geophysics simulation of the 1906 San Francisco Earthquake (actual bi-lateral rupture sequence, M=7.8).

Scientifically, inspection of these simulation results has provided better understanding of the actual magnitude and the fault rupture process that occurred during the 1906 earthquake and has helped delineate the geophysical phenomena causing anomalous ground shaking patterns, as well as indicating localized areas of expected intense motions in future earthquakes. Because the 1906 earthquake occurred before the invention of strong ground motion instruments, the modern simulations provide the first quantitative insight into the amplitudes of ground motions throughout the affected region.

Figure 6 presents results from a second simulation of a M7.8 scenario earthquake along the San Andreas Fault. The magnitude and fault rupture process are identical to the 1906 simulation, except that the epicenter is located near the northern terminus of the fault. Ironically, even though this scenario event begins 200 km to the north of the 1906 epicenter (which was 3 km southwest of San Francisco), the predicted ground motions in the San Francisco Bay Area are significantly more severe than those from the actual 1906 earthquake. The phenomenon exhibited in this simulation scenario is known as directivity of motion. In essence, the seismic amplitudes are significantly increased in the direction of fault rupture and wave propagation. Although directivity is a well-known process, these scenario simulations help demonstrate how very significant this effect can be. In fact, in some San Francisco Bay Area locations, the computed ground motion amplitudes from this scenario simulation are more than five times larger than those for the baseline 1906 earthquake rupture scenario of figure 5, due to the directivity of seismic waves. This indicates that, as destructive as the 1906 earthquake was, it could have been much worse for the San Francisco Bay Area region.

Figure 6. Three-dimensional geophysics simulation of the 1906 San Francisco Earthquake (hypothetical unilateral rupture sequence, M=7.8).

The simulation in figure 6 also illustrates clearly how seismic energy can be redirected and focused toward a particular region. In particular, energy is channeled from the San Francisco Bay Area into the eastern Bay Area and toward California's Central Valley. This energy refocusing is caused by the 3D geologic structure in the region, and indicates that for a San Andreas earthquake rupturing north-to-south, the effect of extreme motions would be felt much more in the eastern Bay Area and further inland than for the 1906 event. The probability of such a large scenario event occurring within the next 30 years is about 1%—small, but not insignificant.

These simulations shed new light on and suggest new interpretations for the 1906 earthquake, and assist in developing a global picture of the regional seismic hazard in northern California. Most importantly, once validated, these simulation methodologies can be adapted to any region of high seismicity and high seismic hazard to assist in understanding the risks to complex structures and systems. The scientific discoveries begun with the extensive field observations just after the 1906 earthquake is echoed today in the process of numerical discovery enabled by HPC platforms that allow dynamic, 3D inspection of these major events for the first time.

End-to-End Simulations: Determining the Risk to Complex Structures and Systems
To fully exploit high-performance geophysics-based simulations of earthquake motions, it is necessary to take the simulation process one step beyond the assessment of ground motions and determine the effects of such motions on major structures. The Holy Grail of seismic simulations is the end-to-end simulation that starts with the initiation of fault rupture, progresses through wave propagation and ultimately to the dynamic structural response indicating structural deformations and forces at the individual component level. Such simulations merge the geophysics models for wave propagation directly with the detailed mechanics models for structural response simulations, for example with the finite models of a bridge structure, as shown in figure 8 (p52).

Such end-to-end simulation processes have recently been completed to assist in understanding the seismic response of long-span bridges that are sited near major earthquake faults. Long-span bridges are flexible structures that can exhibit very long periods of natural vibration (for example, 10-15 s fundamental periods) of the deck system. In the last 10 years it has become clear, based on a sparse set of ground motion measurements with new broad-band accelerometers, that for sites near major earthquake faults (10-15 km from the fault), the ground motions can exhibit long-period components of motion. Such motions are associated with the tectonics of fault movement and the superposition of propagating seismic waves to form a "pulse" movement of the ground. Because of the potential for strong interactions between long-period structures and long-period ground motions, these motions, which have not been included in historical seismic hazard assessments, are potentially damaging to long-period structures. Because of the scarcity of adequate broadband, near-field records from large earthquakes, there is significant motivation and benefit to exploring this phenomenon with HPC simulation models run on the fastest machines.

Because of the scarcity of adequate broadband, near-field records from large earthquakes, there is significant motivation and benefit to exploring this phenomenon with HPC simulation models run on the fastest machines.

In the San Francisco Bay Area, the San Francisco-Oakland Bay Bridge suspension span connects the San Francisco Peninsula with Yerba Buena Island (West Bay Crossing) and continues on its East Bay Crossing to Hayward, CA. It is located approximately 15 km west of the Hayward Fault. This structure, completed in 1936 (well before modern understanding of earthquake effects), is a critical transportation link that carries more vehicles per day than any other U.S. bridge (approximately 280,000 vehicles per day). This bridge was selected for a case study into effects of long-period motions, and a special-purpose nonlinear model was developed to represent the response of the bridge structure (figure 8, p52). This model included specially developed finite elements for the deck and foundation system that were tailored to modeling the overall suspension bridge response with the fewest possible degrees of freedom. The entire nonlinear bridge finite element model in figure 8 (p52) included less than 25,000 degrees of freedom.

Figure 7. Three-dimensional geophysics simulation of the Hayward Fault Earthquake (hypothetical bilateral rupture sequence, M=7.1).

To investigate the effects of near-field motions on this long-span bridge, ground motion simulations were executed using the regional Bay Area geophysics model to simulate the motions from earthquakes on the nearby Hayward Fault. The Hayward Fault extends about 80 km along the eastern side of the San Francisco Bay and the maximum earthquake expected from the Hayward Fault is approximately M7.1. It should be noted that this type of earthquake is very important from a structural design viewpoint, because the probability of occurrence is much higher than a great earthquake like the 1906 event. In fact, the predicted probability of a significant (greater than M6.7) earthquake along the Hayward Fault within the next 30 years is about 25%, making it more likely that existing and new complex structures and systems in the region will be subjected to such an event during their effective lifetimes.

Figure 7 and figure 8 (p52) show the simulated wave pattern and the maximum ground amplitudes from a M7.1 scenario earthquake along the Hayward Fault that ruptures in bilateral fashion (that is, the fault rupture initiates in the middle and propagates in both directions). While the spatial distribution of shaking is not as extensive as that from a larger earthquake along the San Andreas Fault (comparing figure 5 and figure 7 provides a clear visual picture of the significant difference between M7.1- and M7.8-earthquake representations in the logarithmic magnitude scale), the computed strong ground motions near the fault rupture are still quite large and potentially very damaging.

Figure 8. This end-to-end simulation, from fault rupture to structural response, yields synthetic ground motions at the San Francisco Bay Bridge resulting from an M=7.1 earthquake along the Hayward Fault. This demonstrates nonlinear analysis of bridge systems.

Various scenarios of fault rupture were simulated, with the synthetic ground motions at the Bay Bridge site being generated directly from the geophysics model and applied to the base of the structural model of the bridge. Figure 8 (p52) shows graphical plots of the bridge site ground motions associated with one particular rupture scenario of the Hayward Fault. Notable is that the ground displacement and ground velocity exhibit a pronounced pulse indicative of near-field effects. In addition, the ground displacement also exhibits a permanent ground offset resulting from the northward motion of the Earth west of the Hayward Fault. Both types of motions were absent in traditional earthquake strong ground motion measurements because of limitations in the ground motion instruments and the data processing techniques applied. Historically, any motion of 5 s or greater was excluded from the final processed time histories of ground motion. Only with the advent of modern digital ground motion instruments have these long-period waveforms been measurable, so the existing database of such records is still quite small.

To evaluate the significance of the long-period, near-field components of ground motion, the bridge structure response was analyzed in two ways; first with the full synthetic records from the geophysics simulations including the long-period motions, and second with the synthetic records band-passed to remove the long-period components. The results of these two simulations were then compared to assess the impact of the long-period motions on the bridge system. Figure 8 shows a snapshot in time of the bridge displacement under the earthquake ground motion excitations (displacements amplified 10 times).

Evaluation of the bridge element forces and inspection of the graphic, animated response of the bridge system under these two sets of ground motions clearly indicate that the long-period motions associated with the near-field records can have a very substantial contribution to bridge forces. For the Hayward Fault scenarios, the forces in selected bridge components were found to be approximately two times greater when the near-field, long-period components were included. Alternatively, if the near-field ground motion components are not included in the hazard assessment, the bridge forces may be underestimated by a factor of two. These observations have important implications for adequate seismic design and for ensuring that issues are identified in virtual space rather than tragically in the field during the next major earthquake.

This example illustrates the understanding and knowledge gained through simulation of complex processes, with HPC being the great enabler. We have discovered that much can be learned, and significant understanding and intuition increased, by simply inspecting computer graphics animations of these very complex processes.

Challenges and the Path Forward
The power of HPC and the efficacy of advanced simulations in modeling earthquake ground motion have been clearly demonstrated and are supported by robust V&V efforts. The data gained from these simulations are now ingrained in today's earthquake engineering and structural design. One of the next major advances in earthquake hazard assessment will be to develop and implement the appropriate framework for exploiting emerging high-performance 3D geophysics simulations.

High-performance simulations can be leveraged to gain insight and understanding into complex processes of wave propagation and accompanying structural response. This would include the important insights gained by developing an understanding of cause-and-effect between subsurface processes and structural response. For example: Which geophysical and fault rupture parameters have the most bearing on structural response? The answers to such questions will help investigators focus and prioritize future research in seismic hazard assessment.

Future frameworks will address the uncertainties and range of outcomes consistent with known geophysical data. We do not know a priori precisely how a fault is going to rupture for a particular earthquake and the motions for two equal magnitude earthquakes with different rupture scenarios may be substantially different. Decisions with regard to how many rupture scenarios will be "enough" to capture the range of possible outcomes must be made. Investigators must also determine how many high-performance realizations should be executed to constitute a deep enough data bank.

For certain complex structures (nuclear power plants in particular), components of the system will be sensitive to high-frequency motions, in the range of 5-10 Hz. However, available geologic data will not be at fine enough spatial resolution to permit deterministic characterizations of the geologic properties on these fine scales. To fully exploit geophysics models in engineering assessments, researchers must agree upon a strategy for addressing the high-frequency issues. A suite of stochastic geology representations may well be required to permit assessment of the potential range of higher-frequency simulations. Another approach would merge deterministic high-performance simulations over one frequency range (say up to 1 Hz) with motions generated by other means; random vibration at high frequencies or traditional spectrum-compatible motions from a traditional seismic hazard analysis.

The ability to compute seismic wave propagation in 3D on a regional basis is a tremendously powerful new tool. Prospective research will determine the best approaches to implementing this tool in a framework that allows appropriate definition of the seismic hazards, and enables the design of complex structural systems.

In the future, the appropriate combination of geophysics and structural mechanics models that provide an end-to-end simulation framework to assist hazard and risk assessment will provide a powerful new tool for seismic evaluations. Such an approach has the potential benefits of accounting for the site-specific geophysical setting at a particular structural site, and for evaluating cause-and-effect between earthquake fault rupture parameters and structural response. An advanced simulation approach should be developed and implemented initially within the context and framework of traditional hazard assessments.

Thus, much work remains to be done in communicating the details of advanced simulations to the traditional hazard community. HPC machines will play a central role in enabling this future work in earthquake ground motion modeling.

Contributors:
David McCallen and Shawn Larsen, LLNL