Tsunami Analysis Method with High-Fidelity Crustal Structure and Geometry Model

Higher fidelity seafloor topography and crustal structure models have become available with accumulation of observation data. Previous studies have shown that the consideration of such high-fidelity models produces significant effects, in some cases, on crustal deformation results that are used as inputs for tsunami analysis. However, it is difficult to apply high-fidelity model of crustal deformation computations to tsunami computations because of large computational costs. In this paper, we propose a new crustal deformation computation method for estimating inputs for tsunami computations, which is based on a finite element analysis method with remarkable reduction of computation costs by efficient use of the arithmetic space and the solution space. This finite element analysis method enables us to conduct 102−3-times crustal deformation computations using high-fidelity models with a degree of freedom on the order of 108 for the 2011 Tohoku earthquake example. Tsunami computations with typical settings are conducted as an application example to present the advantages and characteristics of the proposed method. Comparisons between results of the proposed and the conventional method reveal that large shallow fault slip around the trench axis may lead to significant differences in tsunami waveforms and inundation height distributions in some cases.


Introduction
Numerical analysis methods for tsunamis have been developed for the purpose of reducing tsunami disasters and further understanding of the geophysics of tsunami phenomena. Such methods have also been applied to reproduce observed tsunamis caused by megathrust earthquakes [e.g., Suppasri et al., 2012;Satake et al., 2013;Ren et al., 2015]. In these studies, researchers typically have used an analytical solution in a homogeneous half-space [e.g., Okada, 1985;Mansinha and Smylie, 1971] to compute crustal deformation, and then have performed hydrodynamic computation of the water level di®erence in propagation from the computed crustal deformation [e.g., Satake, 1995;Titov and Gonzalez, 1997;Hebenstreit, 2001]. On the other hand, previous research [e.g., Masterlark, 2003] has shown that consideration of a heterogeneous crustal structure and geometry of surfaces and layer interfaces has a signi¯cant e®ect on computation of the crustal deformation in some cases. In addition, several studies [e.g., Hyodo and Hori, 2011] have demonstrated the e®ect of detailed geometry around a trench axis, i.e. whether the fault surface reaches the trench axis, on the analysis results. A subduction earthquake may have large shallow fault slip [e.g. Fujiwara et al., 2011]; thus, such an e®ect may have a large impact in analysis for such earthquakes. Therefore, tsunami computations based on computed crustal deformation results, considering high-¯delity crustal structure data that faithfully reproduce observed data, including detailed geometry around the trench axis, is desirable.
Because of the progress made in observational techniques and analysis methods, such high-¯delity crustal structure data have been accumulated. For instance, for the area in and around the Japanese Islands, the plate geometry is estimated [Hashimoto et al., 2004]; the velocity structure of northeastern and western Japan, at 1 km resolution in the most detailed part [Koketsu et al., 2008]; and the estimated trench axis geometry for northeastern Japan [Nakamura et al., 2013] are available. Such high-¯delity crustal structure data and geometry are also being developed in other subduction zones. Considering this advancement, several studies have carried out tsunami analysis based on crustal deformation computation using¯nite element (FE) computations, which can properly consider heterogeneous and complex crustal structures [e.g. Grilli et al., 2013]. However, the resolution of the FE models used in these studies was relatively low because the construction of an FE model with detailed modeling of geometry around the trench axis and computations using the model requires large computational costs. Here, we roughly estimate the computational cost of using the highest-resolution (i.e. 1 km) crustal structure data available for the area in and around the Japanese Islands. The degree of freedom (DOF) of the FE model is at least on the order of 10 8 , if we seek to guarantee the convergence of the solution in a domain of 10 3 Â 10 3 Â 10 2 km, which includes the focal region of an M9-class megathrust earthquake (e.g. the 2011 Tohoku-Oki earthquake). In a tsunami analysis, the input data for the hydrodynamic computations are the temporal history of crustal deformation due to a certain fault rupture process. This history is computed by superimposing Green's functions (GFs) due to unit small faults in the fault surface based on the fault rupture process. The fault surface of a megathrust earthquake is typically divided into 10 2À3 small faults. Therefore, 10 2À3 crustal deformation computations using an FE model with a DOF on the order of 10 8 are required. This approach was previously considered to be di±cult to realize because of the large computational cost.
In this study, we propose a tsunami analysis method enhanced by a high-¯delity crustal structure model and geometry. We compute 10 2À3 GFs with an FE model with a DOF on the order of 10 8 based on high-¯delity crustal structure data. Then, we compute the temporal history of crustal deformation by superimposing the GFs based on the possible fault processes in the target location. Finally, we carry out hydrodynamic computations for tsunami by inputting the computed crustal deformation. The remainder of this paper is organized as follows. Section 2 describes the methodology in detail. In Sec. 3, we provide examples of tsunami analyses for realworld problem settings and show the advantage of the proposed method by comparing the results with those obtained with a conventional method.

