Predictive Simulation for Surface Fault Occurrence Using High-Performance Computing

Numerical simulations based on continuum mechanics are promising methods for the estimation of surface fault displacements. We developed a parallel finite element method program to perform such simulations and applied the program to reproduce the 2016 Kumamoto earthquake, where surface rupture was observed. We constructed an analysis model of the 5 × 5 × 1 km domain, including primary and secondary faults, and inputted the slip distribution of the primary fault, which was obtained through inversion analysis and the elastic theory of dislocation. The simulated slips on the surface were in good agreement with the observations. We then conducted a predictive simulation by inputting the slip distributions of the primary fault, which were determined using a strong ground motion prediction method for an earthquake with a specified source fault. In this simulation, no surface slip was induced in the sub-faults. A large surface slip area must be established near a sub-fault to induce the occurrence of a slip on the surface.


Introduction
Since the occurrence of massive earthquakes in Taiwan and Turkey in 1999, there have been growing concerns regarding potential damage to various infrastructure systems and buildings caused by surface fault ruptures. For on-site fault assessment at nuclear power plants (NPPs), it is important to estimate fault displacements and their impact on facility safety functions. Important NPP facilities are separated from primary (or main) faults, which are direct extensions of earthquake source faults, and detailed geological surveys are conducted. However, there are several cases of secondary faults (or sub-faults) located beneath or close to these important facilities. The effects of sub-fault slip on these facilities are being discussed.
Numerical simulations based on continuum mechanics are promising evaluation methods for the estimation of surface fault displacements. However, there are major difficulties in simulating fault rupture processes. One challenge is the requirement of significant numerical computations to simulate such processes for a target area that is only a few hundred meters in size. In practice, the following two-step simulation method is adopted [1]: (1) evaluating the boundary displacement of an area containing a target facility by determining the crustal deformation based on the elastic theory of dislocation [2] and (2) evaluating the deformation of the target area using a detailed two-dimensional model with high resolution and fidelity. A few researchers have conducted dynamic fault rupture simulations to evaluate ground motions or rupture propagation using the finite difference method, boundary element method, or finite element method (FEM) [3][4][5][6][7]. It was difficult to conduct dynamic rupture simulations at a sufficiently high resolution because it requires a large amount of computation resources. In recent years, however, fault displacement analyses using dynamic rupture simulation with high-resolution detailed analytical models have begun to be reported [8][9][10][11].
Another difficulty is the stability loss in the solution of the initial and boundary value problems to which numerical analysis is applied. Stability implies that a solution does not change significantly when a small variation is added to the problem and stability loss leads to drastic changes in solutions induced by small disturbances in numerical computations.
We overcame these difficulties by applying high-performance computing (HPC) to solve a high-resolution analysis model. We developed an FEM strengthened with the following two functions: (1) a symplectic time integration explicit scheme to conserve the energy of the system and (2) rigorously formulated high-order joint elements. The FEM was also enhanced with parallel computing capabilities [12,13]. Another use of HPC is to perform multi-case analyses with different input conditions in parallel. This enables us to evaluate the variability of the solutions.
We applied the FEM to the simulation of the 2014 Nagano-ken-hokubu earthquake, where surface faults were observed [14]. The primary fault of this earthquake slipped in the form of a reverse fault, including a left-lateral strike-slip component. Through simulations, we confirmed the applicability of the FEM to the evaluation of surface fault occurrence. We also estimated the effect of material uncertainty on surface fault displacement by performing 230 simulations with Latin hypercube sampling [15].
In the paper, we applied the FEM to simulate the 2016 Kumamoto earthquake, where the primary fault slipped in the form of a right-lateral strike-slip fault, including a normal slip component. The applicability of the FEM was further confirmed for this earthquake because the observed fault displacement was reproduced in the simulation. We also examined the potential for the predictive simulation of surface fault occurrence using the FEM. The same analysis model used for the 2016 Kumamoto earthquake was adopted, but the input earthquake fault at the bottom of the analysis model was changed. The conditions under which secondary faults induce surface fault occurrence are discussed based on our predictive simulations.

