#### Abstract

A quasistatic simulation of highly nonlinear problems under fault movements was carried out using the EXPLICIT module of ABAQUS. Combined with the secondary development program of the software, the application of the strain softening Mohr–Coulomb model in the simulation was realized. Free field-fault systems were simulated with two types of fault types (normal and reverse faults), four fault dip angles (45°, 60°, 75°, and 90°), and two kinds of soil (sand and clay). Moreover, the rupture laws and sensitivities of the sand and clay were studied with different soil thicknesses and different fault dip angles in the free field. The results show that the width of the ground zone with obvious deformation, which represents the point of the fault outcrop, the critical displacement of the fault, and the rupture characteristics of the overlying soil are closely related to the fault type and soil parameters. The critical displacement of the reverse fault is larger than that of the normal fault. The width of the ground zone with obvious deformation varies from 0.65 to 1.3 and does not exhibit a regular relationship with the type of soil. Compared with a normal fault, the rupture of a reverse fault is not prone to exposure at the surface.

#### 1. Introduction

Understanding the failure mechanisms of the free field under fault movements can provide a reference to analyse complex structures, indicating that the free field is an important part of engineering design and research. The simulation of the free field is simple, and there are no additional structural interactions with the fault, making it a relatively simple case of fault movement. When a bedrock fault is covered with a certain depth of soil, the dip angle, the displacement of the fault, and the soil parameters can be determined. However, there are still two problems worthy of our attention: (1) whether the fault rupture zone can be exposed at the surface along a fault outcrop and (2) what the deformation of the ground is and what causes significant deformation of the ground caused by a rupture of the overlying soil. These problems in practical applications of engineering and research are particularly important, and thus, many scholars’ attention has been directed toward this issue in recent years.

At present, many scholars have studied the rupture laws of the overlying soil due to fault movements, and most studies adopted numerical simulation methods [1–5]. Previous studies have found that when the thickness of the overlying soil reaches 30 m∼50 m or even up to 75 m, the vertical displacement of the fault will reach 3%∼5% or 7% of the soil thickness, respectively, and the overlying soil layer will rupture. When the thickness of the soil layer is more than 100 m, it is difficult for the overlying soil to rupture. If the overlying soil layer on a reverse fault is in a compressive stress state or if the shear modulus of the soil is high, the soil layer can easily rupture. Before the fault rupture zone develops at the ground surface, the propagation of seismic waves can rupture the ground surface [6]. When the soil layer contains weak layers, it can mitigate the effects of the fault movement on the overlying soil. If the fault displacements are the same, the destruction of a reverse fault will be the largest, and the rupture zone of the overlying soil will finally tilt toward the footwall. The destruction caused by a normal fault would be the second largest, and a strike-slip fault is the least destructive [7]. For the case of overlying soil failure due to vertical movement, as the dip angle becomes larger, a greater vertical displacement is required for the surface to rupture. Due to the influences of inertia forces, the loading rate of the fault movement has a certain effect on the displacement required for the ground to rupture and for the fracture angle of the soil layer. The fault type also has a great influence on the relationship between the ground surface fracture angle and the dilation angle and between the ground surface fracture angle and friction angle [8]. When the overlying soil thickness is the same, as the fault displacement changes, the ground rupture zone will be different. However, the permanent deformation characteristics of the overlying soil are almost the same. If the shear wave velocity is low, the ground rupture zone will be wide. With an increase in the soil thickness, the fracture zone also widens. Compared with a clay layer, the effect of a sand layer is very small [9].

