Stress Heterogeneity and Slip Weakening of Faults under Various Stress and Slip

The rock mass of deep underground engineering is in the complex geological environment of high stress, high temperature, and high water pressure. In the process of deep mining and underground space development, the fault-slip seismic source may cause engineering accidents with strong destructive capacity. An in-depth study of fault slip characteristics is very important in the engineering disaster prevention and control. In this paper, a slip model was established based on the finite element software ABAQUS. A total of 20 loading ways are set for various stress and slip, which include the possible slip conditions of fast slip, slow slip, and critical state. By comparing the simulation diagrams and collecting the data of representative grid elements on the loading surface and slip surface, the slip characteristics such as stress heterogeneity under different loads are analyzed. The results show that the increase of slip velocity will make the slip unstable, and the local stress and deformation will become irregular. The spatial stress heterogeneity and the resulting local high working rate will lead to the decrease of the friction strength and the slip weakening. These results can provide some useful suggestions for the research of seismic activities caused by fault slip.


Introduction
Nowadays, due to the saturation of the space capacity, the development of underground space has gradually become the focus of the world [1]. At the same time, with the development of the economy, people's demand for resources is greater than in the past. As shallow resources tend to be exhausted, deep mining has become a new normal for the development of resource-based industries (minerals, gas, oil, etc.) [2]. It is necessary to analyze the underground rock mass whether it is developing underground space or deep mining [3,4]. The deep underground engineering rock mass is in a complex geological environment of high geostress, high temperature, and high water pressure. With the increase in depth and ground stress, mining conditions have become more complex and the rock will be unstable [5], leading to frequent engineering accidents. Among them, the engineering accidents such as rock bust [6] caused by fault-slip seismic sources due to larger magnitude and stronger destructive capacity, resulting to serious casualties and eco-nomic losses [7,8]. The fault-slip source refers to the source type that triggers the surrounding seismic activity due to the relative movement of the fault along the fault plane. In deep mining, the corresponding source types include cavity collapse type and pillar bursts type [9].
The tectonic movement of the crust will cause rock formations to break and form various types of faults. When a fault is released, the kinetic energy may increase the pressure on the surrounding faults and cause vibrations. Such a cyclical process leads to frequent occurrences of various micro and mega earthquakes [10]. In different natural environments, the slip mechanism of fault is different and confusing. Therefore, an in-depth study of fault slip mechanisms is crucial in preventing and controlling earthquakes and engineering disasters. Problems related to faults have attracted much attention in tribology, fracture, and seismic mechanics [11,12]. In the research of fault slip, slip mode [13], slip stability [14,15], the influence of fluid [16,17], slow slip [18][19][20], friction strength [21], etc., have been widely discussed [22]. There are now more experimental methods for small-scale slip experiments [23][24][25]. Centimeter-level rock samples and fault gouges are used in combination with various shear and deformation equipment. However, the sample size is small and the slip process is relatively fast, so the data and conclusions are not enough to represent a wide range of natural faults. Although several large-scale slip experiments have been carried out [26][27][28], the slip velocity can be realized in a certain range, and the change of stress and strain is limited.
The change of normal stress and loading rate will affect the fault response and slip process [29,30]. In this paper, a slip model was established based on the finite element software ABAQUS. Different loading ways are set in consideration of the complex crustal movement and engineering disturbance. The model is loaded according to different stress conditions and loading rates to achieve fast and slow slip. Compared with the existing large-scale fault slip experiment results, the numerical model results are verified to be referential. The stress, strain, displacement, and reaction force diagrams of the slider and the fixed block are obtained. Based on these results, we analyzed the changes in stress, strain, work rate, etc. Finally, the slip characteristics under different stress and slip rate are summarized.

Methods
In order to explore the impact of stress on slip and friction, we conducted a numerical simulation based on the finite element method. Finite element method is widely used to solve the interaction of geometric structure in geological system [31]. Here, we used ABAQUS software to carry on the modeling calculation. The model is two mutually contacting fault blocks. The model refers to the basic information of Indian metagabbros [26] (located in the seismic zone). After consulting the data, the value range of its Young's modulus and Poisson's ratio is obtained, and the size and number of grid cells are set to facilitate the visibility of the output results. The basic parameters of the model values are shown in Table 1.: In order to simulate fault motion similar to the experiment, the loading of the upper fault block is realized by applying consistent vertical pressure and lateral velocity constraint, and the fixation of the lower fault block is realized by applying completely fixed displacement constraints in three directions. The final initial model before calculation is shown in Figure 1.: During the relative motion of the fault, a single fault may slip slowly or quickly. In experimental studies, the slip events are usually divided into slow slip events and fast slip events [28]. In 2015, Yamashita et al. [26] used the secondgeneration large-scale friction device at the National Institute of Earth Sciences and Disaster Prevention (NIED) in Japan (using the vibration table as the driving force to make capable of high loading rate and large slip distance loading) to carry out experiments with meter-sized rock. In order to compare the simulation results with the experimental results, this simulation will be based on the same 12 loading ways as the experiment. In addition, 8 new loading ways are added in order to simulate various slow and fast slip. Finally, A total of 20 loading ways are analyzed and compared together.
The specific loading parameters and model numbers are shown in Table 2.