Method
We¯rst compute the temporal history of sea°oor elevation due to crustal deformation, which is used as an input for hydrodynamic computations of tsunamis. The crustal deformation due to fault dislocation in location x at time t is written as: where u i , S, g ij , and u f j are the ith components of response displacement, fault surface, GF, and slip on the fault surface in the jth direction, respectively. Here, we describe the formulation more generally to consider temporal change for earthquakes with large focal regions, in which the fault slip does not occur uniformly, but gradually propagates through the fault surface. Because crustal deformation can be considered as a quasi-static process in tsunami analysis, g ij is a static GF.
As mentioned in Sec. 1, conventional tsunami analysis involves the computation of g ij in Eq. (1) by simplifying the crustal structure and geometry by means of an elastic homogeneous half-space. The FE method can reproduce the complex geometry in the analysis and automatically satis¯es traction-free boundary conditions. Therefore, the FE method is a proper tool for crustal deformation computation that considers the heterogeneity of crustal structure. We seek to realize Eq. (1) by computing the response displacements G ij due to unit small faults in the FE analysis and superimposing the small faults. Here, G ij is due to unit fault slips in two directions, the strike-slip direction and the slip direction perpendicular to it, in each small fault, and the entire fault surface is divided into n f small faults. In this case, n f Â 2 sets of linear equations, as shown below, must be solved: where K, G ij , and F ij are the global sti®ness matrix of n Â n, obtained by the FE formulation for linear elastic material; the n-dimensional vector for displacement response; and the n-dimensional vector for outer force, respectively. The subscription ij stands for the case with unit slip in the jth direction in the ith small fault. F ij is computed by inputting the unit fault slip described by a bicubic B-spline shape, obtained using the split node technique [Melosh and Raefsky, 1981]. To compute GFs in this manner, 10 2À3 FE analyses must be carried out. If we use relatively lowresolution FE models, in which the DOF n is about 10 6 , Eq.
(2) is not very di±cult to solve. In fact, some researchers have performed crustal deformation and tsunami computations using this approach. On the other hand, the construction of a highdelity FE model with detailed geometry around the trench axis while simultaneously maintaining the mesh quality is di±cult. Even if such a model is constructed, its DOF is on the order of 10 8 , and the computational cost of solving this large problem (10 2À3 times) is a di±cult issue. Ichimura et al. [2013] proposed an FE model construction method using highdelity crustal structure data and a fast crustal deformation computation method using the FE model. This study extends their method to the construction of a highdelity FE model with the detailed geometry around the trench axis. We also present a method to compute 2n f di®erent G ij rapidly enough for practical use. First, we present a method to construct a high-¯delity three-dimensional (3D) FE crustal model directly from the layered crustal structure data, which includes elevation and material parameters, without using a computer-aided design model or manual tuning. The method is illustrated in Fig. 1. This method generates the mesh in two ways: voxel elements are generated using the background grid, and tetrahedral elements are generated via grid-wise Delauney triangulation. Thus, a robust adaptive mesh can be generated, and control over element quality can be maintained. The method can be used to automatically generate a hybrid mesh with tetrahedral and voxel elements with a good aspect ratio. Finally, a node list sub i is generated for the computation of F ij . In¯nite elements [Zienkiewicz et al., 1983] are generated at the bottom and along the side plane of the constructed FE model to approximate the in¯nite boundary condition. Although solving Eq. (2) 10 2À3 times with a DOF on the order of 10 8 requires a large computational cost, the use of a supercomputer leads to low usability. We seek to e±ciently use a general distributed-memory PC cluster to solve the equations. Two possible approaches to the solution of Eq. (2) in such an environment are available: solving KG ij ¼ F ij 2n f times sequentially using message passing interface (MPI) parallel computation for the entire computational environment, and solving KG ij ¼ F ij independently using each computation node. In this study, we adopt the latter approach to avoid construction of FE models for MPI parallel computation and the degraded performance caused by such computation. Thus, we use a commonly used iterative solver, the conjugate gradient (CG) method, which enables computation in a more memory-e±cient manner compared with direct solvers. In this approach, we should note that the available amount of computational memory is small and that less CPU power may lead to an excessively long computational time. We use the fast and memory-e±cient solver suggested by Ichimura et al. [2013] to avoid these problems. In the solver, the matrix-vector multiplication consisting of K and a temporal n-dimensional vector t require large computational memory if the entire K is stored in the memory. Therefore, we instead replace the multiplication with an equivalent on-the-°y computation: where Kt l e , Kv l e , and Ki l e are an element sti®ness matrix in the lth linear tetrahedral element, the trilinear voxel element, and the in¯nite element, respectively. ne t , ne v , and ne i is the number of linear tetrahedral element, the trilinear voxel element, and the in¯nite element, respectively. The element sti®ness matrices are computed on the Notes: (a) Let a uniform background cell cover the entire target domain. We can set the necessary minimum resolution ds. (b) The geometry of the ground surface and interfaces is simpli¯ed slightly to maintain good element quality. In the target region, we generate tetrahedral and voxel elements using the background grids. Tetrahedral elements are generated in grids that intersect with the ground surface or interface, and these are generated using grid-wise Delauney triangulation. If the voxel elements are unnecessarily small, we enlarge the elements by merging eight small voxel elements into one. (c) We insert tetrahedral elements between the elements of di®erent sizes to maintain the integrity of the mesh. In¯nite elements are generated on the sides and bottom of the model. We pick up the nodes that are included in subfault i (node list sub i is, respectively, generated for n f subfaults). The black points in the¯gure are the nodes in a subfault, assuming that a subfault consists of four nodes.
°y, without storing them in the memory. t l e is part of the vector t that is included in the lth element. The convergence characteristic of the CG solver is improved by e±cient use of the arithmetic space and the solution space; namely, adaptive preconditioning, mixed precision operations, and the geometric multigrid method. In this manner, the computational time required to solve a problem with about 1.6 Â10 8 DOF is reduced by about 1/7 compared with the use of a conventional method with e±cient usage of memory (for information about the details of the methods and the improvement of performance, see Ichimura et al. [2013]). Algorithm 1 shows an algorithm to obtain G ij by solving KG ij ¼ F ij for n f Â 2 times, which is an extension of this fast and memory-e±cient solver.
By using the computed displacement response G ij due to each unit small fault, we compute the temporal history of crustal deformation as where u sf ij ðtÞ is the temporal history of fault slip in the ith unit small fault in the jth direction. By adding the vector of temporal history of crustal deformation to the original sea°oor elevation data, we compute the temporal history of the change in sea°oor geometry. As the temporal history includes not only the vertical, but also the horizontal, displacement of the sea°oor, the e®ect of translational motion of an inclined plane [Tanioka and Satake, 1996] is also considered. Because Kajiura [1963] theoretically proved that short-wavelength components of crustal deformation do not contribute to the generation of tsunami, we used a¯lter proposed by Kajiura [1963] to attenuate the short-wavelength components of crustal deformation.
Algorithm 1 Algorithm for computing G ij . Here, upper and lower are the set of element numbers located over and under the fault surface, respectively, and δut l e is the displacement associated with the lth element.
for i = 1 to nf do for j = 1 to 2 do 1. compute F ij by split node technique δut is set as half of the dislocation in j direction on the nodes of ith subfault (sub i , see Fig. 1 We input the computed crustal deformation data to the nonlinear shallow water equation below and perform hydrodynamic computation of tsunami [Saito and Furumura, 2009;Baba et al., 2015].
where h is the water height from the sea surface at rest, t is time, is the longitude, is the colatitude, R is the Earth's radius, d is the water depth, and the variables M and N are depth-integrated quantities along longitude and latitude lines, respectively. g is the gravitational constant and n is Manning's roughness coe±cient. In this paper, we use the¯nite di®erence scheme with the spherical coordinate system, leapfrog scheme, and staggered grid for the discretizing of Eq.
(2) [Baba et al., 2015]. The code is parallelized by hybrid parallelization of MPI & OpenMP, and it enables the rapid performance of tsunami analyses.

Numerical Experiment
We apply the tsunami computation method described in Sec. 2 to tsunami in realistic problem settings. By comparing the results with those of a conventional method (i.e. a tsunami analysis method using crustal deformation computation in a homogeneous half-space), we present the advantages of the proposed method. Speci¯cally, we carry out tsunami computations in a domain of the Tohoku region of Japan, shown by the square in Fig. 2, using two examples of fault-slip distribution estimated for the 2011 Tohoku-Oki earthquake with di®erent characteristics. Because the purpose here is to investigate the e®ect of introducing the proposed method to problems in realistic settings and we do not aim to¯t the results to observations, calibration of parameters and comparison with observational data are not performed. We¯rst perform crustal deformation computation in the domain of the black square in Fig. 3. A data set of a 3D four-layered crustal structure by Koketsu et al. [2008] is used following Ichimura et al. [2013], who carried out crustal deformation computations in almost the same domain. Table 1 shows the material properties of each layer. In the FE model used in Ichimura et al. [2013], the top surface of the Paci¯c Plate, which corresponds to the fault surface, did not reach the trench axis, i.e. the geometry around the trench axis was not modeled in detail. This study introduces detailed modeling of the geometry around the trench axis that is consistent with the observational data [e.g. Nakamura et al., 2013]. The new geometry is constructed by referring to the results of a sea°oor topography survey conducted in the area and interpolating the original data in the same manner as in past studies (the adjustable-tension continuous-curvature surface-gridding algorithm [Smith and Wessel, 1990] with a tension factor of 0.25 is used similarly as in previous research [e.g. Ide et al., 2010]. Figure 4 shows the 3D FE model constructed using the proposed   To resolve these slip distributions, we divide the fault surface into 180 unit faults and pick up sub i . The displacement response G ij due to the unit fault slip in two directions, the strike-slip direction and the slip direction perpendicular to it, is computed in each small fault. As a result, we carry out 360 crustal deformation computations. Figure 3(c) shows the central points of the unit small faults. Although the proposed method enables the computation of single G ij in a short time, computation of all of G ij is accelerated by performing input/output operations e±ciently with the algorithm of Agata et al. [2016]. With the presented improvements, we computed all of the G ij in 23.7 h using 16 computer nodes of a PC cluster (each node had 12 computer cores and 48 GB computer memory; 16 nodes of IntelR Xeon X5680 [6-core 3.33 GHz/12 MB/QPI-6.4 GT/s]; Â 2 and 8 GB [DDR3-1333 Registered ECC] Â 6). We assume here that the fault rupture propagates from the rupture initiation point in a concentric circular shape at 2562 m/s. The fault rupture velocity is 0.72 times the shear wave velocity [Geller, 1976], and the shear wave velocity is the average of those in the hanging wall and the footwall. The rupture initiation point is located at the central point of the small fault which includes the epicenter of each fault model. The temporal history of sea°oor elevation is computed by superimposing G ij based on Eq. (4). Figures 5 and 6 show snapshots of the temporal history of 1750018-10 J. Earthquake and Tsunami 2017.11. Downloaded from www.worldscientific.com sea°oor elevation. We hereafter call the high-¯delity model used in the abovedescribed analysis the \HFM". For comparison, we also show the results obtained using a simpli¯ed model (SM), in which the surface geometry is°at and the material properties of all layers are the same as those of layer 1. The results obtained with the SM are equivalent to those obtained using an analytical solution in a homogeneous half-space. The settings of the fault surface and the unit small faults are equivalent to those in the HFM. The ground surface elevation is the same as the average of that in the HFM. We use the fault surface geometry from Koketsu et al. [2008] without interpolation in the SM because consistent connection of the fault surface and the trench axis is di±cult. Following the manners in past studies [e.g. Satake et al., 2013], the sea°oor displacement is calculated by mapping the deformation computed in the SM to the sea°oor geometry. We consider four cases consisting of di®erent crustal and fault models, which are listed in Table 2. The SuzukiHFM and SuzukiSM results in Fig. 5 show large di®erences in deformation À À À not only in the¯nal displacement, but also in the deformation process À À À whereas the overall di®erence between the OzawaHFM and OzawaSM results in Fig. 6 is smaller. To see the di®erence   Fig. 2. The horizontal axis is the x coordinate (km) for the region A-A 0 . SuzukiHFM and SuzukiSM produce large di®erences, particularly in the short-wavelength components near the trench axis. OzawaHFM and OzawaSM consist mainly of long-wavelength components, and their di®erence is relatively small. CR, UM and PAC stands for Crust layer, Upper-mantle layer and Paci¯c Plate, respectively, as shown in Table 1.  Fig. 2 in Fig. 7. It is natural that fault slip close to the sea°oor is a®ected largely by the approximation in the crustal model. Increasing the reliability of crustal deformation computations by consideration of shallow slips is also important, as improved observational techniques have revealed that fault slip reaches the shallow fault surface in some earthquakes, e.g. Fujiwara et al. [2011] showed that the fault slip in the 2011 Tohoku-Oki earthquake reached the trench axis by comparing bathymetric data acquired before and after the earthquake. In addition,°attening the surface geometry or decreasing the geometric resolution of the HFM may lead to larger di®erences in the crustal deformation results compared with those obtained with the full HFM because the e®ect caused by such approximation depends on the boundary conditions and the material properties. To make full use of the valuable data obtained by observation, the use of FE models of higher¯delity than the observational data in crustal deformation is important. The tsunami hydrodynamic computations is performed in the domain of 36. 368-40.6678 N and 140.5140-150.0000 E, which includes the focal region and the east coast of the Tohoku region (Fig. 2). The grid size of the¯nite di®erence scheme is 2 s (in the longitudinal and latitudinal directions). About the see°oor and groundsurface topographic data used in the tsunami computations, see Data and Resources. The temporal history of sea°oor elevation is inputted to the tsunami computation, which is computed by Eq. (4) with G ij . A¯lter that attenuates the short-wavelength components of crustal deformation proposed by Kajiura [1963] is used here. The transparent boundary condition is used at the boundaries of the target domain in the o±ng, whereas the inundation of a tsunami is computed at the land-sea boundary, introducing the wet-dry condition. Manning's roughness coe±cient, which is related to the friction of bottom surfaces, is 0.03 m 1=3 s in the terrestrial area and 0.02 m 1=3 s on the sea°oor. The tsunami propagation and inundation process for 6 h is computed with the time-step length of 0.1 s, which satis¯es the stability condition for the numerical scheme (the CFL condition). We used 144 computer nodes of the K computer [Miyazaki et al., 2012], and the execution time was about 4.5 h.
We compare the tsunami analysis results at 43 points, which include 35 points where tsunami observational data are available (see Satake et al. [2013] for detailed information) and 8 points near the tsunami source (Fig. 2). Figure 8 picks up the waveforms at 41, 42, and 43, each of which illustrates the typical characteristics of tsunami response. 41 (near the tsunami source) and 42 (near the coastal area) show a maximum 50% di®erence in amplitude between the SuzukiHFM and SuzukiSM results, and the waveforms also di®er signi¯cantly. In contrast, only small di®erences are seen in the results of the OzawaHFM and OzawaSM. This tendency is common at 23 of 43 points. On the other hand, the di®erences between the SuzukiHFM and SuzukiSM results and between the OzawaHFM and OzawaSM results are small at the other 20 points, as in 43. Actually, it is di±cult to estimate this tendency beforehand due to the complexity of conditions in these computations.  Figure 9 shows the inundated area and the inundation height in a larger domain of the dashed square in Fig. 2. Figure 10 shows the inundated area and the inundation height in a smaller domain of the dashed square in Fig. 9. The di®erences between the OzawaHFM and OzawaSM results are small in the entire domain, but those between the SuzukiHFM and SuzukiSM results for the inundated area and inundation height are signi¯cant in some areas, as in Fig. 10. The source of the di®erence is only the di®erence in inputted crustal deformation. Additionally, as seen in Fig. 7, the Oza-waHFM and OzawaSM, which consist mainly of long-wavelength components, do not show much di®erence because the fault slips occur only in the deep part of the fault surface. Such small di®erences in tsunami analysis results are reasonable when the di®erence in the inputted crustal deformations is small. The short-wavelength components of crustal deformation due to shallow slips produce a large di®erence in the SuzukiHFM and SuzukiSM, particularly near the tip of the trench. As a result, we expect the SuzukiHFM and SuzukiSM results for the short-wavelength components in tsunami analysis to di®er. Thus, we expect that the introduction of the HFM to tsunami analysis can increase reliability in some cases.

Concluding Remarks
We proposed a method enabling tsunami analysis using a high-¯delity crustal structure and geometry model and showed the advantages of this method by comparing the results with those of a conventional method. The proposed and conventional methods produced signi¯cantly di®erent results when fault slip reaches the shallow fault surface, but they show smaller di®erences for deep fault slip. We can expect more tsunami data and crustal structure data to become available in the near future because of the advancement of observational techniques. In such a situation, the proposed tsunami computational method, which considers the e®ects of high-  Fig. 9. The di®erence in the crustal structure models, HFM and SM, solely produces di®erences in di®raction, which lead to large local di®erences in inundation height. ¯delity crustal structure and geometry derived from observations, is expected to become increasingly important. Because this study aimed to develop the method and compare its results with those of conventional methods, we considered only typical problem settings and did not compare the results with tsunami observational data. In future work, we would like to compare the results with observational data and investigate the e®ect of introducing the proposed method using an updated highdelity model with soft sedimentary layers on the sea°oor.

Data and Resources
In tsunami computations, data from the M-7000 series of the Japan Hydrographic Association (http:==www.jha.or.jp=shop=index.php?main page=advanced search result&categories id=1308), General Bathymetric Chart of the Oceans (GEBCO) by British Oceanographic Data Centre (http:==www.gebco.net=data and products= gridded bathymetry data), and submarine topographic surveys conducted by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) [Kido et al., 2011] are used to construct the sea°oor tomography data in the target domain. M-7000 is a digital data set of sea°oor topography that covers the domain from the coastal area to 70 miles o® the coast. It contains information on the locations of the isobathic lines, and the spatial resolution is 50 m in the most detailed part, near the coastal area. GEBCO is a topography mesh data set with a 30-s grid size that covers terrestrial and the oceanic areas around the world. The submarine topographic survey data of JAMSTEC has high resolution in the deep sea, but it does not cover the entire domain of analysis. We construct smooth sea°oor topography data with a 2-s grid size by using mainly M-7000 and JAMSTEC survey data for the coastal area and the deep sea, respectively, while interpolating GEBCO data for the region without those detailed data. We use 50-m grid data from the Geospatial Information Authority of Japan (GSI) for the topography of terrestrial areas, and we assume the locations in which GSI data exist to be terrestrial, i.e., the coastline is de¯ned by the GSI data.