Estimation of Fault Displacement
For fault displacement simulations, we constructed a continuum model of the ground and faults (Figure 1). An input slip ∆ was defined at the bottom of the primary (or main) fault and the surface slip δ was obtained after calculating the spread and dispersion of the slip on the fault plane. For on-site fault assessment at NPPs, the surface slip on branch or secondary faults (or sub-faults) must be estimated, rather than the slip on the main fault. A fundamental difficulty in simulating the fault rupture process is the stability loss in the solutions of the initial and boundary value problems, which leads to significant changes in the solution as a result of small disturbances.
The treatment of uncertainties stemming from limitations in the quality and quantity of available relevant data for underground structures, stress states, and source fault dynamics is important. Numerous simulations must be conducted under different conditions (capacity computing).
When significant uncertainty exists, the surface fault displacement δ calculated using capacity computing is distributed over a wide range. In this case, δ was unpredictable. It is important to evaluate the critical input slip ∆ c , which is the minimum input slip that causes a surface slip, before evaluating δ. This concept is illustrated in Figure 2. Quasi-static simulations facilitate the determination of the relationship between the input and surface slips. This is extremely important for the estimation of ∆ c . Sawada et al. [13] demonstrated that the calculated ∆ c and δ from quasi-static simulations were close to those calculated using dynamic simulations with damping and proposed the use of quasi-static simulations for capacity computing. Therefore, quasi-static simulations were adopted in this study. Therefore, the symplectic time integration described above was not used and seismic waves were not included in the simulations.

FEM for the Fault Displacement Problem
We implemented a rigorous joint element for modeling the slip behavior of faults and a symplectic time integration method with a Hamiltonian formulation using the opensource FEM software FrontISTR [16] to develop a numerical tool for fault displacement simulation [17].
In joint elements that are widely used in the field of rock mechanics, nodal forces are directly calculated from relative nodal displacements on a discontinuity [18]. By using this method, element stiffness matrices can be easily calculated. However, they differ from the standard method for constructing finite elements. Joint elements for fault displacement simulations should be defined to ensure proper convergence for the chosen element sizes and applicability to curved faults. We developed a rigorous joint element in which element stiffness matrices are derived from the Lagrangian on a discontinuous surface using an isoparametric formulation [12]. The following is a formulation for a second-order triangular joint element.
The variant of the Lagrangian δL J is given by the following area integral on the discontinuous surface S: where u(= u + − u − ) is the relative displacement; u + and u − are the displacements on the top and bottom surfaces of S, respectively; n is the unit normal vector of S; t is the unit tangential vector of S; and k n and k t are the normal and tangential spring constants, respectively. A ⊗ B denotes the tensor product of vectors A and B.
Considering the case that S is a curved surface, an isoparametric formulation is adopted. The three-dimensional global coordinate system x of S is given by the following function of the two-dimensional local coordinate system ξ (see Figure 3): where x a (a = 1, 2, . . . , 6) is the nodal coordinate of the joint element and N a is the shape function. The area integral on S is transformed into the area integral using ξ, and the Gaussian integral is applied: Figure 3. Three-dimensional second-order triangular joint element and its local coordinate system.
The unit normal vector n is formulated as follows: where A × B denotes the vector product of vectors A and B. The differential global coordinates s is defined as follows: The unit vectors have the following relationship: where e is the identity matrix. The element stiffness matrix E ab J is obtained as follows: We implemented joint elements using a nonlinear spring-type constitutive equation to represent fault movement [12,13]. The constitutive equation for the shear movement of the joint element is defined as follows (see Figure 4): where τ is the shear stress and u is the slip on the fault plane. The spring coefficient per unit area (shear stiffness) κ is described by the following function for a slip u: where κ 0 is the initial shear stiffness, κ d is the final shear stiffness, and u cr is the critical slip. We considered the effects of the normal stress on the shear stiffness of the faults and defined the following equation: where κ 0 is the shear stiffness per unit of normal stress, σ n is the normal stress on the fault plane, and b is a constant. We assumed that the first peak of the slip-shear stress relationship shown in Figure 4 is identical to the shear strength. The parameters κ 0 and b in Equation (9) can be determined such that the shear strength satisfies Coulomb's law of friction as follows: where c is the cohesion and φ is the friction angle. The earthquake occurred along the Futagawa and northern Hinagu fault zones [19]. Surface faults of 28 and 6 km in length were observed along the Futagawa and Hinagu fault zones, respectively (see Figure 5). Vertically, the north side of the Futagawa fault zone sank ≥1.0 m and the south side rose ≥0.3 m. Horizontally, the north side was displaced by more than 1.0 m to the east and the south side was displaced by more than 0.5 m to the west. The surface faults were investigated and the fault displacement distribution along the faults was obtained [20]. A right-lateral strike-slip fault of up to 2.2 m was observed at Dozon in the eastern part of Mashiki. The surface fault branched into two faults from Dozon to the west. Although these two branched-out faults exhibited a right-lateral strike-slip fault, a left-lateral strike-slip fault was observed between them. The target area of our numerical simulations was a 5 × 5 km square area containing the faults around Dozon, as shown in

