A Feasibility Study of an FEM Simulation Used in Co-Seismic Deformations: A Case Study of a Dip-Slip Fault

For this study, we conducted a numerical simulation on co-seismic displacement for a dip-slip fault in a half-space medium based upon a finite element method (FEM). After investigating technical problems of modeling, source and boundary treatment, we calculated co-seismic deformation with consideration to topography. To verify the numerical simulation results, the simulated co-seismic displacement was compared with that calculated using a dislocation theory. As a case study, considering the seismic parameters of the 2008 Wenchuan earthquake (M 8.0) as a source model, we calculate the co-seismic displacements with or without consideration of the terrain model in the finite element model to observe terrain effects on coseismic deformation. Results show that topography has a non-negligible effect on co-seismic displacement, reaching from -11.59 to 4.0 cm in horizontal displacement, and from -3.28 to 3.28 cm in vertical displacement. The relative effects are 9.05 and 2.95% for horizontal and vertical displacement, respectively. Such a terrain effect is sufficiently large and can be detected by modern geodetic measurements such as GPS. Therefore, we conclude that the topography should be considered in applying dislocation theory to calculate co-seismic deformations.


INTRODUCTION
Science dislocation theory was first introduced in seismology by Steketee (1958). Many scientists developed analytical expressions to calculate co-seismic deformations for different earth models and fault types. These theories are applied widely in interpreting geodetic observations. Okada (1985Okada ( , 1992 presented a set of mathematical formulas based upon previous studies. His expressions are used widely in calculating co-seismic deformations for a half-space medium model. The Earth is not a perfectly elastic medium, but instead is a viscoelastic body. When an earthquake occurs, the earth not only undergoes elastic deformation but also experiences postseismic deformation, which refers to viscous relaxation at long time scales. To study co-seismic and post-seismic deformations, Wang et al. (2003) and Wang (2005aWang ( , b, 2006) developed a dislocation theory for a lay-ered viscoelastic half-space earth model. His computing code is useful in calculating co-seismic and post-seismic deformations of an elastic Maxwell, or Kelvin medium model (Zhang et al. 2008).
Furthermore, to consider the effect of curvature and layer structure, Sun et al. (1992Sun et al. ( , 1996Sun et al. ( , 2003Sun et al. ( , 2004 and Sun and Okubo (1993, 1998, 2002 proposed a dislocation theory for a spherically symmetric, non-rotating, perfectly elastic and isotropic (SNREI) Earth model to compute co-seismic deformations, such as potential, gravity, geoid, and strain changes. Recently, Fu (2007) and Fu and Sun (2007) investigated co-seismic deformations for a three-dimensional inhomogeneous earth model.
In sum, dislocation theories have been developed for different earth models. However, these dislocation theories sometimes cannot explain geodetic observations well. The reasons might be regarded as two: the error of geodetic observation or, the earth model used in dislocation theories because all the dislocation theories above are based on a regular geometric earth model, without consideration of the surface topography. In fact, it is difficult (even impossible) to model topography in any dislocation theory presented above (Okada 1985(Okada , 1992Okubo 1993, 1998;Fu and Sun 2007). To date, the topography effect on co-seismic deformation remains to be an open question (Wang et al. 2003;Wang 2006). Therefore, in this study, we strive to investigate topographic effects on computing co-seismic deformations, using the FEM method. Finite element numerical simulation is an effective method because it can accommodate the Earth medium with any geometrical shape because of the development of computers and computing technology. It can easily build a real earth model to facilitate calculations to include the effects of topography. In dealing with a fault in the FEM method, Melosh and Raefsky (1981) presented a split node method. In this method, the value of the displacement at a single node point shared by two elements depends upon which element it is chosen, thus introducing a displacement discontinuity between the two elements. Melosh and Raefsky (1981) show that the degrees of freedom do not increase by splitting. Later, Zeng and Song (1999) and Zhu and Wang (2005) presented a double-node technology in dealing with the fault in the FEM method. The double-node method is formally different from the split node method since one node is doubly numbered at the fault plan. However, it is essentially the same as the split node method, because both methods give the same description of the relative movement at the fault plane. However, the double node method is convenient in practical FEM computation and result analysis in comparison with the split node method. In this study, we adopt the double-node technology, to investigate the terrain effect on co-seismic deformations and apply it in the 2008 Wenchuan earthquake (M 8.0).

FINITE ELEMENT TECHNIQUE OF THE FAULT MODEL
When the finite element method is applied in computing co-seismic deformations, the important steps are the division of the finite element and the focal fault handling. In this section, we discuss them with different modeling and numerical tests. To verify the correctness and accuracy of the finite element numerical simulation results, we compare numerical results with that of Okada's analytical solution. For simplicity, the simulation is performed for a two-dimensional finite element, and three different inclinations of the fault: horizontal fault, 33° dip-slip fault and 90° vertical fault.

Fault Geometry Model and Element Division
We first consider a fault buried in a geometric model region with a range of 1000 km × 500 km (width × height).

Boundary Conditions and Fault Displacement Load
The lower edge of the computational domain is a fixed boundary. On the upper boundary is the free boundary. The free boundary condition is set around the X direction of nodes along both sides of the border of the model, but the fixed boundary condition is set around the Y direction of nodes, which is the zero displacement constraint. For any computational domain (in theory) there must exist a boundary effect: because of the selection of the border at the lower boundary and both ends, it must have a certain error. In general, the boundary effect decreases with distance from the border. It can usually be ignored for a fault scale that is greater than 10 to 30 times larger. In the regional model of Fig. 1, the fault to transverse area ratio is 1:50. Its boundary effects are expected to be minimal. The detail will be discussed in the section 5.
Addressing the fault is the most critical part of using finite element method to study the dislocation problem. As described above, we considered a fault plane in the middle of a region with 20 km to the surface. To construct a mutual dislocation model, we make the fault plane process into an internal boundary plane that comprises a number of doublenodes. They are given in the opposite direction of the displacement constraints, as shown in Fig. 1b. The left side of the fault plane is numbered as 1, 2, 3, ..., 13; the right side is correspondingly numbered as 26, 25, 24, ..., 14. The first and the 26 th stand for the left and right of one point, which are respectively shown on the left plate boundary surface of the plane and the right plane boundary surface. When dislocation is imposed, the left side of the nodes has an upward load displacement, whereas a downward load displacement occurs on the right side of the nodes. Thereby, the study area of the dislocation model is established.

Selection of The Model Medium Parameters
Current research is limited to the elastic medium model, so the selected model medium parameter is only Poisson ratio υ. For simplicity and in contrast to the simulation results and the theoretical results presented by Okada (1985), all calculations in this study use the same Poisson ratio: υ is equal to 0.25.

ANALYSIS OF SIMULATION RESULTS
For the 90-degree inclination of the dip-slip fault model, we first give the displacement to the double-node fault. The X direction is fixed, with zero displacement constraints; the Y direction is applied to 2 m, but in the opposite direction of the displacement constraints. The surface displacement field of the finite element simulation is shown in Fig. 2. Comparison with the analytical formulas of Okada (1985) was completed. The results of these two methods shown in Fig. 2 are not consistent, and reveal a great mutual difference. From the figure, it is apparent that the surface displacement deformation in the horizontal is almost completely inconsistent. Near the fault the vertical displacement is consistent, but a large difference exists in the epicentral distance of 25 to 400 km range. In the finite element numerical simulation, the boundary condition of the displacement constraint is correct, but an error exists in the load on the fault at the nodal displacements.
Next we sought to make the finite element simulation results and Okada theoretical solution results consistent. After long-term exploration, we found a problem in dealing with the node displacement load of the fault. The X direction dislocation nodes should not be fixed, but the fault model should be treated as a relative sliding displacement. Therefore, the X direction is free, with no displacement constraint. The Y direction was given a relative displacement slide at the double-nodes in the fault plane. These new simulation results and Okada's analytical solution are shown in Fig. 3. It is apparent from Fig. 3 that these finite element simulation results and Okada's theoretical solution are completely consistent, which is readily apparent at the nodal displacement. The imposed conditions of the node are entirely correct.
The inclination of the 90° slip fault is a special case presented to verify the correctness of the finite element model and boundary condition imposed, constructing two new dislocation models as follows: a horizontal fault whose fault width is 20 km (Figs. 1c and d)  The different inclinations of dip-slip fault numerical simulation results show that it is necessary to impose the correct boundary conditions in the finite element simulation. The double-node processing of the fault plane and the relative displacement of the boundary conditions are correct.

CORRECT CONSTRAINTS OF THE BOUNDARY CONDITIONS
Although the results of the two methods presented in the previous section show good agreement and the absolute displacement is smaller, the boundary effect remains and the relative error between the numerical solution and theoretical method is very large. Careful comparison of the calculated region at both ends of the distribution difference shows that when the model is selected, the finite element numerical  solution and theoretical solution are still not completely consistent in the finite element model boundary region. To eliminate the boundary effect of the finite element model, we use the prior calculation of the theoretical value as the new constraints of the two boundaries, which will greatly reduce the boundary effect. Therefore, we use the values calculated using mathematical formula of the co-seismic deformation field of the Earth's interior as the initial values of the two boundaries of the finite element model. Results obtained using the two methods are completely consistent, the three dip-slip fault model calculations are shown in Figs. 6, 7, and 8.

SELECTING THE AREAL SIZE OF THE FINITE ELEMENT MODEL
Selection of the computational domain size is an important issue when calculating the dislocation deformation using the finite element method. This reduces the boundary effect and is an important prerequisite for ensuring the accuracy. Therefore, in this section, we discuss a means to ensure sufficient calculation accuracy in selecting the size of the computational domain.
To calculate the co-seismic displacement in this study, we examine a case with 10 km fault width, with a fault distance of 20 km to the top of the model, for example, a 2 m slip on the 90 degree fault, by increasing the computational domain with a change of length and width at the same time. Then we observe the convergence of the calculation results to determine the correct selection of the area size of the finite element model. Calculated results of the horizontal and vertical displacement are shown respectively in Figs. 9 and 10. The horizontal axis represents the ratio of the horizontal length of the model region and the fault width. The vertical axis represents the contrast between the calculated minimum or maximum value and the analytically obtained result. As Figs. 9 and 10 show, the finite element numerical solution has been fully consistent with the analytical solution after the ratio is 30 times between the length and the fault width in the finite element model. Therefore, to ensure the calculation accuracy when resolving the seismic dislocation deformation problem, the correct size of the area of the finite element model must be 30 times larger than the fault.

CASE STUDY: WENCHUN EARTHQUAKE
As described earlier, when we use a comprehensive dislocation theory to calculate the co-seismic displacement, its effectiveness is only acceptable for a regular geometric Earth model such as a spherical or infinite space, without considering topography. However, because of the presence of the terrain, its impact to theoretical calculation is unclear. For example, the surrounding geological environment of the 2008 Wenchuan earthquake showed that the Qinghai-Tibet Plateau to the Sichuan Basin has a large undulating terrain. The terrain profile is shown in Fig. 11. It is worthwhile to study the terrain using the classical dislocation theory to calculate the co-seismic deformation. Therefore, the impact of the terrain effect can be discussed as the following. Select the parameters of Wenchuan earthquake: 33° dip angle, 40 km fault width, and a 4.06 m slip, which is the outcrop to the surface (Ji and Hayes 2008). For comparison, we considered two models. One does not incorporate the terrain model; the other incorporates topography based on the model. After calculating the co-seismic displacement in the two models, the difference is the topographic effect. Results of the co-seismic difference displacements are shown in Fig. 12, with a horizontal difference displacement on the left and vertical difference displacement on the right. The figure shows that topography has a great impact on co-seismic displacement. In the horizontal direction, the displacement range is from -11.59 to 4.0 cm. The range of displacement in the vertical direction is small from -3.28 to 3.28 cm. Such a great surface deformation is detected sufficiently by modern space-based observations such as GPS. Figure 12 further indicates that the magnitude of terrain effect is basi-cally proportional to the epicentral distance, i.e., the nearer the epicentral distance, the larger the terrain effect. On the other hand, as is shown generally in Fig. 12, the terrain effect shows a negative terrain effect in the high topographical side, but a positive effect on the lower topographical side. These behaviors should be further investigated for more topographical cases and seismic events so that we can obtain more general conclusions on the relation between the topography and co-seismic deformations.
Similarly, results of the percentage of co-seismic deformation in the topographic effect are shown in Fig. 13. The percentage is calculated as presented below.  (1) Therein, i D stands for the absolute value of the co-seismic difference displacement considering the topography effect or not; j max signifies the absolute value of the maximum co-seismic displacement considering the topography effect. Results show that the relative impact is 9.05% of the horizontal direction and 2.95% of the vertical direction. Hsu et al. (2011) also studied the impact of topography on surface displacement using a three-dimensional FEM calculation for the 2005 M w 8.7 Nias-Simeulue, Sumatea earthquake. They found the maximum difference is about 5% and mainly localized near the trench where the topographic gradient is highest. Basically the present conclusions coincide with that of Hsu et al. (2011), although both studies are made for different targets and different topographical models. Thus, we conclude that the terrain effect must be considered in the application of dislocation theory to calculate the co-seismic and post seismic deformation.

SUMMARY AND DISCUSSION
(1) As described in this paper, through an introduction of double-node technology, imposing relative displacement in the fault plane, and making a theoretical value to use the constraint in the region at both ends, a numerical simulation elucidated an arbitrary dip-slip fault combination of dislocation theory and the finite element method. Comparison revealed that examination using the finite element method to calculate the co-seismic displacement is correct. (2) Discussion and simulation of size selection of the computational domain were presented in this paper. Results show that, to ensure the calculate accuracy and boundary effect error when resolve the seismic dislocation deformation problem, the correct size of the area of the finite element model must be 30 times larger than the fault. (3) The calculation method was applicable to the 2008 M 8.0 Wenchuan earthquake. The terrain influenced calculations of co-seismic displacement. Results show that the topographic effect might be readily apparent when applying dislocation theory to calculate the co-seismic and post-seismic deformation. The impact of topography must be considered.
In sum, this paper presents a program and method that are useful for a finite element technique to calculate the coseismic deformation, and to support preliminary research and discussion of the topographic effect. To assess the effects of the Earth's curvature accurately and to present more detailed discussion of topographic effects, a follow-up study will also be conducted to expand these ideas for application to three-dimensional space. Moreover, better exploration of the layered effect of the Earth model, selection of the medium parameters, and consideration of the influence of gravity and other issues should be attempted in future studies, because the layer structure, material parameters, fault slip model and internal interface variations are important in computing precise co-seismic deformations, as reported by Hearn and Bürgmann (2005), Moreno et al. (2009), andSato et al. (2010), respectively.