Impact of Ambiguity of Physical Properties of Three-Dimensional Crustal Structure Model on Coseismic Slip and Interseismic Slip Deﬁcit in the Nankai Trough Region

: Since huge earthquakes are expected along plate subduction zones such as the Japan Trench and Nankai Trough, the estimation of coseismic slip and interseismic slip deﬁcit is essential for immediate response and preliminary measures to reduce damage. Recently, analysis considering the complex topography and underground structure of the plate subduction zone has been performed for improving the estimation performance. However, the three-dimensional (3D) crustal structural model needs to be improved continuously. In this paper, we obtained Green’s functions for 3D crustal structural models with ambiguity by 3D crustal deformation analysis, and the coseismic slip and interseismic slip deﬁcit were estimated. Here we enabled the calculation of many Green’s functions with different physical properties of the 3D crustal structure by utilizing a GPU-based 3D crustal deformation analysis method that signiﬁcantly reduces the analysis cost. The physical properties on the upper plate’s side, which are located above the plate boundary fault, were changed. We found no signiﬁcant difference in the estimation performance, except for the upper crust, which most of the fault slip area is in contact with, in the case of coseismic slip estimation. In contrast, the coseismic slip estimation when the properties of the upper crust was changed had a signiﬁcant error, and a negative slip was estimated at the deep part of the plate boundary where no slip was originally given.


Introduction
Many large earthquakes are expected to occur along plate subduction zones, such as Cascadia, Alaska, Chile, Kurill, Japan, and Nankai. When such a large earthquake occurs, the damage is expected to be enormous due to shaking and tsunamis. This damage can be reduced by grasping the phenomenon immediately after the earthquake and providing the immediate and efficient evacuation and emergency response. Important proactive measures can also be implemented, such as considering the rupture scenarios of a massive earthquake before the occurrence, quantitatively assuming the ground motion, crustal deformation and the tsunami caused by those scenarios, and identifying and preparing the vulnerable parts of the city. For the former, estimating the distribution and amount of fault slip is necessary when a large earthquake occurs. For the latter, it is necessary to improve the estimation performances of interplate couplings during the interseismic period. Recently, to improve the estimation performance, analysis using crustal deformation observation data on the seafloor just above the epicenter area of the giant earthquake along the trench is performed [1], and analysis considering the complex topography and underground structure of the plate subduction zone has been performed (e.g., [2]). These results suggest that the effects of realistic topography and subsurface structure cannot be ignored in the slip estimation. Precisely, in the Nankai Trough area, in addition to the on land observation network, an advanced seafloor observation network has been developed [3], and plans for further expansion are underway. Additionally, a three-dimensional (3D) crustal structure model to assume long-period seismic motion has been developed [4,5], and has been used to compute the Green's function of crustal deformation [6]. However, the crustal structure model needs to be improved in the future, as shown by the fact that the underground structural exploration in the Nankai Trough region has been performed even after the 3D crustal structure model was constructed. In fact, a 3D crustal structure model focused on the offshore has recently been proposed [7]. Therefore, it is essential for the purposes mentioned at the beginning to understand how the model error affects the estimation results quantitatively. However, since the analysis cost of 3D crustal deformation analysis using 3D crustal structure is high, quantitative evaluation of coseismic slip and interplate coupling considering the error of physical properties and the shape of the crustal structure has not been conducted extensively. Few examples have been conducted by our research group, which evaluated the effects of errors on the physical properties and plate boundary shapes on the slip of the Tohoku-oki Earthquake in the Japan Trench area [8,9]. These studies have also shown that ambiguities in the crustal structure model can affect the estimation results, but only a simplified assessment has been made. Therefore, we have developed a crustal deformation analysis method to perform inverse analysis using a highly detailed 3D crustal structure model, which is essential for the goal of the accurate estimation of highly detailed slip distribution and slip deficit. We integrate GPU-based 3D crustal deformation analysis based on high-performance computation with practical inverse analysis to address more realistic issues which were unexamined in previous studies, such as evaluating which part of the crustal structure mainly affects the estimation slip performance. Specifically, many Green's functions are obtained by 3D crustal deformation analysis using the proposed 3D crustal structure, and coseismic slip and interseismic slip deficit were estimated considering practical observation errors. Furthermore, many 3D crustal structure models with different physical properties were constructed to consider the ambiguity of the physical properties of the 3D crustal structure, and the effect of ambiguity on the estimation performance of the coseismic slip and the interplate coupling was evaluated. Section 2 shows the 3D crustal deformation analysis and estimation methods of the coseismic slip and interseismic slip deficit, and Section 3 shows the estimation performance evaluated by numerical experiments. Section 4 summarizes the paper.