Above all, most previous numerical analyses were concerned with the free field-fault system and were focused on studying the fracture characteristics of the overlying soil. However, few studies have been conducted on obvious ground deformations or the sensitivity of the overlaying soil deformation and rupture with different faults and soils. To further study the response characteristics of the free field under fault movements and understand the influences of fault movements on the overlaying soil, this paper comprehensively considered four key factors through numerical simulations, including the quasistatic analysis, mesh size, strain softening, and material damping. The ABAQUS/EXPLICIT software was used for the simulations and calculations; the simulations used the Mohr–Coulomb yield criterion in consideration of the properties of material strain softening. In addition, this study focused on two questions. The first is whether the fault rupture surface can be exposed at the ground surface and the outcrop position, and the second addresses the ground deformation mode and the zone of obvious ground deformation due to overlying soil fractures. Therefore, the fracture deformation characteristics of the overlying free field under a reverse fault and a normal fault can be obtained.

#### 2. Numerical Simulation Methods and Conditions

##### 2.1. Finite Element Model and Method

Combined with the quasistatic method, the numerical simulation of the EXPLICIT module in the finite element software ABAQUS is used for this numerical simulation. The quasistatic method has been widely used in previous studies to solve problems of fault movement and propagation [10]. This method can reflect the dynamic characteristics of loading to a limited extent and indicate the response of the overlying soil under fault movements. Since the problem involves relatively large displacements, the effect of large deformation is considered to obtain improved simulation results.

The explicit dynamics analysis procedure is based on the implementation of an explicit integration rule together with the use of diagonal or “lumped” element mass matrices. The equations of motion for the body are integrated using the explicit central difference integration rule:where is the velocity, is the acceleration, is the increment number, and is the midincrement value.

A small amount of damping is introduced to control high frequency oscillations. With damping, the stable time increment is given bywhere is the fraction of critical damping in the highest mode.

The explicit integration rule is simple but cannot provide the computational efficiency associated with the explicit dynamics procedure. The explicit procedure requires no iterations and no tangent stiffness matrix. A special treatment of the mean velocities, including and , is required for the initial conditions as well as for certain constraints and presenting the results. For the presentation of the results, the state velocities are stored as a linear interpolation of the mean velocities:

The central difference operator is not self-starting because the value of the mean velocity needs to be defined. The initial values (at time *t* = 0) of the velocity and acceleration are set to zero unless otherwise specified by the user. Therefore, the condition is set as follows:

Substituting (4) into the updated expression for yields the following definition of :

The selected soil constitutive model is the Mohr–Coulomb model, which can realize the numerical simulation of most geotechnical engineering problems and achieve good simulation results. To consider the strain softening characteristics, this paper utilizes the user subroutine USDFLD (VUSDFLD) to supplement the Mohr–Coulomb model. The main purpose is to change the field variables at the material points to alter the properties of the material. The following interface programs can be customized: SUBROUTINE USDFLD(FIELD, STATEV, PNEWDT, DIRECT, T, CELENT, TIME, DTIME, CMNAME, ORNAME, NFIELD, NSTATV, NOEL, NPT, LAYER, KSPT, KSTEP, KINC, NDI, NSHR, COORD, JMAC, JMATYP, MATLAYO, LACCFLA) C INCLUDE 'ABA_PARAM.INC' C CHARACTER^{∗}80 CMNAME, ORNAME CHARACTER^{∗}3 FLGRAY(15) DIMENSION FIELD(NFIELD), STATEV(NSTATV), DIRECT(3,3), T(3,3), TIME(2) DIMENSION ARRAY(15), JARRAY(15), JMAC(^{∗}), JMATYP(^{∗}), COORD(^{∗}) *user coding to define* FIELD *and, if necessary*, STATEV *and* PNEWDT RETURN END SUBROUTINE VUSDFLD(nblock, nstatev, nfieldv, nprops, ndir, nshr, jElemUid, kIntPt, kLayer, kSecPt, stepTime, totalTime, dt, cmname, coordMp, direct, T, charLength, props, stateOld) C INCLUDE 'VABA_PARAM.INC' C CHARACTER^{∗}80 CMNAME, ORNAME CHARACTER^{∗}3 FLGRAY(15) DIMENSION FIELD(NFIELD), STATEV(NSTATV), DIRECT(3,3), T(3,3), TIME(2) DIMENSION ARRAY(15), JARRAY(15), JMAC(^{∗}), JMATYP(^{∗}), COORD(^{∗}) *user coding to define* FIELD *and, if necessary*, STATEV *and* PNEWDT RETURN END