Analytical Model
The target of our simulations was a 5 × 5 km square area containing the branched-out faults in Mashiki town. The target depth from the ground surface was approximately 1 km. The National Research Institute for Earth Science and Disaster Resilience provides a database of the crustal structure density and elastic wave velocity covering all of Japan on the website of the Japan Seismic Hazard Information Station (J-SHIS) [22]. In this database, elevation data regarding the ground surface and boundaries of elastic wave velocity and density are provided at the center of every grid of 30 arcseconds of latitude and 45 arcseconds of longitude. By using this database, we constructed an analytical model that includes the shape of the land and crustal structure. We determined the strike and dip angles of the Futagawa and Hinagu faults by referring to the paper by Asano and Iwata [21]. The dip angles of the Futagawa and Hinagu faults were 65 • N and 72 • N, respectively. Although the upper edges of the faults did not reach the ground surface, we extended them to the ground surface with the same dip angle. Figure 7 presents the surface faults along the Kiyama river in Mashiki town, where the Futagawa fault branches out into 2 faults. We treated the line labeled "Main", which is one of the two branched-out faults, as the main fault, which is identical to the Futagawa fault. The other faults, which are labeled as "SubN", were included in the analytical model as sub-faults. Three faults, namely Sub2, Sub3, and Sub5, were observed during the 2016 earthquake, with Sub2 being the other branched-out fault. Sub3 was a left-lateral strike-slip fault bridging the main fault and Sub2. Sub1 and Sub4 were faults known to exist before the 2016 earthquake. However, surface faulting was not observed along them during the earthquake.   [23]. Sub1, Sub2, and Sub4 are vertical. The ground is divided into three layers based on the J-SHIS database. The ground is discretized using second-order tetrahedral elements and the fault planes are discretized using second-order triangular joint elements. The fault plane element size is approximately 50 m and the total number of degrees of freedom is approximately 3.76 million. The white dots from a bird's-eye view in Figure 8 indicate the locations of the output points for the fault displacement results. "Y2780" indicates that the y coordinate of these plots is 2780 m. We assumed that the ground was a linear elastic body and derived the Young's modulus and Poisson's ratio of each layer from the density and elastic wave velocity values in the J-SHIS database. We conducted an initial stress analysis to determine the values of the parameters κ 0 and b on the fault planes using Equations (10) and (11) before fault displacement analysis. Normal displacements were fixed at the bottom and side boundaries in the initial stress analysis. We assumed a cohesion of 0.025 MPa and friction angle of 25 • on the fault planes by reference to the experimental value of the crushed zone in NPP sites [14]. The ratio between the initial and final shear stiffness κ d /κ 0 was assumed to be 0.01 in reference to the analyses conducted using a simple ground/fault model and a real earthquake ground model [13,14]. The critical slip u cr was 0.1 m in reference to the value used in the dynamic rupture simulation of the fault [4]. The material parameter values are summarized in Table 1. The size of the analytical model is presented in Figure 8 and is smaller than the length of the main fault. The boundary conditions of the analytical model were determined using large-scale analysis. We used the elastic theory of dislocations [2] with a slip distribution on the main faults, as shown in Figure 6, as the input. The upper edge of the fault is at a depth of 2 km. We extended the faults to the ground surface and assumed that the slip on the extension was the same as that just below the extension. The ground displacement was calculated and the forced displacement at the bottom of the FEM model was determined. The maximum slip on the bottom of the main fault was 5.2 m, which is twice the estimated 2.6 m slip. The other outer model boundaries were traction free. We refer to this simulation as case1. Figure 9 presents the displacement magnitude contour and deformation of the analytical domain. The deformation is presented with 300-times amplification. When ∆ is 2.6 m, which is the slip estimated by Asano and Iwata [21], there is a displacement gap along the main fault, Sub2, and Sub5. Additionally, there is a displacement gap along Sub1 and Sub4 when ∆ is 5.2 m, which is twice the estimated slip value. This gap indicates a surface slip on each fault.   Figure 11 presents the profiles of the surface slip for the strike-slip and vertical components along Main, Sub2, Sub3, and Sub5. Figure 12 presents the observed slip distributions along the Futagawa fault and Hinagu fault [20]. Comparing the calculated and observed strike-slip, the calculated slip is smaller on Sub2. On the other hand, the surface slip on Sub5 distributes in a longer section in the calculation. The sum of the calculated surface strike-slip values (approximately 2.0 m) is in good agreement with the total displacement of the southwestern section of the Futagawa segment in Figure 11. The vertical component of the calculated surface slip on the main fault is intermediate to the observed ones, but the sum of the vertical slips on the main fault and Sub5 mostly exceeds the observed ones.  Table 2. ∆ c is small for Sub2 and Sub5, where the surface displacement was actually observed in the 2016 earthquake. The sense of the slip at each location is presented in Table 2. The sense of the slip in the simulations was remarkably similar to that of the observed surface faults, as shown in Figure 12.