Methods
Herein, we propose an inverse analysis method for evaluating the coseismic slip and interseismic slip deficit on the plate boundary fault using many Green's functions calculated using a GPU-based 3D crustal deformation analysis method that reduces the analysis cost caused by large-scale 3D crustal deformation analysis. By using the proposed method, we evaluate the estimation performance. The details of the GPU-based 3D crustal deformation analysis method and the inversion analysis method for estimation of the slip and slip deficit are described in this section.

GPU-Based 3D Crustal Deformation Analysis
Herein, to accurately model the complex crustal structure and accurately consider the stress-free boundary condition at the ground surface, Green's functions are calculated using the finite element (FE) method with unstructured tetrahedral quadratic elements. Here, the crustal deformation of coseismic slip and interseismic slip deficit due to fault slip is regarded as a linear elastic deformation, and analysis is performed based on the FE method. Using the summation convention, the governing equation is given as: where, Here, σ and f represent the stress and external force. λ, , δ, µ are the first Lame coefficient, strain, Kronecker delta and shear modulus, respectively. By discretizing the governing equation using the FE method, a linear system of equations, is obtained. Fault slip is evaluated using the split-node technique [10]. Dirichlet boundary conditions are imposed on the bottom and sides of the FE model. Generally, it is difficult to construct a large and accurate 3D FE model that accurately reflects a complex crustal structure. Herein, we used the method reported in Ref. [11] to automatically and robustly generate 3D FE models using unstructured tetrahedral quadratic elements. The mesh around the fault, where localization of deformation is expected to occur, is further refined to generate fine elements without a significant increase in the total FE model size. The computational cost of the FE analysis for large-scale models becomes large. The method reported in Ref. [12] has been proposed for reducing the computational cost of large-scale FE analysis. This method was developed to enable high-speed analysis on CPU based large-scale computing environments. This method was ported to GPU environments using performance portable OpenACC for further enhancing the computational speed. Since the solver is designed to have high parallel performance and reduced memory access, it is possible to demonstrate high performance even on GPU environments. In OpenACC programming, it is possible to port computation to GPU by adding an OpenACC directive clause without making significant changes to the original program. In the operation of the matrix vector product based on the element-by-element (EBE) method [13], it is impossible to allocate temporary arrays for matrix-vector multiplication results considering device memory capacity. Thus, atomic operation is used to add the local matrix-vector multiplication results to a global result vector in EBE kernels on GPUs. Additionally, it has been shown that the performance is sufficiently high even when compared with a program well-tuned by CUDA [14]. For more detailed GPU implementation, see Ref. [15].