The vertical movement velocity of the fault is calculated and compared with = 1 m/s, 0.5 m/s, 0.1 m/s, and 0.05 m/s. The ground vertical displacement curves and ground tilt curves are shown in Figure 1.

**(a)**

**(b)**

According to Figure 1, the high sliding velocity of the fault will concentrate the strain in the lower part of the overlying soil, while the peak value of the ground tilt is smaller than the conditions with a relatively low velocity. However, when the velocity decreases to 0.1 m/s and continues to be reduced, the plastic strain and ground displacement of the soil will no longer show clear changes. Instead, the time required for the calculation will increase. Therefore, the fault sliding velocity used in this numerical simulation is 0.1 m/s, and the quasistatic method is selected for the simulation.

The ground vertical displacement curves and ground tilt curves with different grid sizes are shown in Figure 2. For the finite element calculation, the grid division directly affects the accuracy of the calculation. The simulation compares the conditions where the size of the grid in the middle of the model is *L* = 0.1*H*, *L* = 0.075*H*, *L* = 0.05*H*, and *L* = 0.025*H*, where *H* is the depth of the soil. As seen from Figure 2, the ground deformation is different due to the different grid sizes. With a decrease in the grid size, the maximum inclination of the ground increases gradually. In addition, the calculation results may be small in the large grid and will negatively impact the understanding of the rupture characteristics of the overlying soil. Therefore, the grid size of *L* = 0.025*H* is adopted in this model.

**(a)**

**(b)**

The overlying soil will produce a relatively large deformation and strain, and the reduction of the soil strength should be considered at this point. If the numerical simulation does not consider the impact of this factor, the results may be different from those in practical situations. In addition, previous scholars have made a good comparison of this result [11–13].

Strain softening and nonstrain softening are both carried out in this simulation. When considering strain softening, the residual friction angle and residual dilation angle are and , respectively, and the critical value of the plastic shear strain is . The plastic strain zones of the two conditions both pass through the overlying soil layer and have uniform distributions. However, when strain softening is not considered, the angle between the rupture zone and the horizontal plane is obviously smaller, and the maximum inclination angle of the ground is lower [14, 15]. The ground displacement in consideration of strain softening is shown in Figure 3.

**(a)**

**(b)**

It can be concluded from the comparison of the plastic strain zone and ground deformation under the two conditions that considering strain softening for this model is necessary to fully estimate the effects of fault movements on the overlying soil.

For the quasistatic problem, the method can greatly reduce the influence of dynamic waves on the model, but the damping of the material cannot be ignored. In view of this, Rayleigh damping is considered in the model [16], and numerical simulations both with and without material damping are performed. When considering material damping, the plastic strain zone breaks through the overlaying soil layer with a uniform distribution. However, when material damping is not considered, not only is the overlying soil layer penetrated by the rupture zone but also a local plastic strain zone appears on the left ground side of the fracture zone. In addition, the rupture of the overlying soil is irregular [17]. Similarly, results from considering the material damping of the ground displacement are shown in Figure 4.

**(a)**

**(b)**

From Figure 4, it can be seen that an upward sharp angle appears on the left side of the vertical displacement curve related to the region of local plastic strain concentration. The ground tilt value fluctuates when material damping is not considered. Therefore, although the quasistatic method is used for the numerical simulation, the damping of the material is still considered.

In this model, the overlying soil is a homogeneous single body. The fault plane is flat and penetrates the bedrock to reach the lower part of the soil. Therefore, all of the attention can be focused on the internal response of the overlying soil due to the fault movements. The structure and parameters of the finite element model are shown in Figure 5. When the fault fracture zone appears at the ground surface, the vertical displacement of the fault is the critical fault displacement, which is recorded as *d*_{0}/*H*. *S*/*H* and *P*/*H* are expressed as the corresponding values for the completion of the fault displacement, namely, *d*/*H* = 5%.

