U21A-01 INVITED 08:30h
Modeling of global seismic wave propagation on the Earth Simulator
Accurate modeling of seismic wave propagation in fully three-dimensional (3-D) Earth models is of considerable interest in seismology in order to determine both the 3-D seismic wave velocity structure of the Earth and the rupture process during a large earthquake. We use the Earth Simulator to model seismic wave propagation resulting from large earthquakes on a global scale. Simulations are conducted based upon the Spectral-Element Method (SEM), a high-degree finite-element technique. We include the full complexity of the Earth, i.e., a 3-D seismic wave velocity and density structure, a 3-D crustal model, the elliptical figure of the Earth as well as topography and bathymetry. We use 507 nodes of the Earth Simulator (4056 processors) and model the three dimensional Earth with 13.7 billion grid points. A total of 7 terabytes of memory is needed. We have reached a vectorization ratio of 99.3%, and attained a total performance of 10 teraflops (32% of the peak performance). The very high resolution of the mesh allows us to calculate theoretical seismic waves for realistic fully 3-D Earth model, which are accurate at periods of 3.5 seconds and longer. Usual seismological algorithms, such as normal-mode summation techniques that calculate quasi-analytical synthetic seismograms for one-dimensional (1-D) spherically symmetric Earth models, are typically accurate down to 8 seconds. In other words, the SEM on the Earth Simulator allows us to simulate global seismic wave propagation in fully 3-D Earth models at periods shorter than current seismological practice for simpler 1-D spherically symmetric models. We have modeled seismic waves generated by recent large earthquakes and show the importance of including 3-D structure in the seismic wave simulation. The results demonstrate that given a detailed earthquake source model, precise models of the Earth's mantle and crust, a precise numerical technique, and a large computer, seismic waves generated by large earthquakes propagating through heterogeneous Earth models, which span an amplitude range that covers several orders of magnitude and a few decades in frequency, can be accurately modeled. In particular, we have successfully attempted for the first time an independent validation of an existing 3-D Earth model. These theoretical seismograms calculated for fully 3-D Earth models should improve our understanding of earthquake source mechanisms and the Earth dynamics.
U21A-02 INVITED 08:55h
Forecasting Crustal Activities at Plate Subduction Zones Through Advanced Computing
Plate subduction at ocean trenches causes diverse crustal activities such as earthquakes, volcanism, and mountain building. Among them, earthquakes may be the most interesting and serious phenomena for us. In general, the entire process of earthquake generation consists of tectonic loading due to relative plate motion, quasi-static rupture nucleation, dynamic rupture propagation and stop, and fault strength restoration. Recent progress in the physics of earthquake generation enables us to quantitatively describe the entire process of earthquake generation by a coupled nonlinear system, which consists of a slip-response function and a fault constitutive law. The driving force of this system is observed relative plate motions. The system to describe the earthquake generation cycle is conceptually quite simple. The complexity in practical modeling mainly comes from complex structure of the real Earth. The aim of our research program (CAMP) is to develop a physics-based, predictive simulation system for the entire process of earthquake generation in and around Japan, where the four plates of Pacific, North American, Philippine Sea, and Eurasian are interacting with each other in a very complicated way. The total simulation system consists of a crust-mantle structure model, a tectonic loading model, and a dynamic rupture model. First, we constructed a realistic 3D model of plate interface geometry in and around Japan by applying an inversion technique to ISC hypocenter data. A 40 km thick elastic surface layer overlying a Maxwellian viscoelastic half-space models the crust-mantle rheological structure. With this structure model, given relative velocities among the plates, we can compute long-term crustal deformation due to steady plate subduction. The results demonstrate that the steady subduction of oceanic plates brings about steep uplift at island arcs, sharp subsidence at ocean trenches, and gentle uplift at outer rises. Second, by integrating microscopic effects of abrasion and adhesion at contacting rock surfaces, we theoretically derived a rational fault constitutive law, called the slip- and time-dependent law. Incorporating this constitutive law into the steady plate subduction model, we can develop a quasi-static model of stress accumulation driven by relative plate motion. By applying the boundary integral equation method, we have also developed a simulation model for dynamic rupture propagation on a 3D curved plate interface. Thus, combining the quasi-static stress accumulation model and the dynamic rupture propagation model on the same structure model, we can finally construct a unified simulation system for the entire process of earthquake generation. Outputs of the simulation system are the crustal deformation and internal stress change associated with seismic and/or aseismic slip at plate interfaces. From comparison of these outputs with observations, we can extract useful information to estimate past slip history and the present stress state by using an inversion technique. Given the past slip history and the present stress state, we can predict the next step fault-slip motion and stress change through computer simulation. As an example, we show the quasi-static stress accumulation process at the 1968 Tokachi-oki seismogenic region, northeast Japan, and the subsequent process of rupture initiation, propagation and stop there. In this simulation we forced dynamic rupture to start by giving artificial stress drop. Then, unstable rupture started, but it was not accelerated. This means that the dynamic rupture is not accelerated, if the stress state has not reached to a certain critical level.
U21A-03 INVITED 09:20h
Numerical Simulations of Convective Dynamos and Electromagnetic Induction in a Three-Dimensional Heterogeneous Earth on the Earth Simulator
Numerical simulations of convective dynamos and electromagnetic induction in a 3-D heterogeneous earth have been performed by using the Earth Simulator. For the geodynamo simulation, three types of simulation code have been developed on the Earth simulator. In the STM (spherical harmonic Spectral Transformation Method) dynamo simulation, dynamo action was confirmed at the Ekman number, E = 10-5, 10-6, and the Rayleigh number (Ra) up to 10Rc, where Rc is the critical Rayleigh number for the onset of convection. In order to perform low Ekman and high Rayleigh number geodynamo simulation, we ensure high spatial resolution in terms of spherical harmonics up to degree 255 and in the radial direction up to 256 grid points. These results were obtained using 64 nodes of the ES at a performance of 1.18 TFLOPS. Results of the simulation indicate the transition of the dynamo state from no magnetic field to unstable polarity reversing dynamos through stable dipole-dominant states as Ra increases. With the FEM (Finite Element Method) simulation code, dynamo simulations were performed with E as low as 5 x 10-5. Besides these two approaches, a new method (Fourier spectral Transformation Method) was developed in order to further decrease the Ekman number. The simulation code for the new method was transferred to the ES and tested with the benchmark test codes for the geodynamo simulation. For the Electromagnetic induction modeling, several simulation codes were installed to the ES. These induction models calculate the electromagnetic induction in a spherical earth with a 3-D electrical conductivity structure in time and frequency domains. Frequency domain responses are commonly used to analyze a steady state part of the geomagnetic variations to detect the conductivity structure of deep interiors of the earth, whereas time-domain approaches are also useful for analyzing the transient responses of the induction due to the external magnetic disturbances observed during large magnetic storms. Among the installed codes, the time-domain approach based on the formulation given by Hamano(2002) achieved a high performance on the ES. This method utilizes the standard decomposition of the magnetic field into toroidal and poloidal parts, and spherical harmonic expansions of both the magnetic fields and the conductivity heterogeneity. A finite difference approximation was used to solve the set of diffusion equations for the spherical harmonics. This method can be efficiently used to analyze transient geomagnetic variations to estimate the 3-D conductivity structure of the earth. Since the method directly solves the temporal variation of the expansion coefficients corresponding to the source field, which is also expanded by the spherical harmonics, the method can work together with the results of the spherical harmonic analysis on the observed record of the geomagnetic field. In the highest resolution model, spherical harmonics up to degree 128 are used and the computing performance is about 60 percent of the peak performance on 160 nodes of the ES. In the frequency domain analysis, a non-linear iterative inversion code on the ES was developed to estimate the 3-D electrical conductivity structure. We applied these codes into the experimental EM data measured in the north Pacific region and estimated the electrical conductivity structure in the mantle transition zone. The conductivity structure indicates that the several features of the large-scale anomalies exist such as a conductive region beneath Hawaii and a resistive region beneath Philippine, and some features are consistent with the seismic wave velocity structure.
U21A-04 09:45h
A Large-Scale Simulation of Earthquake Cycles Along the Nankai Trough, Southwest Japan: Lateral Variation in Dip Angle of the Slab to Control the Nucleation Position
The Pacific and the Philippine Sea plates subduct beneath northeast and southwest Japan, which produces considerable lateral variations in structure and disastrous great earthquakes. Under a collaboration project on the Earth Simulator entitled "Simulation of Earthquake Generation Process in a Complex System of Faults", we are constructing regional 3-D heterogeneous viscoelastic models in northeast and southwest Japan to simulate cycles of great interplate earthquakes along the Japan trench and the Nankai trough, respectively. At the meeting, first we will introduce the present status and the plan of this project. Then we will report an interesting result of a simple but large-scale quasi-static simulation of earthquake cycle along the Nankai trough in southwest Japan as described below. The occurrence pattern of historical great interplate earthquakes along the Nankai trough indicates the existence of rupture segments correlated with 100 km scale structural heterogeneity in the configuration of the subducting Philippine Sea plate. At least, the latest earthquakes, the 1944 Tonankai and the 1946 Nankai earthquakes, initiated their ruptures off the Kii peninsula, where the dip of the slab is steeper than that in the adjacent areas. In order to investigate how such heterogeneity in the slab geometry affects the occurrence of earthquakes, a three-dimensional numerical simulation of earthquake cycles based on a rate- and state-dependent friction law is conducted for the large source area of about 700 km along the trough and 300 km in the dip direction. We simply use a flat dipping plane as a model plate boundary and divide it into 147,456 equal-area cells, each of 1.20 km x 1.22 km. The slip history on each cell is numerically calculated by assuming equilibrium of the frictional stress with the elastic stress. Because of the large number of cells, a large computation power is required, which forces us to use the Earth Simulator. The depths of the plate boundary are determined from seismic surveys and seismicity data. Frictional constitutive parameters and effective normal stress on the plate boundary are assumed to be dependent on the depth of the plate boundary thus determined, and the depth-dependent frictional properties are projected on the model plane. In the simulation, we take into account also the different convergence rate along the Nankai trough. According to a recent GPS data analysis, the rate is 6.5 cm/yr in the western portion and is decreasing eastwards from the Kii peninsula to 2 - 3 cm/yr in the Tokai region. The simulation results show that preslip occurs and rupture starts off the Kii peninsula where the dip angle is higher than that in the surrounding area. In the present model, the width of locked zone during an interseismic period becomes narrower in the area with a higher dip of the slab, and the stressing rate is higher there. Although the slab is dipping most steeply in the Tokai area, the plate convergence rate is so low that the stressing rate is not high. Thus, the stressing rate is the highest off the Kii peninsula and simulated ruptures start there. Our simple simulation suggests that the lateral variation in dip angle of the slab is one of controlling factors for the nucleation of great interplate earthquakes.