Inverse Analysis of Slip and Slip Deficit on the Plate Boundary
Fault slip is estimated using geodetic data on the ground surface. The slip distribution s(x) is expressed as: where X i represents the distribution function of the fault slip in the unit fault i, and a i is the model parameter. The Green's function, which represents the displacement of the observation point corresponding to each unit fault, is calculated by FE analysis, and the observation model is expressed as: Here, d represents the observation point data, and H represents the Green's function. It is assumed that error e follows the normal distribution of the variance-covariance matrix V with a mean value of 0. In addition to the observation errors, e includes errors, such as the uncertainty of the crustal structure model and the high-frequency component of the slip distribution that cannot be expressed by the unit fault. However, we assume that the error corresponding to the observation accuracy independent for each observation point is included in e. Due to the smaller number of observation data compared to the number of unit faults, and errors in the observation data and crust model, inversion may become unstable. To stabilize the inversion problem and obtain a reasonable solution, it is necessary to perform regularization using a priori information on the slip distribution. Here, the model parameters are determined by minimizing the objective function, where a T Ga is a term used to constrain the smoothness of the slip distribution [16].
Since it is impossible to assume the area of seismic slip in advance, unit faults are set wider than the area in which the seismic slip actually occurs. As unit faults are defined as a local fault slip, and the actual slip distribution is estimated as a superposition of unit faults, a sparse solution is expected. Therefore, by adding a Group Lasso term to Equation (6), the objective function considering sparsity is expressed as: Here, g represents a group of variables that correspond to unit faults with different strikes in the same region. This objective function is minimized by alternating direction method of multipliers [17]. As shown later in the results, negative slips appear in some cases for the coseismic slip estimation. Therefore, we also estimated the coseismic slip as the purely opposite direction of the plate convergence by the objective function with non-negative constraints added: Here, s represents variables for the opposite direction of the plate convergence. For all the above objective functions, hyper parameters of the weight among the terms ν and ξ are determined by k-fold cross-validation (CV) [18]. Here, the displacement data in all directions measured at one observation point are regarded as one observation data point. In the iteration in CV, the sum of the weighted squared residuals of the validation set and the predicted value based on the slip distribution estimated by the validation set are used as the CV score.