##### 2.2. The Model Grid and Boundary Conditions

The model grid and boundary conditions with a 20 m deep overlying soil layer are shown in Figure 6. The bottom boundary of the model is the interface between the overlying soil layer and the fault assuming that the ground and the interface between the soil and bedrock are both horizontal. To minimize the boundary impact of plastic strain accumulation on the central fault grid, the model width should be approximately 3 to 4 times larger than the height, and thus, the simulation adopts a factor of 4. The central area of the fault should be approximately 1.5 to 2 times larger than the height, and thus, a factor of 2 is used here. The model is composed of plane strain quadrilateral elements. The grid should be as refined and regular as possible to ensure the accuracy of the calculation. Therefore, the grid size of the model’s central zone is 2.5% of the model depth *H*, while the percentage gradually increases from 2.5% to 5% from the middle to the outside of the other parts.

The lower boundary is divided into two parts. The left side of the footwall is established with fixed constraints in the horizontal and vertical directions. A displacement is applied to the right side of the hanging wall equivalent to a displacement applied to the fault. The transformation of the fault dip angle is simulated by adjusting the ratio of the horizontal and vertical displacement components. The horizontal displacement on the left side of the hanging wall is equivalent to that of the bedrock. The right side of the footwall is established with a fixed constraint in the horizontal direction, and both sides of the boundary are unconstrained in the vertical direction. The contact between the overlying soil layer and the bedrock is considered to be fully complete. In addition, since this hypothesis is more reasonable than establishing a rough interface between the rock and soil layers, relative slip at the interface can be completely avoided.

In this case,

for (open) and for (closed), where is the contact pressure between two surfaces at a point and is the interpenetration of the surfaces.

The contact constraint is enforced with a Lagrange multiplier representing the contact pressure in a mixed formulation. The virtual work contribution isand the linearized form of the contribution is

#### 3. Analysis of the Displacement Deformation and Parameters of the Overlying Soil

##### 3.1. Model Parameters

This simulation considers only normal and reverse faults; therefore, strike-slip faults are not studied. The fault vertical dislocation velocity *V* is 0.1 m/s, and the grid size of the model *L* is 0.025*H*. A total of 7 fault dip angles are considered; 3 angles of the normal fault are defined as *β* = 45°, *β* = 60°, and *β* = 75°, and 4 angles of the reverse fault are defined as *β* = 45°, *β* = 60°, *β* = 75°, and *β* = 90°. To compare the normal and reverse fault parameters directly, the fault dip angle is defined as *β* = 180° − *ß*_{reverse fault}, namely, *β* = 135°, 120°, 105°, and 90°. This calculation considers both dry sand and clay, and both soils adopt the Mohr–Coulomb elastoplastic constitutive model. The dry sand considers strain softening, while the clay does not have this characteristic. The elastic modulus *E* is a function of the depth *H* according to *E* = *E*_{0}*H*^{1/2}, where the elastic modulus of the dry sand is *E*_{0} = 20 MPa and that of the clay is *E*_{0} = 5 MPa. The parameters of the two types of soil are shown in Table 1. The thickness of the overlying soil is 20 m in the main research. To study the influence of the thickness on the response of the overlying soil, *H* = 5 m, *H* = 10 m, and *H* = 40 m are also considered.

##### 3.2. Displacement and Deformation Analysis

If the numerical parameters are in accordance with the established model, the simulation results are in good agreement with the actual investigations of fault movements [18]. When the overlying soil thickness is 20 m, the resulting ground deformation curves with different fault displacements and different types of overlying soil (dry sand and clay) are shown in Figures 7–10. Two nondimensional parameters are adopted here, namely, *y*/*H* and *x*/*H*.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