Results and Discussion
3.1. Normal Stress Change. In order to explore the influence of normal stress changes, we analyze the calculation results when the loading rate is 1 × 10 −2 m/s and the normal stress range is 1.3 MPa to 6.7 MPa. The overall stress distribution of the model can be directly seen through the output Mises stress diagram. It is shown in Figure 2. Mises stress is an equivalent stress based on shear strain energy. It is the equivalent stress of the fourth strength. When Mises stress reaches the yield limit, the material yields. It can be expressed as Mises stress is an equivalent stress; it considers the first, second, and third principal stresses. It can be seen from the figure that the normal stress ranges from 1.3 MPa to 6.7 MPa, and the maximum equivalent stresses are 9.14 MPa, 16.09 MPa, 19.88 MPa, and 41.89 MPa.
With the increase of normal stress, the distribution of equivalent stress is obviously more uniform. The stress value of the upper slider is generally larger, while the stress value of the lower fixed block is generally smaller. For the overall model, the part of the lower fixed block that is not in contact with the upper slider has a smaller stress value, and the side where the upper slider is tangentially loaded has a larger stress value.
For the upper slider, the bottom surface is the slip surface. The tangential loading surface has greater stress on the side of the slip surface. The wear damage of natural rock joints mainly occurs in the local large dip angle area [32]. For this simulation, with a normal stress of 6.7 MPa, the maximum value of the overall model is the corners of the loading surface. The slip surface of the upper slider is the most vulnerable place for small damage. In the actual slip experiment, due to excessive stress and friction, the slip surface will be locally damaged. In the photos taken after the experiment, a series of grooves can be seen 2 Geofluids [26]. The grooves start with very sharp edges and become wider with continuous slip. The grooves are full of gouges and seem to be severely crushed and broken. In order to more clearly see the stress value change during the slip, we select the grid element on the lower right side of the tangential loading surface to calculate the stress components under different normal stresses. The number of the grid element is 30000. The statistics value is shown in Table 3. The slider slips along the positive direction of the z -axis. The values shown in the table are all negative; that is, the unit is subjected to compressive stress, and it increases with the increase of normal stress. Because of the selected grid unit at the coordinate origin, the stress component in the x-axis direction is also negative, and the stress component in the y-axis direction is positive or negative.
The output result of the calculation is LE (true strain/logarithmic strain), and the strain is close to the actual strain.
Generally, the relative linear strain is used in the analysis of small deformation. However, in actual deformation, l 0 is changed to l n , takes many steps. Therefore, the actual strain should be expressed as, Here, l 0 is the initial distance between the two mass points, and l n is the distance after deformation. The strain distribution is shown in Figure 3. The normal stress is from 1.3 MPa to 6.7 MPa, and the maximum strain values are 1:285 × 10 −3 , 1:098 × 10 −3 , 1:024 × 10 −3 , and 9:268 × 10 −4 , respectively. When the loading rate is the same, as the normal stress increases, the strain value will become smaller. When the normal stress is 1.3 MPa, for the overall model (except for the tangential loading surface), the difference of strain is not obvious. As the stress gradually increases, the difference becomes obvious. For strain, when the normal stress is 1.3 and 2.7 MPa, the strain value on the entire tangential loading surface is within a smaller range. When the normal stress is 3.3 MPa, the strain value range begins to be divided into two parts. When the normal stress is 6.7 MPa, it is distributed into four parts. When the normal stress is small, the overall strain of the loading surface is in a larger range, and there is no obvious difference due to the proximity to the upper surface or the slip surface.
Similarly, we select the lower right grid element on the tangential loading surface to calculate the strain components under different normal stresses. The number of the mesh element is 30000, and the statistical values are shown in Table 4. The initial line segment length of the unit in the three directions is less than the length after slip. In addition, on the xy plane, in terms of the right angle between the x-axis and the y-axis, there are larger and smaller values. Comparing the values of strain components under different normal Y Z X Figure 1: Initial model, constraints, and loading mode. The upper part is the slider, and the lower part is the fixed block. The vertical and tangential directions are loaded at the same time to make the slider slip.
3 Geofluids stresses, their magnitudes are not directly affected by changes in normal stress.