Discussion
We constructed a 5 × 5 × 1 km analytical model and derived the boundary conditions from the elastic theory of dislocation using a slip distribution in the main fault obtained from seismic source inversion analysis. When forced slip occurred at the bottom of the main fault, we observed a steep increase in the surface slip on the main and sub-faults. On the sub-faults that appeared in the 2016 earthquake, the surface slip appeared at a smaller input slip than on the other sub-faults. The sense of slip on each fault in our simulations was remarkably similar to the observations during the real earthquake.
As shown in Figure 12, the surface fault slip is distributed in an uneven manner. Although there are locations where significant slip occurs, there are many locations where the slip is small or does not occur in the entire section. The calculated surface slip distribution is continuous and agrees well with the measured total slip distribution (enveloped slip distribution curve). This indicates that the proposed numerical method enables a conservative estimation of surface fault slip.
The arrows in Figure 10 represent the Sub1, Sub2, and Sub5 contact edges with the main fault. The slipped sub-fault areas propagate from these contact edges. In the 2014 Nagano-ken-hokubu earthquake simulation [14], the slip on the sub-fault was larger near the surface than near the contact edge with the main fault. The slip distribution on the sub-faults was remarkably different between the two earthquake simulations.
However, some phenomena were not reproduced in our simulations. For example, the surface slip on Sub3 requires a large input slip. In the simulations, the small faults tended to be less likely to slip.
In general, faults are identified from ruptures left on the ground surface. It is difficult to include faults that cannot be identified from the ground surface in the model. Moreover, it is difficult to fully capture the subsurface structure and connectivity of faults. In the analysis of this paper, faults were extended to the model boundaries or other faults. The larger the fault plane, the greater the possibility that large displacements will be calculated. On the other hand, if many faults are identified, it is difficult to represent all of them with joint elements. As sub-faults, it is reasonable to include in the model faults that are relatively large or that affect the focused facilities. In general, if the number of faults is limited rather than considering many faults, the fault displacement for each single fault will be calculated as larger. However, this needs to be further investigated. For predictive simulation, we were not able to use the results from the seismic source inversion analyses shown in Figure 6 and we had to assume the slip distribution on the main fault. In this paper, we propose a method to set the slip distribution on the main fault and apply it to predictive simulation using an analysis model of the 2016 Kumamoto earthquake. We ignore the Hinagu fault and consider a scenario in which only the Futagawa fault ruptures in the following predictive simulations.