From Figures 7–10, it can be seen that the critical displacement of the reverse fault is larger than that of the normal fault. For the normal fault at 45°, the ground deformations of the sand and clay both show a concave shape. For clay, the width of the ground effect caused by the phenomenon is larger, but the concave shape is not as obvious. This may be due to the presence of a second fracture zone that is nearly symmetrical across the middle of the model, perpendicular to the main fracture zone in the overlying soil. In addition, this phenomenon is eliminated with an increase in the fault dip angle. It can be concluded that when the dip angle of the fault is small, a second fracture zone may appear in the overlying soil in the free field. Therefore, when the dip angle of the fault is small in an actual engineering project, it is necessary to perform a comprehensive analysis according to the site conditions and have an exhaustive understanding of the various situations that may be caused by fault movements. Thus, a small dip angle of the fault is closely related to different parameters of the overlying soil; that is, the extent of the dip angle that can cause a second fracture zone depends on the soil parameters.

The ground incline curves are shown in Figures 11–14, and the curves of the maximum absolute inclination angle of the ground (*dy*/*dx*) are shown in Figure 15.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

As indicated in Figure 15, the maximum inclination of the ground corresponding to sand is larger than that corresponding to clay regardless of the fault dip angle [19]. In the reverse fault, the internal compression of the sand is smaller than that of the clay. Therefore, the response of the overlying soil to the surface deformation is greater due to the fault. However, in the normal fault, the uneasy, tensile clay can reduce the soil stretching and shearing caused by the fault due to its cohesion and tensile strength. Therefore, the greater the strength of the overlying soil, the more severe the response will be with the fault, and this violent reaction will be accompanied by a greater ground deformation.

##### 3.3. Parameter Analysis

With an increase in the fault movement, the maximum inclination of the ground increases gradually; in addition, the increase in the speed is slow at the beginning, for example, when *d*/*H* = 0.5%. When the fracture zone reaches the ground, the increase in the speed of the maximum inclination becomes more rapid. The relationship curves between *d*_{0}/*H* and the fault dip angle with the different soils are shown in Figure 16. The following conclusions can be drawn from these curves. First, under a normal fault (*β* < 90°), *d*_{0}/*H* increases with an increase in the fault dip angle *β* for the different soils with values distributed between 0.2% and 0.8%. Second, under a reverse fault (*β* > 90°), the value of *d*_{0}/*H* increases and reaches 1.2% when *β* = 135°. Obviously, the critical displacement of the fault *d*_{0}/*H* increases with an increase in the fault dip angle *β*, and the critical displacement of the reverse fault is larger than that of the normal fault. In general, the *d*_{0}/*H* value of clay is larger than that of sand, and clay requires a larger value of *d*_{0}/*H* than sand to make the fracture zone appear on the ground.

There are two approaches to determine the location of the fault outcrop. One approach is to find the intersection between the main rupture zone and the ground, and the other is the maximum inclination point of the ground. The shortest horizontal distances from the above two points to the fault are expressed as *S* and *P*, respectively. *S*/*H* and *P*/*H* have the same laws; that is, both of them decrease gradually with an increase in the fault dip angle for the normal fault, and they have an increasing tendency with an increase in the fault for the reverse fault (where the fault dip angle *β* = 180° − *ß*_{reverse fault}). In addition, under a normal fault, the *S*/*H* and *P*/*H* values of sand are less than those of clay, but they are larger than those of clay under a reverse fault. The outcrop point of the fault and the intersection between the fault and soil layer are connected by segments, and the angles between each segment and the horizontal plane (referred to as the “fracture zone inclination”) are shown in Table 2. From Table 2, it can be seen that the fracture zone inclination of sand is higher than that of clay in a normal fault, while the opposite is true in a reverse fault.

The relationship between *L*/*H* and the relative displacement *d*/*H* of the lower bedrock is illustrated in Figure 17. It can be seen that the *L*/*H* value does not exceed the given critical value on the ground before the movement of the bedrock, after which *L*/*H* increases at a rapid rate. However, the change rate of *L*/*H* decreases obviously and tends to be stable when *L*/*H* reaches 1%. It can be seen from the trend line that *L*/*H* is directly proportional to the fault dip angle and that different soils correspond to different values of *L*/*H*. Moreover, when the corresponding *d*_{0}/*H* value is higher, the value of *L*/*H* is higher. This trend illustrates that the magnitude of *L*/*H* is closely related to the accumulated deformation in the soil before the fracture zone reaches the ground; this explains why the *L*/*H* value of clay is larger than that of sand. Therefore, the value of *L*/*H* in the reverse fault is larger than that in the normal fault. This is probably due to the effect of *d*_{0}/*H*, but the difference in *L*/*H* between the reverse and normal faults is smaller than *d*_{0}/*H* because the change in *L*/*H* is momentarily halted when *d*/*H* is larger than 1%.

**(a)**

**(b)**

To illustrate the influences of different soil depths *H* on the various parameters, *H* is set to 5 m, 10 m, 20 m, and 40 m for the calculations. The relationships between *P*/*H*, *S*/*H*, *L*/*H*, and *d*_{0}/*H* and the depth *H* are shown in Figure 18.

**(a)**

**(b)**

**(c)**

**(d)**

From Figure 18, it can be seen that the performance characteristics are the same regardless of the soil type. Namely, the *P*/*H* and *S*/*H* values are nearly uncorrelated with the depth *H* and have some independence. The *L*/*H* and *d*_{0}/*H* values increase with an increase in the depth *H*. With variations in the depth *H*, the changes in *P*/*H* and *S*/*H* are not large and tend to be smooth. Therefore, this method can be used to estimate the location of the fault outcrop point and provide convenience for engineering design and research purposes. The *L*/*H* and *d*_{0}/*H* values increase gradually with a change in the depth; that is, the obvious surface deformation zone increases and the critical displacement required for the fault outcrop increases.

#### 4. Conclusions

The finite element software ABAQUS is utilized to simulate the free field-fault system and to study the fracture mode, ground displacement, and deformation characteristics of the overlying soil in the free field.(1)The direction of a fault in the fracture zone of the overlying soil may deflect or bend with respect to the fault dip. When the dip angle of the normal fault is less than or equal to 45°, a second fracture zone may appear in the overlying soil that may lead to a downward movement of a block of the triangular fracture, which is similar to Coulomb’s Earth pressure theory of soil mechanics. However, this phenomenon will be eliminated with an increase in the fault dip angle.(2)Compared with a normal fault, the fracture zone in a reverse fault cannot easily outcrop at the ground surface; that is, a larger critical fault displacement is required before the outcrop is achieved. Moreover, with an increase in the fault dip angle, the critical fault displacement increases gradually, and the width of the obvious surface deformation zone at the ground surface also increases gradually from the normal fault at 45° to the reverse fault at 45° (*β* = 135°).(3)Under the conditions with the same fault dip angle and displacement, the dynamic response of sand is greater than that of clay. In addition, sand exhibits a larger deformation at the ground surface in close relation to the soil strength. Sand is a compressive and nontensile material, and the compression of sand is less than that of clay due to its high strength. Therefore, the response of the overlying soil surface to the movement of the fault is greater, and the fracture zone dip angle of the sand is smaller than that of clay. In a normal fault, the cohesion and a portion of the tensile strength of clay can reduce the effects of the extension and shearing coincident with the fault movements. Therefore, the fracture zone dip angle of sand is larger than that of clay.(4)The width of the obvious deformation zone at the ground surface and the critical displacement of the fault are directly proportional to the depth of the overlying soil. However, the proportional relationship between the horizontal distance of the two outcrop points to the fault and the soil depth H is not obvious.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was supported by the National Natural Science Foundation of China under Grant no. 41602332 and the Student Innovation and Entrepreneurship Training Program of China under Grant no. 201710615051.