Numerical Experiment
Herein, we investigate the estimation performance of the coseismic slip and interseismic slip deficit off the coast of Nankai using the proposed method. Since the Nankai Trough earthquake has not yet occurred since the advanced observation network was installed, we used artificial observation data generated by the virtual coseismic slip based on a large event following Hyodo and Hori [19] (Figure 1). Specifically, a 3D FE model was generated based on the proposed Japan integrated velocity structure model version 1 (JIVSM) [4,5], and artificial observation data were generated by inputting a virtual Nankai Trough coseismic slip and interseismic slip deficit. Here, we used the displacement data at the observation station as the observation data. After adding plausible observation errors to the artificial observation data, we confirmed whether a coseismic slip and interseismic slip deficit could be estimated by inverse analysis [20]. Next, we confirmed the effect of the proposed ambiguity of the physical properties of the crustal structure on the estimation performance. First, we explain the 3D FE model. A numerical simulation model for computing Green's function is constructed based on the JIVSM seismic velocity structure model. The coordinate transformation method of the crustal structure model and the simplification of the layer structure is conducted following Hori et al. [6]. The crustal structure data are projected onto the Cartesian coordinate system with its origin at east longitude 135 • , latitude 33.5 • , and the ellipsoidal height 0 m. The target area domain is set between −1248 km < x < 1248 km, −1248 km < y < 1248 km, and −1100 km < z. Figure 2 shows the FE model generated by discretizing the target crustal structure model with the minimum element size in the horizontal direction ds = 2000 m. From the examination result of Hori et al. [6], an element size of around 500-m is required for numerical convergence of the results; thus, elements are refined such that the element size in the horizontal direction is 500 m in the range of −448 km < x < 448 km, −448 km < y < 448 km, and (h(x, y) − 10) km < z < (h(x, y) + 10) km. Here, h(x, y) is the z-coordinate of the plate boundary (for region beyond the trench axis, h(x, y) represent the z-coordinate of the ground surface). The generated FE model comprises 5.7 × 10 8 tetrahedral quadratic elements with 2.3 × 10 9 degrees of freedom. The unit fault defined by [6] is used as the basis function. Four hundred and sixty subfaults are set on the plate interface, and two types of unit faults, one in the plate convergence direction and another in the orthogonal direction component, are set respectively. Therefore, it is necessary to compute 460 × 2 = 920 Green's functions per crustal structure model. In order to analyze 104 models, including the model of the original physical property values, 104 × 920 = 95,680 Green's functions were computed. Here, a seven compute-node cluster each with a 2.80 GHz 24-core AMD EPYC 7402P CPU and four NVIDIA Tesla A100 GPUs were used in computing the 95,680 Green's functions (i.e., 28 GPUs were used in parallel with 28 MPI processes). Generally it would require a very large computation cost to compute 95,680 Green's functions using the FE model with 2.3 × 10 9 degrees of freedom. However, the analysis time is greatly reduced by using a suitable analysis algorithm and GPU implementation, which allowed for computing all Green's functions in 1.9 × 10 6 s (22.0 days).
Based on the larger event in the study reported by Hyodo and Hori [19], the virtual Nankai Trough coseismic slip and interseismic slip deficit are set. However, it is moved so that the slip reaches the trench axis, and the slip amount is set to 0 at regions where the unit fault is undefined. The coseismic slip is set to the opposite direction of the plate's subduction. The reference interseismic slip deficit distribution is set in the opposite direction of the coseismic slip based on Hyodo and Hori [19], with its maximum value scaled to 6.5 cm (Figure 1). is based on the study by Yasuda et al. [21] and Tadokoro et al. [22], DONET real-time is based on the study by Wallace et al. [23], DONET post-processed is based on the study by Machida et al. [24], and others are based on the study by Murakami et al. [20]. The ocean bottom pressure data of N-net under construction from off Shikoku to off Hyuga-nada is also added for the observation data during earthquakes, with the observation error levels set based on the similar S-net system [20]. Observation data with 1000 patterns of independent observation noise is generated for considering the wide variety of observation errors.  To evaluate the influence of the ambiguity of the crustal structure data to the estimation accuracy, we generated models in which the physical property values of the five layers presented in Table 2 (i.e., sedimentary layer 1, sedimentary layer 2, upper crust, lower crust, and continental mantle) were changed, and the coseismic slip and interseismic slip deficit were estimated based on the computed Green's function for each model. The range of p wave velocity v p is set based on JIVSM and the survey results of Nakanishi et al. [7] ( Table 2). As the s wave velocity v s and density ρ follows a linear trend with respect to v p when similar layers are grouped in JIVSM, we set v s and ρ from the value of v p based on the linear regression by the least squares method. Furthermore, the sedimentary layers were grouped into two layers in this model, and a grid search was performed for 100 cases. Regarding the upper crust, lower crust, and continental mantle, the lower limit of the range of assumed physical property values was set as the original physical property values. Three additional cases, where the values of the upper crust, lower crust, and continental mantle were changed to the upper limit value of the assumed property range, were computed.  For the ambiguity of the sedimentary layers, for sedimentary layers 1 and 2, 10 candidates are set so that they are equally spaced from the range, and 100 physical property values from the combination of these candidates are adopted. In addition, three cases are considered for the upper crust, lower crust, and mantle, each with an upper limit value. Here, the upper limit for the upper crust is larger than the value for the lower crust when given, but this was set to see the effect of the extreme case.

Layer ID in JIVSM v p (m/s) v s (m/s) ρ (kg/m 3 )
Sediment