Predictive Simulation
We set the slip distribution using a strong motion prediction method called recipe [24], which includes the following steps: (1) set the main fault geometry (length, width, depth, and dip angle), (2) calculate the seismic moment, (3) calculate the average slip from the seismic moment and shear modulus of the crust, and (4) set the asperity number and location to calculate the slip on them and on the background. To calculate the seismic moment M 0 , the following equation is used for relatively large earthquakes [25]: where S denotes the square of the main fault. We adopted the geometry of the Futagawa Fault presented by Asano and Iwata [21]. Because S = 5.04 × 10 8 m 2 was derived from a length of 28 km and width of 18 km, M 0 was calculated as 1.41 × 10 19 Nm. The average slip D was calculated from S and the shear modulus µ using the following equation: when we used µ = 30 GPa, which was based on the assumption of a hard rock mass in the deep underground, D = 0.934 m was obtained. We set a uniform slip distribution of 0.934 m on the entire Futagawa fault, as shown in Figure 14a, and conducted our first predictive simulation. We refer to this simulation as case2. In this recipe, asperities were set on the seismic source fault. Subsequently, we set the slip distribution by considering the existence of asperities. Although it is recommended to use an empirical equation for the relationship between the short-period level and equivalent radius of the total asperity area, we determined the total asperity area to be 1.28 × 10 8 m 2 by referring to research on the ratio of the total asperity area to the total fault area introduced in the recipe. We placed 2 asperities of equal area (64 km 2 ) on the Futagawa fault, as shown in Figure 14b. The average slip of asperity D a was determined to be ξ times the average slip D as: where ξ = 2 was used based on an analysis of recent crustal earthquakes. The average slip of the background D b was derived from the following equations: where M 0a is the total seismic moment of the asperities, S a is the total area of the asperities, M 0b is the seismic moment of the background, and S b is the area of the background. We obtained D a = 1.869 m and D b = 0.616 m. We set these slip values on the Futagawa fault, as shown in Figure 14b, and conducted a second predictive simulation. We refer to this simulation as case3.

