High-fidelity elastic Green’s functions for subduction zone models consistent with the global standard geodetic reference system

Green’s functions (GFs) for elastic deformation due to unit slip on the fault plane comprise an essential tool for estimating earthquake rupture and underground preparation processes. These estimation results are often applied to generate important information for public such as seismic and tsunami hazard assessments. So, it is important to minimize the distortion of the estimation results on the numerical models used for calculating GFs to guarantee assessment reliability. For this purpose, we here calculated GFs based on a numerical model that is of high fidelity to obtain realistic topography and subsurface structural models of the Earth. We targeted two well-known subduction zones in Japan, the Nankai Trough and the Japan Trench. For these subduction zones, databases for realistic topography and subsurface structural models of the Earth are available in the “Japan integrated velocity structure model version 1”, which was proposed for earthquake hazard assessments conducted by the Japanese government. Furthermore, we eliminated the inconsistency in processing calculated GFs and space geodetic observation data for surface displacements, which is often overlooked, by using the same coordinate system. The ellipsoidal shape of the Earth, which is often approximated with a projected plane or a spherical shape, was also incorporated by faithfully following the definitions of the coordinate systems in Geodetic Reference System 1980, which is the global standard for space geodesy. To calculate elastic GFs based on such high-fidelity subduction zone databases with the ellipsoidal shape of the Earth, we introduced the finite element (FE) method. In the FE meshes, the resolution of the topography and subsurface structure is the same as that of the original databases. Recent development of the state-of-the-art computation techniques for the rapid calculation of crustal deformation using large-scale FE models allows for GF calculation based on such a high-fidelity model. However, it is generally not easy to perform such calculations. Thus, we released a library for the GFs calculated with 1-km grid spacing on the ground surface in this study to the geoscience community on a web server, aiming to contribute more reliable seismic and tsunami hazard assessment.