Estimated Result without Ambiguity of Crustal Structure Model
Using the proposed 3D crustal structure model, we evaluated the performance of coseismic slip and interseismic slip deficit estimation. At first we estimated them via optimization expressed in Equation (6). As shown in Figure 4a, the overall behavior could be captured for the coseismic slip. However, the estimation accuracy deteriorates locally far away from the observation point and in places where a steep change is in the reference slip distribution (Figure 4b). The maximum standard deviation is 0.48 m, which is smaller than the maximum slip of 20.58 m and the maximum error of 5.8 m, indicating that the variation in estimation results due to the diversity of measurement errors is small. We can see that regularization suppresses the instability of inversion caused by observation noise. The standard deviation tends to increase in areas where the amount of slip is large, and no observation points exists. In the following, we use the relative error, as an index of the degree of agreement with the reference slip. While the mean value of relative error for the reference slip was 0.103, the standard deviation of relative error had a very small value of 0.006; thus, we can see that a certain estimation accuracy is obtained regardless of the diversity of the observation errors. Figure 5 shows the ground surface displacement distribution predicted from the mean of the estimated coseismic slip distribution. As horizontal displacement cannot be observed in real-time seafloor observation systems, large differences of 4.30 m and 1.51 m occur in the x and y components of surface displacement at observation points when comparing the cases for the predicted and reference coseismic slip. On the other hand, the maximum error for the z component was 0.12 m; thus, we can see that parameters consistent with observation values were obtained. From the above, we can see that it is difficult to accurately estimate the horizontal component of surface displacement from only the vertical component observations. Next, the interseismic slip deficit was estimated using the crustal structure model with the original physical property value. Although the observation accuracy is about the same level between real-time and post-processed data, the amount of interseismic slip deficit per year is 3.15 × 10 −3 times smaller than that of the coseismic slip. Therefore, the importance of the observed values in estimating the interseismic slip distribution is relatively smaller, and the geometric mean of the selected ν was 7.38 × 10 5 times larger than that of the coseismic case, indicating that the regularization of the spatial smoothing worked more strongly. The mean relative error of the estimated interseismic slip deficit is 0.178, which is significantly larger than the estimation error of 0.103 in the estimation of the coseismic slip distribution. The estimated mean of the interseismic slip-deficit distribution shown in Figure 6 has a smoother distribution, and the error is mainly in the low-frequency components. The standard deviation tends to be large in areas near the trench axis with a large slip.