Loading Rate Change.
In order to explore the influence of loading rate changes, we analyze the calculation results when the loading rate is 6.7 MPa and the loading rate is 1 × 10 −4 to 3 × 10 −2 m/s. The stress distribution diagram is shown in Figure 4.: It can be seen that when the normal stress is the same, the smaller the loading rate, the more uniform the equivalent stress distribution. When the loading rate is small, the influence of the part of the lower fixed block that is not in contact with the slider is not obvious. As the loading rate increases, the slip distance becomes larger, and there are obvious stress Normal stress gradually increases       Geofluids tangential loading surface under different loading rates is similar, and the maximum value is also at the two corner points on the side close to the slip surface. It shows that the most prone to local damage is also on the slip surface.
We select the grid element on the lower right side of the tangential loading surface to perform the statistics of stress components at different loading rates. The number of the grid element is 30000. The statistical value is shown in Table 5. Here, it may be due to excessive normal stress. The stresses in the x, y, and z directions are all negative; that is, for the element, the forces are all compressive stresses.
The strain distribution of MX4-1, MX4-2, MX4-3, MX4-4, and MX4-5 models whose normal stresses are all 6.7 MPa is shown in Figure 5. When the loading rate is from 1 × 10 −4 to 3 × 10 −2 m/s, the maximum strain values are 8:656 × 10 −5 , 1:382 × 10 −4 , 9:268 × 10 −4 , 6:882 × 10 −4 , and 1:934 × 10 −4 , respectively. It indicates that the change of loading rate does not directly affect the change of strain value. When the loading rate is small, it should be divided into blocks. When the slider slips faster, the local strain is more uneven; that is, the difference in the magnitude of the local shear stress at adjacent positions may also be larger. Therefore, in actual slip experiments, strain gauges need to be attached to the slider and the fixed block to measure the local linear strain. The local shear strain and local shear stress can be calculated based on the linear strain.
For MX4-1, MX4-2, and MX4-5 models, the place with the smallest strain value is in the center of the loading surface. For MX4-3 and MX4-4 models, the place with the smallest strain value is near the corner of the slip surface. This is consistent with the result that the magnitude of the strain value of the overall model is not affected by the loading rate. However, the places with the greatest strain are basically on the tangential loading surface. On the tangential loading surface, we select the grid element on the lower right side to perform the statistics of strain components at different loading rates. The number of the grid element is 30000, and the statistical values are shown in Table 6. For the element model selected here, there is no obvious relationship between the change of strain component and the change of loading rate.

Stress and Strain
Analyses. From the output equivalent stress map, we can get the stress and strain values of each mesh element and node of the model. Here, the maximum values of equivalent stress and strain under all loading ways are counted, as shown in Figure 6.
Under different loading rates, as the normal stress increases, the maximum equivalent stress gradually increases. In the slip experiment, the smaller the normal stress, the smaller the critical shear stress. That is, the easier the shear stress τ reaches the critical shear stress, and the easier the fault will slip. Therefore, the slider is the easiest to slip when the normal stress is 1.3 MPa and the slip rate is 1 × 10 −4 m/s, and it is the most difficult to slip when the normal stress is 6.7 MPa and the slip rate is 3 × 10 −2 m/s. It is not that the greater the slip rate, the greater the stress on the slider and the fixed block, or the smaller the slip rate, the greater the stress value. However, according to the later changes when the normal stress is large, the stress value during slow slip is mostly larger than that during fast slip. This condition may be because the loading time will be longer in order to reach a certain slip distance in slow slip.
As shown in Figure 6, compared with the equivalent stress, the change of the strain value is rather messy. However, the change law at 1 × 10 −3 m/s and 1 × 10 −4 m/s is the same and the strain values are both smaller. In the previous analysis of the strain, the distribution of the strain is unifoof the strain map, when the slip rate is small, the distribution of the strain is uniform. The normal stress of 3.3 MPa and the slip rate of 1 × 10 −2 m/s can be regarded as critical states. When the normal stress is 3.3 Mpa, if the slider slips at a velocity greater than 1 × 10 −2 m/s, it is in an unstable slip state. Therefore, there are obvious inflection points of the stress and strain values at this time, which is not as stable as the slip velocity is small.