Introduction
Accurately estimating earthquake rupture and the preparation processes underground, including aseismic slips and interplate couplings, is essential for better understanding earthquake generation processes and assessing seismic and tsunami hazards. However, because direct Open Access † Takane Hori and Ryoichiro Agata contributed equally to this study *Correspondence: horit@jamstec.go.jp 1 Research Institute of Marine Geodynamics, Japan Agency for Marine-Earth Science and Technology, Yokohama, Kanagawa, Japan Full list of author information is available at the end of the article measurements of the plate interface underground where the slips take place are limited to a few examples in the shallowest portion from the seafloor by ocean drillings (e.g., Tobin et al. 2009), surface measurements are the only observation available for such estimations in most cases. In general, by assuming linear elastic media, surface displacements due to arbitrary fault slips in a known fault plane embedded underground are calculated by superimposing surface response (Green's) functions (GFs) for elastic deformation due to unit slip on the fault plane. This allows for linear slip inversions of geodetically observed surface displacements and displacement rates to estimate coseismic slips (e.g., Yabuki and Matsu'ura 1992;Chlieh et al. 2007;Iinuma et al. 2012), aseismic slips (e.g., McGuire and Segall 2003;Miyazaki et al. 2004;Nakata et al. 2017;Qiu et al. 2018), and slip deficit rates (e.g., Jin et al. 2007;Hashimoto et al. 2009;Loveless and Meade 2010;Yokota et al. 2016). Tsunami wave form inversion of coseismic slip also uses GFs to relate seafloor deformation to fault slip distribution (e.g., Satake 1987). Thus, a precomputed set of GFs is a useful tool for investigating underground processes regarding earthquakes based on observations on the Earth surface as well as in other fields of geophysics (Pan 2019).
GFs calculated in the elastic homogeneous half-space have been widely used for such purposes because their analytical expressions (e.g., Comninou and Dundurs 1975;Okada 1985Okada , 1992Meade 2007) enable to calculate them using accessible computation programs without requiring high computation costs. Recent geodetic observation techniques, such as continuous nationwide network of Global Navigation Satellite System (GNSS) observation in land areas (e.g., Miyazaki and Hatanaka 1998;Beavan et al. 2016), high-resolution land deformation data provided by interferometric synthetic aperture radar (InSAR) (e.g., Massonnet et al. 1993;Wright et al. 1991;Kobayashi 2017), seafloor deformation data from GNSS-Acoustic (GNSS-A) technique (e.g., Chadwell and Spiess 2008;Sato et al. 2011;Tomita et al. 2015;Yokota et al. 2016) and ocean bottom pressure (OBP) gauges (e.g., Ito et al. 2011;Kaneda et al. 2015), have provided more collective information on crustal deformation. Combined use of these observation data can result in improved analyses of crustal deformation for large regions where topography and elastic properties have significant spatial variations. Notably, several studies reported that a consistent analysis of geodetic data distributed in such a wide region often requires a more realistic Earth model than the elastic homogeneous half-space, by incorporating realistic Earth models based on the finite element (FE) (Masterlark 2003;Aagaard et al. 2013;Ichimura et al. 2016) or spectral element (SE) methods (Komatitsch and Tromp 2002;Gharti et al. 2018). For instance, the effects of topography and elastic properties were discussed in Masterlark (2003), Ichimura et al. (2013), Kyriakopoulos et al. (2013), Williams and Wallace (2015), Williams and Wallace (2018), Langer et al. (2019). These results suggest that GFs calculated by incorporating realistic topography and subsurface structure models are indispensable for estimating fault slips and interplate couplings based on recent geodetic data, which are observed not only on land but also on the ocean floor.
Neglect of the ellipsoidal shape of the Earth, usually approximated with a projected plane or a spherical shape at best, is another factor that introduces errors when comparing calculations to observation data. This approach may be problematic in two aspects. One is the effect of the ellipsoidal shape in calculating the responses due to fault slips, which has been studied by Cheng et al. (2019). The other is disagreement between the coordinate system used to describe observed data and that for the utilized numerical model: for GNSS observation data, a coordinate system defined in Geodetic Reference System 1980 (GRS80) adopted by the International Union of Geodesy and Geophysics (Moritz 2000) and the Japanese Geodetic Datum 2000 (Matsumura et al. 2004) following the International Terrestrial Reference Frame (http:// itrf. ensg. ign. fr/) is usually applied, while the numerical models are usually defined in a locally projected Cartesian or a symmetric spherical coordinate system. Although less attention has been paid to this latter aspect, it is a more fundamental issue that includes the former. If the modeled area is small, no significant distortion will occur, regardless of which coordinate transformation method is used. However, as mentioned above, when a relatively wide area where the ground is considerably curved is targeted, the difference between the coordinate systems may result in a non-negligible error in the calculated response. Such errors can be avoided when the coordinate system used for the model calculation is the same as that for GNSS observation data. Therefore, to calculate GFs for a significantly wide region, it is preferred to consider the Earth as an ellipsoidal-shaped body by following the definitions of the coordinate systems in GRS80. Hereafter, we use the term "high-fidelity model" to refer to a numerical model that faithfully incorporates the database for realistic topography and subsurface structure models of the Earth and the coordinate system used in the geodetic reference systems.
The use of high-fidelity FE models in calculating GFs is not an easy task in terms of computation cost. Recent development of the state-of-the-art computation techniques highly tuned-up for massive parallel computation environments Fujita et al. 2017) can provide a solution to overcome such calculation cost. However, despite the wide availability of environments in which such tools have been well developed in terms of both software and hardware, the human resources (such as training and labor availability) required to use such tools is still not widely available. Therefore, publishing a set of GFs calculated based on a high-fidelity numerical simulation model is a significant contribution to the seismological research community.
Consequently, we developed an elastic GF library based on three-dimensional heterogeneous subduction zone models using the FE method, which is open to the public on a web server. The target regions here are the Nankai Trough and the Japan Trench subduction zones, both of which are known for megathrust earthquakes (Ando 1975;Ozawa et al. 2011). For the Japanese Islands, a unified velocity structure model "Japan integrated velocity structure model version 1" (JIVSM) (Koketsu et al. 2009(Koketsu et al. , 2012, which has been proposed for earthquake hazard assessments conducted by the Japanese government, was deemed suitable for our study. We used a highfidelity FE model that faithfully follows the resolution and layer boundaries in the subsurface structure model in JIVSM. The coordinate transformation required for processing the data faithfully follows definitions of the coordinate system used in the geodetic reference systems. After the coordinate transformation, we constructed a high-resolution FE mesh from the structured data using an automated mesh generator developed particularly for stratified layered structures (Ichimura et al. 2009. Then we perform FE analyses using a high-performance FE solver (Fujita et al. 2017), which allow for calculating GFs incorporating the high-resolution digital elevation model (DEM) data of JIVSM.
In the following sections, we first describe the procedure of constructing FE models from a seismic velocity structure data set defined in the geographical coordinate system. We also explain the calculation method of GFs using advanced computation techniques for FE modeling. Then we describe the details of calculating the GF library in the Nankai Trough and Japan Trench subduction zones. An application example of the GF library in geodetic inversion of the interseismic displacement rates in the Nankai Trough subduction zone is also presented.

Coordinate transformation between the geographic and a Cartesian coordinate system
Coordinate transformation of underground structure data from the geographic coordinate system to a local Cartesian system To calculate elastic GFs, we must assume the elastic properties of the underground structure. JIVSM, as with other seismic velocity structure models that are defined for a wide region including the focal region of a magnitude-9 class earthquake, such as CRUST1.0 (Laske et al. 2013), includes DEM data consisting of multiple layers defined in the geographical coordinate system (i.e., longitude, latitude, and elevation) with seismic wave velocities in each layer. In contrast, calculating elastic deformation is usually performed based on a Cartesian coordinate system. Therefore, coordinate transformation of DEM data from the geographical coordinate system into a Cartesian system is necessary for constructing a numerical simulation model (i.e., an FE model). However, the projection to half-space (e.g., Mercator Projection) or approximation with sphericity, in which some user-dependent parameters are required to define the transformation, are usually performed to construct a numerical simulation model in a Cartesian coordinate system for GF calculation. Acknowledging these considerations, we describe the procedure of DEM data coordinate transformation in this section, faithfully following the genuine definition of the coordinate systems.
As mentioned in Introduction, we performed coordinate transformation based on GRS80. The reference ellipsoid of GRS80 is practically equivalent to that of the World Geodetic System 1984 (WGS84) (https:// www. nga. mil/ Produ ctsSe rvices/ Geode syand Geoph ysics/ Pages/ World Geode ticSy stem. aspx), (Slater and Malys 1998), which is used for data analyses in the global positioning system of the United States. We defined a target Cartesian coordinate system in which the origin is located at a representative (central) point of the region of interest P c . The north direction is the y-axis, and the upward direction is the z-axis (then, the x-axis is automatically determined in the right-handed system). Hereafter, the unit of length is meter (m) unless otherwise noted. We first transformed the coordinate of a point P in the input DEM data to the Earth-centered, Earth-fixed (ECEF) Coordinate System (Fig. 1a), as where ϕ , , and h are latitude, longitude, and ellipsoidal height (the sum of elevation and geoid height), and X, Y, and Z are the coordinates in the ECEF system. N(= a/ 1 − e 2 sin 2 ϕ ), a, and e are the prime vertical radius of curvature, equatorial radius, and first eccentricity of the Earth ellipsoid. For the reference ellipsoid of GRS80, a = 6378137 m, e = 0.081819191042815791 . We use the geoid model in the Earth Gravitational Model 2008 (EGM2008) (Pavlis et al. 2012), a high-resolution global geoid model including both the land and sea area. (1) We then translated and rotated the coordinate in the ECEF system to a local Cartesian system in the region of interest ( Fig. 1b) for convenience of defining the target region for FE mesh generation, as where ϕ c , c , and h c are the latitude, longitude, and ellipsoidal height of P c , X c , Y c and Z c are the coordinates of P c in the ECEF system, and x, y, and z are the coordinates in the local Cartesian coordinate system. The mapping matrix in Eq. 4 is the product of a rotation matrix with π 2 + c around the Z-axis, one with π 2 − ϕ c around the X-axis. In this manner, we performed coordinate transformation of each point in the original DEM data set defined in the geographical coordinate system. The method we adopted to construct FE models requires DEM data whose points are aligned in an equally spaced structured grid on the x-y plane in the local Cartesian coordinate system; this is described in detail in the next subsection. Although the data points of the original DEM data are usually defined in an equally spaced structure grids in the geographical coordinate system, the structured grid does not hold after the coordinate of the data points are transformed into the local Cartesian system based on Eq. 4. Therefore, we generated new data points that are aligned in an equally spaced structured grid in the x-y plane of the local Cartesian system, whose z coordinate we interpolate from the neighboring points in the original data (Fig. 1c). We use the nearneighbor interpolation in which an average value of the nearest four points, which is weighted by the inverse of the distance to each of them, is assigned to each new point.

Coordinate transformation of calculated displacements at observation points
From the response displacement field computed using the FE method (see the next section), we can extract displacements at any points on the ground surface. While the xyz components of the displacements are defined in the local Cartesian coordinate system characterized by P c , geodetic observation data of GNSS and GNSS-A are usually provided as ENU (east-west, north-south, and up-down) coordinates. Therefore, we cannot directly compare the observed displacements with the calculated ones. The displacements calculated Fig. 1 Schematics of the coordinate transformation. For simplicity, 1/8 portion of the GRS80 reference ellipsoid is visualized. a Transformation of the coordinate of P from the geographical coordinate system to the ECEF system. Notably, P is located at the ellipsoidal height h. Because the region for the input data is defined in a rectangle area in the geographical coordinate system, it appears to be distorted in the ECEF system, as shown by the area surrounded by the dotted-dashed lines. b Transformation of the coordinate of P from the ECEF system to a local Cartesian coordinate system in which the origin is located at P c . Assuming h c = 0 , P c is located on the surface of the GRS80 reference ellipsoid. c Generation of new data points denoted by the black circles that were aligned in an equally spaced structured grid in the x-y plane in the local Cartesian coordinate system. White circles are the data points in the original DEM data that are not aligned in an equally spaced structured grid in the local Cartesian system. The z components in the points of the black circles are interpolated from those of the white circles by using the near-neighbor algorithm, as denoted by arrows in the local Cartesian coordinate system are required to be transformed into ENU components at P. In a consistent manner with that presented above, we perform the coordinate transformation of the calculated displacements based on the GRS80 ellipsoid. The coordinate transformation is defined as where u E , u N , u U , u x , u y , and u z are displacements in the ENU and xyz components at P. This manner first transforms the components in the ECEF system and transform it again to the ENU components in P. The surface tilts and strains, which are also calculated for GFLSZ (see Appendix), are transformed into the ENU system in the same manner.

FE mesh construction
In developing the GF library, we aimed to use a numerical simulation model of high fidelity to the reference underground structure database, i.e., FE mesh with the layer-boundaries in the same resolution as that of JIVSM, which is approximately 1 km. Suppose that we construct such a mesh in a large domain in which the length of one side is a few thousand kilometers with more than 10 layers of different seismic velocities, the typical degreeof-freedom (DOF) is expected to be in the order of billions. Robust construction of high-quality mesh with complex shapes of the layer-boundaries in such a largescale model is not an easy task. In general, producing the FE mesh itself sometimes takes much more time than FE analyses for large-size problems (Fujita et al. 2016). We constructed a high-resolution FE mesh from the DEM data after the coordinate transformation using a robust, fast, and automated tetrahedral mesh generator developed specifically for stratified layered structures , which is an improved version of the method proposed by Ichimura et al. (2009). In the method, the mesh is generated based on a uniform background cell covering the entire targeted domain, which defines the resolution of the layer interfaces ds. The geometries of the ground surface and the layer interfaces were simplified slightly to maintain good element quality. Unnecessary elements were merged to generate larger elements elsewhere. This method enabled the automated and robust construction of high-resolution tetrahedral mesh directly from digital elevation model (DEM) data of crustal structure without creating a computer-aided design (CAD) model.

Calculation of GF due to unit slip on a subfault based on FE modeling
For the basis function of unit slip, we used the normalized third-order bicubic B-splines (Yabuki and Matsu'ura 1992) (please see De Boor 1972 for the definition). We input B-spline-shaped slips on subfaults located on the plate boundary along a structured grid (also see Fig. 5 in the next section). We considered slip in two directions: that of subduction set based on past studies (Heki and Miyazaki 2001;Sella et al. 2002) and that perpendicular to it. The slip in the location beyond the trench or trough axis was set to zero (see also Fig. 5a in the next section). We used the FE method to compute the response displacement due to unit slip on the subfault. By applying FE formulation to the governing equation of elastic deformation and arranging it with boundary conditions, the problem becomes the following linear equation: where K , u , and f are the global stiffness matrix, the displacement vector, and the outer force vector. Fault slips are converted to the force vector using the split node technique (Melosh and Raefsky 1981). The Dirichlet condition is applied on the side and bottom of the modeling domain. Specifically, displacement in the direction of the normal vector of the surface was fixed at zero, and there is no imposed displacement in the other directions (Aagaard et al. 2013). As aforementioned, an FE model with high fidelity of the data of the seismic velocity structure typically results in billions of DOF. Therefore, we must solve a large linear system to obtain u . For this purpose, we introduced a state-of-the-art FE solver for computing crustal deformation (Fujita et al. 2017). The basic algorithm consists of the conjugate gradient (CG) method with memory-efficient calculations of matrixvector products. The convergence of the CG solver is improved by using an adaptive pre-conditioner that combines the geometric and arithmetic multigrid method. As (6) Ku = f , a result, the solver obtained a speedup of 20 times from a widely used solver in some cases, being scalable for 10 5 of computation cores (Fujita et al. 2017).

Calculation of Green's functions in the Japanese Islands
Overview of "Japan integrated velocity structure model version 1 (JIVSM)" Following the method described in the previous sections, we constructed a numerical simulation model of the Nankai Trough and the Japan Trench subduction zones based on JIVSM. JIVSM is a seismic velocity structure model that covers a great portion of the Japanese Islands, which was developed to simulate long-period ground motion and its seismic hazard assessment. A structure including shallow soils to oceanic and continental mantles is modeled as a 23-layered structure, whose geometry and elastic properties are provided in the data set. Because the shallowest 14 layers in the model are thin, and seismic-velocity variation among these layers is assumed to be insensitive to crustal deformation at a relatively long wave length, we merged the 14 layers into two for our calculation: one that includes the accretionary wedge and layers shallower than that and another that includes the remaining layers. Because the JIVSM is basically defined in the old Japan Geodetic System (Tokyo Datum) based on the Bessel ellipsoid, we first converted the coordinates to GRS80 and then applied them to the aforementioned procedure.
We consider the deep portion of the plate boundary, where interplate fault slips are expected to take place, is located at the boundary between the top of the subducting oceanic crust and the bottom of the continental crust or mantle in JIVSM. However, in the shallow portion of megathrusts, fault slips do not always take place in the boundary between the continental and oceanic crust, nor are they necessarily to occur on a single plane. Because this study does not aim to address such complexity for calculating GFs, we simply assume that the relative motion between two plates is accommodated on a single fault. We generated a new surface to apply shallow fault slips by connecting the trough or trench axis and the boundary between the continental and oceanic crust smoothly. The connecting points on the boundary of crusts are located at the depth of 15 km, which allows for the generation of a smooth surface. Following studies that generated the surface using point constraints (Ide et al. 2010), we generated the surface using an adjustabletension continuous curvature surface gridding algorithm (Smith and Wessel 1990). The newly generated curved surfaces are shown in Fig. 2.
The parameters for mesh construction, such as the region size and mesh resolution, were chosen based on the result of evaluating numerical solutions of elastic deformation due to an M w 9 earthquake in the Japan Trench region (Agata et al. 2015). As a result, the element size in the horizontal directions in the layer interfaces is ds = 1 km (ds corresponds to the side length here) and the region size is 2500 km × 2500 km × 1100 km. The model domain is taken to be large enough that the approximation error due to the boundary condition on the model boundaries negligibly affects the calculation of displacements in the observation points in the area of interest, i.e., those located in the area of −1000 km ≤ x ≤ 1000 km and −1000 km ≤ y ≤ 1000 km with P c as a b Fig. 2 Seismic velocity (P-wave) profiles on cutting planes. The profiles on the cutting plane along the red lines in Fig. 4 for the a Nankai Trough and b Japan Trench subduction zones are shown. The above panels are the overviews. The bottom panels are close-up views in the areas around the trough and trench axes denoted by the red rectangles in those above. The purple lines denote the newly generated curved surfaces by interpolating between the top surface of the subducting plate and the trough or trench axis the origin. The vertical thickness of each layer is set to be zero in the FE model if it is less than 0.125 ds in the original data because existence of such thin layers leads to difficulties in FE mesh generation. Based on the results of some preliminary calculations with different element size ds, we further refined the elements in the sediment layers in the hanging wall to ds = 0.5 km.

Calculation of GFs in Nanakai Trough and Japan Trench
We calculated GFs in the Nankai Trough subduction zone in southwest Japan (Fig. 3a). We located P c in 80 km off Cape Muroto (Table 1). Figure 4a shows the FE mesh for the target region. the meshed model surface well expresses the ellipsoidal shape of the earth centered on the east of Cape Muroto. Figure 2a shows the cross section of the velocity structure incorporated in the FE model. On the plate boundary of the constructed model, we located 460 structured grid points to introduce subfaults with 3rd-order B-spline function of 20-km node spacing, where we set the depth limit of slip as 50 km. Considering two directions of slip on each subfault, we a b Fig. 3 Overview of the target regions for calculating GFs. a the Nankai Trough and b the Japan Trench subduction zones. Gray dots denote the central point of B-spline-function-shaped slip distribution applied in a subfault. Triangles and inverse triangles denote the GEONET and the representative seafloor stations (Kaneda et al. 2015;Sato et al. 2011;Yokota et al. 2016) of geodetic observation, respectively Table 1 Location of P c , the area for calculating GFs and the subduction direction in the Nankai Trough and the Japan Trench subduction zone, respectively The subduction directions are set based on Heki and Miyazaki (2001) for the Nankai Trough and (Sella et al. 2002) for the Japan Trench, respectively calculated 920 GFs in total. Figure 5a shows an example of calculated GFs (surface displacements) due to a B-spline-shaped slip applied on a subfault on the plate boundary in the Nankai Trough subduction zone. Note that unit slip in the location beyond the trench or trough axis is set forcibly to zero in this example.

Region for GFs
In the same manner, we also calculated GFs in the Japan Trench subduction zone in northeast Japan (Fig. 3b). We located P c near to the epicenter of the 2011 Tohoku-Oki earthquake (Table 1). Figure 4b shows that the meshed model surface well expresses the ellipsoidal shape of the earth centered on the epicenter. Figure 2b shows the cross section of the velocity structure incorporated in the FE model. In the constructed model, we located 405 structured grid points with 20-km spacing, setting the depth limit of slip as 100 km. In total, we calculated 810 GFs. Figure 5b shows an example of the calculated displacements (surface displacements) due to a B-spline-shaped slip applied on a subfault on the plate boundary in the Japan Trench subduction zone.
See Appendix for the information on the file contents of the GF library.

Application to slip deficit rate estimation in the Nankai Trough
Using the GF library, we estimated slip deficit rate (SDR) during the interseismic period in the Nankai Trough subduction zone and compared the result with the one obtained based on a homogeneous elastic half-space Yokota et al. (2016). Yokota et al. (2016) estimated interplate SDR of the Nankai Trough subduction zone using the interseismic displacement rate data: seafloor geodetic data of GNSS-A and continuous land data of GNSS. They performed geodetic slip inversion of these data using the method of Yabuki and Matsu'ura (1992) Examples of GF calculated based on FE modeling. a A unit slip and its response displacement in the direction opposite to the plate convergence in Nankai Trough and b one in the direction perpendicular to the plate convergence in the Japan Trench subduction zone. Note that unit slip in the location beyond the trench or trough axis is set forcibly zero in a. Gray dots denote the central point of B-spline function shaped slip applied in a subfault. The color map shows slip distribution. The black arrows denote the horizontal (top) and vertical (bottom) surface displacements due to the fault slip. The white arrows denote fault slip directions as prior information and the hyperparameter on the constraint was determined using Akaike's Bayes information criterion (ABIC). We used the same data set and method to estimate the interplate SDR of the Nankai Trough subduction zone, except that we incorporated 3D heterogeneous elastic property and the ellipsoidal shape of the Earth in the analysis using the GF library. Figure 6a shows the estimated SDR using the GF library. Large SDR is estimated from Cape Ashizuri, which is similar to the preferred model of Yokota et al. (2016) (Fig. 7a) and other SDR models (Watanabe et al. 2018;Noda et al. 2018). However, the area of the large SDR is significantly smaller in our model than in the model of (Yokota et al. 2016), which appears to be more consistent with the fact that the maximum SDR in this subduction zone is expected to be no more than the reported plate convergence rate, i.e., around 7 cm/ year (Heki and Miyazaki 2001). Particularly, our model effectively explains the observed displacement at the station ASZ2 without requiring large slip deficit on the plate boundary right below the station as seen in the preferred model of Yokota et al. (2016). Comparing Figs. 6b and 7b, we also observed better fit to the data in the station TOK1 and TOK2 in our SDR model. This difference in the fit to the observation may be attributed to the consideration of consistency between the fault geometry, the trough axis, and surface geometry in our numerical model. Such consistency is expected to significantly increase the accuracy of the model prediction if the database for the underground structure incorporated in the analysis is enough reliable. Furthermore, this consistency in modeling the structure around the trough axis allows for locating fault patch close to the trough axis. This contrasts with the conventional approaches based on the combination of a planar fault and flat ground surface, which does not allow for locating such patches when the curved trough axis and the surface topography are to be considered. It is natural that a fault slip model which considers such fault patches better explains the observation close to the trough axis.
In summary, introducing a more realistic Earth model gives a more consistent estimation result to both the knowledge on plate tectonics and the geodetic data. Further discussion on the SDR distributions should be made after considering other effects such as in-land rigid block motion (Nishimura et al. 2018;Noda et al. 2018) and viscoelastic relaxation in asthenospheric mantle (Noda et al. 2018;Johnson and Tebo 2018;Watanabe et al. 2018), which are beyond the scope of this study.

Conclusion
We propose the concept of matching the coordinate system and geometry of the calculation model and the observation data when calculating the GF. In other words, the world geodetic system (ITRF and compliant ellipsoid GRS80) used in GNSS for geodetic data observation was introduced. In addition, three-dimensional topography and underground heterogeneous structure, which has recently been reported as necessary, have been introduced. Because the calculation is expensive, we published the results of the calculation using the large-scale fast finite element algorithm on the web as a library. We expect that our library will help researchers incorporate a more realistic Earth model in their analyses without laborious numerical model construction, numerical simulations, and data processing.
We used FE models with high fidelity to a unified velocity structure model in Japan; the resolution of the layer boundaries of the underground structure is the same as that of the DEM data in "Japan integrated velocity structure model version 1" (Koketsu et al. 2009(Koketsu et al. , 2012. The coordinate transformation of the DEM data of the velocity structure from the geographical coordinate system to a local Cartesian coordinate system was performed faithfully following the definition of each coordinate system, instead of performing projection to half-space or approximation with a spherical shape. Computation cost for FE mesh generation and FE analyses associated with the high-resolution FE models was overcome by introducing a robust, fast, and automated tetrahedral mesh generator developed particularly for stratified layered structures ) and a state-of-the-art FE solver (Fujita et al. 2017) to the FE analyses, respectively. We also presented a demonstrative application to SDR estimation in the interseismic period in the Nankai Trough subduction zone using the proposed GF library. Comparing the estimation result with that based on the simplified Earth model indicated that the proposed GF library provides a more consistent estimation result to both the geodetic data and the understanding of plate tectonics. More detailed discussion on the estimated slip will be accomplished in future work.
We expect that the prevalence of the GF library, which is calculated using numerical simulation models with high fidelity to a database of underground structures and definitions of the widely accepted coordinate systems, helps reduce the inconsistencies in various studies of estimating the underground processes. Continuous effort for developing such tools is expected to eventually contribute to more reliable seismic and tsunami hazard assessments conducted by using those estimation results.
DEM data of the underground structure in the other subduction zones around the world, such as CRUST1.0 (Laske et al. 2013), are also applicable to the procedure of FE mesh generation described in this study, as long as the data sets are defined in the geographical coordinate system.
Our GF library should be revised in accordance with major updates in the reference seismic velocity structure model in the future. In addition, GF libraries for other subduction zones in Japan, such as the Kuril subduction zone, in which The Headquarters for Earthquake Research Promotion in Japan Government forecasted the high probability of a magnitude-9 earthquake (The Headquarters for Earthquake Research Promotion 2017), and the Boso region in east of Tokyo, which is associated with regularly occurring long-term slow slip events (Sagiya 2004), are also expected to be calculated in the next step. These calculations can be performed based on the same framework presented in this study. There is also a demand for viscoelastic GF to analyze geodetic observation data over longer periods. The computation cost required to produce a viscoelastic GF library is expected to be significantly higher than those for elastic one. Therefore, introducing more efficient algorithms of FE analyses for a newer computational environment (e.g., Ichimura et al. 2019) are essential.