Estimated Result with Ambiguity of Crustal Structure Model
The lowest relative error when changing the properties of the sedimentary layers was 0.097 for (v p1 , v p2 ) = (2400, 3940) m/s, while the largest error was 0.11 for (v p1 , v p2 ) = (1700, 5500) m/s. In the range of v p considered, the relative error decreased as v p1 for sedimentary layer 1 increases, while the relative error was minimized around v p2 = 3940 m/s for sedimentary layer 2 (Figure 7). A similar tendency also exists for the CV score when the optimum hyper parameters are selected, which deviates from the original physical property values (red dot in Figure 7). Therefore, it is impossible to predict the correct physical property values of the sedimentary layer from the estimation results. We can see that the effect of ambiguity of the physical properties of the sedimentary layer on the estimation results of coseismic slip is smaller than the level that can be constrained by the observation data. The relative error when different physical property values were applied to the upper crust increased to 0.127, indicating that the estimation accuracy was significantly reduced. In contrast, estimation accuracy did not change significantly when the physical property values of the lower crust and the continental mantle were changed. Figure 8 shows the plate convergence direction component of estimated slip distribution and the reference solution. In the case of changing the physical property values of the sedimentary layer, there is no remarkable tendency regarding the estimation distribution and estimation error even for the worst parameter case. On the other hand, when the properties of the upper crust are changed, slip in the opposite direction was estimated in a region where slip was zero in the reference solution. In any noise pattern, these negative slips systematically appear only when the physical properties of the upper crust are changed ( Figure 9). Additionally, a lower frequency component appears in the difference in the slip distribution in a wide region when the properties of the upper crust are changed (Figure 10). This is due to the inconsistency of the high-frequency components of the Green's function, which affects the inverse analysis results. In the case where the physical properties of the upper crust is changed, the geometric mean of ν, which is related to spatial smoothing, is 4.32 times larger than that of the model with correct physical properties; thus, we can see that strong regularization is applied. This trend occurs because it is difficult to reproduce the slip distribution when the noise different from the assumed observation noise (i.e., the error in Green's functions) is included in the CV. The mean CV score when using the adopted ν was significantly worse in the case where the physical value of the upper crust was changed to 4880 when compared to the score when using the correct value (2760), which indicates that the observed value is unpredictable due to the error in the Green's function (Table 3).   Next, we investigated the improvement in the estimation accuracy in estimating coseismic slip when the physical properties of the upper crust, which had a large estimation error, were changed. To suppress the occurrence of slip in the direction opposite to the correct slip direction, estimation was performed by adding a sparse constraint (Equation (7)). Although slip in a part of the seafloor is reduced due to the effect of the Group Lasso constraint, parameters making a sparse estimation was not selected, and the reverse slip remained unsuppressed ( Figure 11). We consider the reason is that erroneous estimation of the negative slip does not necessarily reduce the explanation of the observation data when the Green's functions used are incorrect. In order to confirm our idea, a non-negative constraint condition on the slips in the plate convergence direction component was applied to forcibly obtain a coseismic slip distribution without negative slip, and we examine whether the relative error increases or not. When the model with the original physical property values was used, the negative slip was hardly estimated without non-negative constraints. Therefore, the influence of non-negative constraints should be weak, and actually, the estimated coseismic slip distribution and relative error remain unchanged (Figure 12a). In contrast, when the physical property values of the upper crust were changed, the model parameter in which the negative slip had occurred was suppressed to zero, and the range of slip distribution was more accurately estimated ( Figure 12b). However, the relative error increased to 0.142 from the relative error of 0.127 in the case without the non-negative constraints. This confirms that the model error is the reason why the negative slip did not disappear even if the sparse constraint was introduced. From the above, we can assume that the physical properties of the crustal structure model may be set incorrectly if a negative slip is estimated in actual data analysis, as negative slip is not anticipated in an actual coseismic slip. Therefore, we can expect improvement in the estimation accuracy of physical properties by estimating the slip (without nonnegative constraints) using Green's functions of models with varying material properties within the permissible range, and selecting the physical property value that results in non-negative slip.   In the case of interseismic slip deficit distribution, no significant change was obtained in the relative error when changing the physical properties, even for the case when changing the properties of the sedimentary layer and the upper crust ( Table 3). The interseismic slip deficit distribution also shows similar results regardless of changes in the physical property values of any layer ( Figure 13). This is expected to be due to the lower S/N ratio in the interseismic observation and effect of regularization, which suppresses the effect of change in the high-frequency component of the Green's functions.

Conclusions
We investigated the effect of the ambiguity of the crustal structure on the estimation performance of the coseimic slip distribution and the interseismic slip-deficit distribution by calculating multiple Green's functions with different physical property values of the 3D crustal structure by utilizing the GPU-based 3D crustal deformation analysis, which can significantly reduce the analysis cost. Here, we changed the physical property values of the sediment, upper crust, lower crust, and mantle on the upper plate side, which are located above the plate boundary fault that causes coseismic slip and interseismic slip deficit. For the case of interseismic slip deficit estimations, no effects exceeding the observation error levels were found, while for the case of a coseismic slip, no significant difference in estimation performance was obtained except when the physical property values of the upper crust, with which most of the fault slip area is in contact, were changed. However, the estimation of coseismic slip when the physical property values of the upper crust were changed had a large relative error, and the estimation results were systematically biased. Specifically, slip distribution, including components opposite to the direction of the reference slip, was estimated in the deep part of the plate boundary where slip was not originally given. The negative slip is unsuppressed even when using sparse constraints in the estimation, and the relative error became even larger in the case when imposing non-negative constraints to prevent negative slip. From these results, it is suggested that incorrect physical properties are set in the upper crust if a negative slip is estimated using coseismic data, and that a more accurate physical property value of the upper crust can be estimated by finding a physical property value that is unlikely to cause negative slips. While only the observation data of the ground surface displacement is used in this study, the volumetric strain and the inclination change are observed in real-time at several points below the seafloor. Th evaluation of the improvement in the estimation accuracy of coseismic slip and interseimic slip deficits by adding other physical quantities indicating seafloor crustal deformation is our future work.