Local Shear Stress Calculation.
The existing fault slip simulation based on the finite element method mainly discusses fault rupture, fracture propagation, and the spatio-temporal evolution characteristics of pore pressure and stress field [33][34][35]. However, the local shear stress is not involved in the discussion of stress. In the analysis of the equivalent stress diagram, the local stress in each part of the fault block is different during the slip process. Here, we calculated the local shear stress of the MX4-3 model in the slip direction. Since the front of the model is on the yz plane, we query the value of shear strain to calculate the local shear stress. The calculation formula is τ yz = Gγ yz , where G is the shear modulus; for this simulation, its value is 41200 MPa.
In laboratory experiments, the measurement of local shear strain requires the strain gauge to be attached to the rock in the form of a strain rosette for measurement. There is no way to attach it to the edge of the slip surface. Generally, it will be offset a short distance up or down. According to the position of the strain gage attached in the experiment, onepoint node is taken every 80 mm on the fixed block 20 mm from the slip surface, taken from just below the initial position of the slider. The starting point coordinates is (0, -0.02, 0), the ending point coordinates is (0, -0.02, 1.52). The local shear stress change diagrams of the slider and the fixed block along the z-axis direction are shown in Figure 7. Among all the picked nodes, the local shear stress is the largest at 7 Geofluids 0.32 m from the tangential loading surface and the smallest at 0.56 m from the tangential loading surface. Moreover, along the long side of the slip surface, the value of the local shear stress is always constant. Large changes occur up and down, but from 0.56 m onwards, the overall local shear stress value is smaller than the range of 0.16 m~0.48 m.

Geofluids
Even in a straight line close to and parallel to the slip surface, the value of the shear stress varies greatly. Because the tangential load acts on the slider, the shear stress of the slider near the loading surface will be greater than that of the fixed block. As it gets farther away from the loading surface, the local shear stress of the slider starts to be lower than the fixed block. As far as the upward and downward trends in the middle are concerned, the slider and the fixed block are still similar, but it may be due to the slip of the upper slider. Therefore, the position where the shear stress of the slider and the lower fixed block has the same changing trend will lag slightly behind. Although the local shear stress value on the fixed block is a bit larger, it is more similar to the local shear stress change law on the slip surface than the slider. Therefore, when the local shear stress is measured by the strain gauge, it is better on fixed blocks.

Friction Strength Change.
According to the reaction force value calculated by simulation, we calculate the stress in the tangential loading direction (the area of the tangential loading surface is 0.01 m 2 ). The calculation results are shown in Table 7. In the slip, there will be areas of stress concentration on the fault, and more weak material (gouge) is produced on the fault than in the outer area, resulting in further stress concentration in these areas. The shear stress on the fault is mainly borne by the stress concentration area, which experiences a high work rate (the product of the shear stress and the slip rate per unit area). That is, the higher the work rate, the lower the strength. These areas should weaken quickly and cause a sudden decrease in the macroscopic friction strength [26]. Many laboratory experiments have studied the potential mechanism of reducing friction (weakness) under high slip rate. Among them, the work rate may be a key parameter.
We have drawn a scatter plot of the calculated friction coefficient under different normal stresses as a function of the work rate, as shown in Figure 8. As the work rate increases, the friction coefficient will attenuate to a certain  9 Geofluids extent, the purple-dotted line in the figure. It represents a friction coefficient of 0.75. During slip, high normal stress and high slip velocity lead to high work rate and high productivity of wear materials. The gouge caused by the dislocation of the fault surface is dragged, resulting in further extension of the groove, which will further produce stress heterogeneity. This is a positive feedback effect, which strengthens the heterogeneity of spatial stress.
Due to the high work rate in the local area, the friction strength decreases with the increase of the slip velocity, and the decrease of the friction strength slows down the generation of gouge, which is a negative feedback effect. These two opposite feedback effects compete during fault slip, so friction behavior depends on competition.

Conclusions
In order to study the impact of stress on slip and friction, we conducted a numerical simulation based on the finite element software ABAQUS. On the basis of the 12 loading ways performed in the experiment, 8 new loading ways are set for small slip and large slip. The stress, strain, displacement, etc., of the fault are calculated through numerical simulation. Conducted analysis of deformed graphics and data are as follows: (1) When the slip velocity is the same, with the normal stress increases, the value of the local stress gradually increases and the distribution along the model is more uniform. However, the strain gradually decreases, and the distribution along the model becomes more uneven (2) When the normal stress is the same, with the increase of the slip velocity, the stress value of each element on   (3) Regardless of the stress value or the strain value, the largest value place is on the tangential loading surface. The stress value on the slip surface is relatively large. The slip surface is prone to local damage and grooves during slip (4) During slip, high normal stress and high slip velocity lead to high work rate. At the same time, high productivity of wear materials further leads to stress heterogeneity. The friction strength decreases with the increase of slip rate, reducing the generation of wear material.
Due to the load conditions, the fault slip can be switched between slow slip and fast slip. The increase of slip velocity will cause the slip to become unstable. The local stress and deformation will also become irregular. The spatial heterogeneity of stress and the resulting local high work rate leads to a decrease in friction strength. For example, if a large rock loses its friction strength at a lower work rate, the fault slip may be rapidly weakened in the early stage of the slip period.

Data Availability
Data is available upon request.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.