Large Surface Slip Area
When we set the slip distribution on the main fault based on the recipe, the surface slip was uniform and smaller than that observed during the 2016 Kumamoto earthquake. The slip along the surface fault was not generally uniform. Therefore, we set a large slip area on the main fault based on case2, as shown in Figure 15. The following empirically obtained equation proposed by Matsuda et al. [26] was used to set the maximum surface slip: where D max is the maximum surface slip and L is the surface rupture length. Here, we used the 28 km fault length for L and obtained D max = 2.8 m. Other equations have been proposed to relate surface fault displacement to the length of faults [27]. We can also use these equations. We assumed that the length of the large surface slip area was 5 km, which is the same as the size of the FEM model, and that the peak slip position was identical to the center of the FEM model. We refer to this simulation as case4. Additionally, we set up 2 cases in which the large surface area was shifted by +1 km and −1 km in the y direction compared to case4, which are referred to as case5 and case6, respectively.
In simulations case2 to case6, the maximum forced slip at the bottom of the FEM model was twice that of the slip distribution.   In Sub2, Sub3, and Sub5, where surface fault displacement can be observed, ∆ c is less than 2.8 m. This tendency is similar to the case1 results. The surface slip is slightly smaller in case4 than in case1. Although the slip distributions at depth are significantly different between case1 and case2, the surface slips on the sub-faults are similar.  Figure 18 presents the variation in the surface net slip δ on Sub2, Sub3, and Sub5 in case4, case5, and case6 as the input slip ∆ increases. The results for Sub5 are similar between the three cases. However, the ∆ c values on Sub2 and Sub3 are much smaller in case5 than in case4 and case6. As shown in Figure 19, the peak of the large slip area is close to the intersection line between the main fault and Sub2 fault. In contrast, in case6, no surface slip occurs on sub-faults Sub2 and Sub3 because the location of the maximum slip on the main fault is far from the intersection between the main fault and Sub2 fault. These results indicate that the slips on the sub-faults are strongly affected by the slip distribution on the main fault close to the sub-faults.

Discussion
We propose a method for setting the slip distribution on the main fault for predictive simulations based on the method of constructing characterized source models for the prediction of strong motion. In cases without asperities and with two asperities, surface slip on the sub-faults did not appear. These results suggest that surface displacement is less likely to occur on sub-faults. In the case with two asperities, a large slip was inputted into the asperities, and it had little effect on the surface slip on the sub-faults.
We also set up cases in which a large surface slip area was located near the sub-faults. We adopted an empirical equation for the relationship between the maximum surface displacement and fault length. We then placed the large slip area in a section of the FEM model. From this simulation, we obtained results similar to those of the simulation with the inversely analyzed input slip. The location of this large slip area had a strong effect on the surface slip of the sub-faults. In other words, for surface slip to occur on a sub-fault, a large slip must occur on a nearby main fault. The probability of a large slip being very close to a sub-fault is relatively small and this condition is extreme.

Conclusions
We developed a parallel FEM program to estimate surface fault displacement. In this study, we applied a numerical method to simulate the 2016 Kumamoto earthquake, where the seismic source fault ruptured in the form of a right-lateral strike slip fault, including a normal fault component. We modeled a 5 × 5 × 1 km domain, including branched-out faults and other sub-faults in Mashiki town. We applied forced displacements on the bottom surface of the model based on the slip distribution on the main faults and elastic theory of dislocation. As the input slip increased, surface slips appeared on both the main and sub-faults. In the simulations, the surface slip appeared at a smaller input slip on the sub-faults compared to the observed values during the 2016 earthquake. The sense of slip on each fault in the simulations was in good agreement with that observed during the earthquake. Furthermore, the calculated surface slips agreed well with the envelope curve of the measured surface slip distribution. This was the second application of the numerical method to real earthquakes. Because our simulations reproduced the main features of two different earthquakes with reasonable accuracy, the proposed numerical method should be applicable to surface fault displacement estimation.
We also proposed a method for setting a slip distribution on main faults for predictive simulation based on a method for strong motion prediction. However, surface slip did not appear on the sub-faults, even in the case with two asperities. We then set up a large slip area near the sub-faults based on an empirical equation of the relationship between the maximum surface slip and the length of a fault. In this simulation, we obtained results similar to those of the simulation with the inversely analyzed input slip. We consider that the slip distribution on the main fault with a large slip area is an extreme case in which surface slip on the sub-faults is likely to occur.
In this paper, we presented fault displacement evaluation using quasi-static simulation. As a next step, we aim to evaluate the superposition of fault displacement and earthquake ground vibration by dynamic analysis that takes into account the rupture propagation of the fault.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.