Excess Pore Water Pressure within a Deep-seated Bedding Fault: Understanding of Earthquake Induced Large Landslide Initiation

The 2008 Wenchuan Ms8.0 earthquake triggered a landslide (i.e., Daguangbao landslide) with a volume of near 1.2 billion m3. This landslide involved the Paleozoic carbonate strata and the geological background of the sliding zone is a pre-existing bedding fault with a maximum thickness of 5 m and depth of 400 m in the slope. The bedding fault was the groundwater channel in the slope with saturated rock mass, according to the latest survey. To reveal the initiation mechanism of Daguangbao landslide (hereinafter short for DGB landslide) related to groundwater, a kind of simplified model of hard carbonate rock-cut slope with weak layer has been proposed to generalize bedding fault to weak layer in hard carbonate rock layer. And response characteristics of the model have been simulated by fluid-solid coupling algorithm in the FLAC3D program. The results showed that the hard carbonate rocks above and below the weak layer were different in deformation during strong earthquake, which led to the dynamic behavior of vibratory stamping-tension to the weak layer. And then, the discordant deformation of the two layers resulted in the excess pore water pressure featured with instant amplification and cumulative growth in the saturated water weak layer. So it was inferred that the rapid drop of inside effective stress caused by excess pore water pressure in the weak layer during strong earthquake would give rise to the sudden shearing of stress in the front locked segment of the slope after rapid stress concentration, and consequently trigger the DGB landslide.

3 simulation and establishing real-size model. The fluid-solid coupling numerical method was used in our study and the groundwater conditions before and during the earthquake were taken into consideration.

DGB landslide
DGB landslide involving the Sinian Dengying Formation dolomite with the volume of nearly 1.2 billion m 3 was the largest landslide ( Fig. 1) triggered by Wenchuan Earthquake and even an uncommonly giant landslide in the world. The source area was 2.5 km long and the back scarp was at most 650 m in height, involving clear sedimentary strata from the Sinian to Triassic System. The 1.8 km long sliding face with the horizontal projected area of 0.3 km 2 was exposed at the south flank of the source area. The giant landslide mass moved downward from its original slope, leaving a 2.    Fig.2 Water-saturated bedding fault (a)The tunnel reveals ground water in the bedding fault , (b)Characteristics of water seepage in the bedding fault rock, (c)Sketch of geological cross section of the bedding fault Evidence proving that there was abundant groundwater in slope before the DGB landslide found by Pei et al (2018) mainly included the following aspects. Firstly, the Daguangbao was subject to sub-tropical humid monsoon climate, blessed with plenty rain with the mean annual precipitation of 1260 mm. Secondly, Changshiban Gully, Menkanshi Gully and Heigou Gully developed in landslide area were with flowing water all the year round, and converge at the front part of Huangdongzi Gully. Thirdly, Dayanwo Spring was developed at the head of Changshiban Gully at an elevation of about 2600 m before landslide. The spring belonging to bedrock fissure water has never ever been dried up. Fourthly, Qingwayan Spring developed at the foot of landslide at an elevation of about 1550 m belonged to karst fissure water. These groundwater points have helped us infer the groundwater level. Finally, Huangdongzi Gully was the main gully in the zone with the flow rate in dry seasons of 0.3 m3/s and the maximum flood discharge of 60 m3/s. Further studies showed that the rock mass in bedding fault was saturated [21] and groundwater flew out from the pores of mylonitic texture zone and breccia zone as well as the cracks of the cataclastic rock zone and cracked rock zone (Figs. 2a and b). It was believed that the about 5 cm-thick clay zone at the bottom of bedding fault was an impermeable boundary that provided continuous flow conditions for groundwater in bedding fault (instead of seeping downwards). So the bedding fault was the water channel developed within the slope before the earthquake.

Basic Theory
FLAC 3D program was based on the coupled analysis of power and seepage by explicit finite difference-based completely non-linear method. Wherein, Finn equation, the upgraded version of Mohr-Coulomb equation with dynamic pore water pressure rising program related to increment in plastic volume strain, can simulate the accumulation of pore water pressure in materials and generation of excess pore water pressure (Chen et al., 2009) under dynamic actions. In this paper, Fluid-solid coupling analysis has been made for the DGB landslide by the FLAC 3D program and Finn equation, in which groundwater conditions were taken into consideration.
Assuming that effective stress is ' 0  , the one-dimensional rebound modulus of material is r E , the relational expression between pore water pressure increment u  and plastic volume strain increment vd   under undrained conditions is: As the complete fluid-solid coupling method takes a lot of time, incomplete coupling method was used during calculation by analyzing the relationship between dynamic time scale and spreading time, applying mechanical and fluid disturbance and selecting suitable fluid-solid stiffness ratio (ratio between fluid modulus and solid modulus) according to the method of Chen et al (2009). This was because that mechanical calculation time was pretty shorter than the time for fluid coupling and spreading. And seepage influence can be ignored. So it was considered as undrained condition. However, if mechanical disturbance caused system imbalance, coupling between fluid and mechanical progress shall be made based on the fluid-solid stiffness ratio. If solid stiffness was large but the fluid-solid stiffness ratio was far less than 1, it was called rigid skeleton issue (on the contrary, it was called flexible skeleton issue). As DGB landslide was dolomite landslide showing large solid stiffness, it can be assumed as a rigid skeleton. In this case, complete fluid-solid coupling calculation could be ignored.

Modeling and Monitoring Point Arrangement
Numerical dynamic calculation by establishing 3D geological model of real lithological association took a long time for giant landslide like DGB landslide. The purpose of this paper was to disclose the response characteristics associated with groundwater and possible contributions to landslide initiation during strong earthquake, instead of showing the whole process of DGB landslide initiation, so it was necessary to simplify the geological prototype. As bedding fault was a weak layer in slope, it was directly called as weak layer in following paragraphs.
According to geological section map ( Fig. 1), numerical model was 4600 m in vertical length, with the maximum height of 2500 m, weak layer was 5m in thickness and rock stratum inclined by 18° (  Fig. 3). As DGB landslide was mainly composed of dolomite, the numerical model of landslide was overall considered as dolomite. The weak layer was simplified as porous medium with parameters of breccia zone and its bottom was set as impervious bed involving real clay zone. Hydrological conditions of the lower stratum of week zone was ignored. The groundwater level of the numerical model was defined based on the groundwater level inferred by the spring position on slope before landslide. The area above clay zone of weak layer and below groundwater level was set as saturated zone. The part of numerical model with the most sparse grids was the dolomite at the bottom of weak layer, with the maximum size of grid of about 160 m. And the part with the most dense grids was the weak layer, with the minimum size of grid of about 1 m.
To study the slope instability mainly controlled by bedding fault and groundwater, monitoring points were centrally arranged in and around the weak layer. A total of 15 monitoring points were arranged for the model. Wherein 6 points (No. J1-J6) were in the weak layer at the interval of 600 m. While the space between Point J3 and two adjacent points is 300 m. Point J1 in rear part of weak layer, Point J3 in the middle and Point J6 in the front part were arranged vertically at a vertical interval of 20 m. Among the three points, one was in the lower stratum of weak layer and the other two were in upper stratum, respectively remarked as K1-K9. To get the spatial variation of pore water pressure in the weak layer, monitoring points were divided into three parts. Part A was on the tension-induced cracking boundary of slope, close to the inner side of landslide, Part C was around the locked segment of landslide and Part B, at the position of the largest initial pore water pressure, was between Part A and Part C. Please refer to Fig. 4 for the details of monitoring point arrangement.  (2), refer to Engineering Geology Manual (Edition V) for the standard penetration hammering number 1 N of 30 provided that the weak layer was composed of medium dense-dense sand. Thus the values of C1 and C2 can be obtained, 0.124 and 3.228, respectively. While the elasticity modulus of dolomite in upper and lower hard layers was up to 60000 MPa, flexible constitutive model was adopted. During seepage calculation, it was assumed that material particles were incompressible and slope was a model with the permeability irrelevant to seepage direction. Volume modulus of pure water at room temperature, 2×109 Pa, was taken as the water volume modulus. And it was assumed that initial saturation below water level was 1.0, namely, pores were full of water. Permeability coefficient in soil mechanics could be converted to the permeability coefficient in FLAC3D via Equation (3). 2 6 (m / Pa sec) (cm/ s) 1.02 10 k k Material volume modulus and shearing modulus were obtained by Equations 4 and 5.
Volume modulus： Shearing modulus：  (Fig. 5), it was known that pore pressure ratio (ratio of pore pressure and confining pressure) increased gradually with the duration of dynamic loading and became stable at last. The curve of pore pressure ratio was consistent with the curve of laboratory dynamic triaxial test. The pore pressure ratio of numerical test reached 0.84 when stable, close to that of laboratory test (0.83), which meant that the parameters of the weak layer materials we provided were effective.
After parameters of Table 1 were given, the layout of pore pressure field (Fig. 6) under initial seepage conditions of slope before earthquake was obtained by static equilibrium. The layout of pore pressure field was directly related to the water level depth. The results showed that the range of pore pressure of weak layer in slope before earthquake was in the range of 2 MPa -4 MPa, with the maximum around Monitoring Point J3 of weak layer (4 Mpa). Right part of gully was beyond the scope of this paper.

Loading Scheme
For the purpose of reducing seismic reflection at boundary of the model and simulating semi-infinite space conditions (Fig. 3), free field boundaries were applied around model to absorb seismic wave. Local damping was provided to consume massive energy of seismic wave from model and display the damping value of model again. The coefficient of local damping was 0.016. If there was no damping parameter, model oscillation was made (Fig. 7) by solving a certain number of steps under gravity conditions. Now, the reciprocal of displacement oscillation period was the natural vibration frequency. The natural vibration frequency of model in this study was 0.25 Hz.  Fig.7 Displacement oscillation curve of the model In this paper, dynamic load was applied to the bottom of the model, with acceleration time coarse as the input loading. The three-way ground motion data of Qingping Seismic Station, the nearest seismic station to landslide (about 4.5 km), were selected. Seismic waves in EW direction and NS direction were synthesized to the sliding direction of landslide, being taken as horizontal input seismic wave. The monitored vertical (UD) seismic wave at Qingping Seismic Station was taken as the vertical input seismic wave.
Conversion equation for the synthesis of horizontal earthquake wave in sliding direction of landslide is shown as follow:  7), the high frequency part over 3 Hz was removed and the low frequency part was remained to reduce the requirements of the model grids to the maximum frequency in seismic wave as Wenchuan Earthquake was strong, with significant non-stationary characteristics in nature.
To ensure that there was no continued velocity or residual displacement at the bottom of the model after dynamic calculation was finished, acceleration time course was taken as the input loading of our model after base line correction (Figs. 8c and d).
To guarantee the numerical accuracy of wave propagation, the relationship between model mesh size and input waveform frequency is: Wherein: l  -model mesh size;  -corresponding wavelength of the largest input frequency.

Basic Dynamic Response Characteristics
The maximum pore pressure of the weak layer in the slope before earthquake is located at Monitoring Point J3 of Zone B, which has the most intense dynamic response. So Zone B is taken for an example to describe the basic dynamic response characteristics of the slope.
No matter in vertical or horizontal ground motion conditions, the acceleration response peak of the lower hard layer is smaller than that of the upper one (positive value refers to accelerate upwards, negative one refers to accelerate downwards) (Fig. 9), which means that acceleration amplifies in upper hard layer.
Positive and negative marks of vertical displacement, respectively, refer to moving upwards and downwards and those of horizontal displacement refer to moving to the right (outer side of slope) and to the left (inner side of slope). No matter in vertical or horizontal ground motion conditions, the frequency and peak of the displacement signals of lower hard layer are smaller than those of upper hard layer (Fig. 10). The different displacement response between the upper layer and the lower layer was called discordant displacement in this paper. Both upper and lower hard layers show discordant motion.
No matter in vertical or horizontal ground motion conditions, stress value at monitoring point fluctuates greatly, with local response peak even reaching 3 times the stress value before applying dynamic load (negative value refers to pressure stress and positive one refers to tensile stress). It is observed in the model tests of this paper that the pressure stress or tensile stress of weak layer is greater than that of upper hard layer or lower hard layer (Fig. 11), suggesting that the stress of weak layer is amplified.

Response Characteristics of Pore Water Pressure
No matter in vertical, horizontal, or two-way ground motion conditions, the response peak of pore water pressure in upper hard layer is basically consistent, response intensity of pore water pressure in the weak layer is far greater than the upper hard layer (assuming that the lower hard layer is in anhydrous condition, no pore pressure is found) (Fig. 12). In vertical ground motion condition, pore water pressure in the weak layer cumulatively rises by phases, with the peak of pore pressure response is more than 6 times the upper hard layer at most. In horizontal ground motion condition, pore water pressure of weak layer increases instantaneously. While pore water pressure response is smaller vertical ground motion condition and the maximum peak of pore pressure response is twice the upper hard layer. In two-way ground motion conditions, pore water pressure in weak layer shall be provided with the characteristics of pore water pressure in both vertical and horizontal ground motion conditions. In addition, peak of pore pressure response of the weak layer increases, with a maximum of 32 MPa.
Peak of pore pressure response is mainly in the mid-rear section of the weak layer. Vertical ground motion condition is taken as an example here for the description (Figs. 12a and 13). The pore pressure of the upper hard layer fluctuates slightly. However, the pore pressure of the weak layer increases rapidly under the effect of dynamic loading. In the first 3 s, the pore pressure of weak layer is in an accelerated accumulation stage and reaches to a large peak after 3 s. Then it fluctuates in the way of accumulating by stages. The initial pore pressure of weak layer at Monitoring Point J3 (4 MPa) reaches the peaks (18 MPa,26 MPa,27 MPa) at 9 s, 12 s and 14 s, respectively.
In the FLAC3D numerical simulation software, the limit value of negative pore water pressure in the rock mass is defined as fluid tensile strength. In this paper, we mainly focus on studying how is excess pore water pressure produced in the weak zone under dynamic effect. Negative pore water pressure is not the key point of study. In addition, the tensile strength of water in the rock mass is not taken into consideration. Its value is set at 0. Therefore, the negative pore water pressure is ignored.   Fig.13 Pore water pressure distribution during vertical vibration

Discussion
Pore pressure of the upper hard layers fluctuates slightly, while that of the weak layers fluctuates greatly during strong earthquake. While the rising and dropping of pore pressure are both completed within a short time and able to reach the peaks (Fig. 12) for a brief time, which shows the transient feature of pore pressure response of weak layer.
The working mechanism of excess pore water pressure of weak layer on the ground motion acceleration amplification is still unclear. Analysis is only made on the cause of transient excess pore water pressure in weak layer: (1) When the pore pressure of the weak layer increases transiently, stress in the weak layer is shown as pressure stress (Fig. 14). The upper hard layer moving downwards to the lower hard layer, or the lower hard layer moving upwards to the upper hard layer what happens now makes the weak layer " compressed " . (2) On the contrary, if the pore pressure of the weak layer decreases transiently, stress in the weak layer is shown as tensile stress (Fig. 15). The upper hard layer moving upwards to the lower hard layer, or the lower hard layer moving downwards to the upper hard layer what happens now pulls weak layer apart. So it can be deemed that excess pore water pressure of the weak layer is caused by the incongruent deformation of the upper and lower hard layers. However, it should be noted that pore pressure dropped to 0 during tension process for more than one time, as negative pore pressure cannot be detected by software. In other words, it will be recorded as 0 in the program even though there is a negative pore pressure.
Although pore pressure is of the features of transient rising and dropping, it is observed that pore pressure extrema shows a tendency of cumulative increase (Figs. 12 and 13). This is due to both tensile and pressure stress in weak layer can result in the damage and fracture of the weak layer, the increasing of the fracture connectivity and the decreasing in the proportion of rock bridges. As a result, the weak layer expands and material porosity increases. As the FLAC3D analysis in this paper cannot show the damage and fractures, the change in porosity of weak layer cannot be obtained. While the result can be obtained through reasonable speculation. When the weak layer is compressed, pore and material are also under compression. When the weak layer is pulled, the original pores and damaged space are in tension. Groundwater will flow into pores as they expand, which will lead to a cumulative increase in the water content of weak layer.  As vibration damage and excess pore water pressure of weak layer increase, shearing resistance of weak layer decreases, effective stress transfers to locked segment of slope (Fig. 16). It can be reasonably speculated that this process is likely to make locked segment be suddenly snipped off. Thus landslide is initiated in this way.
(a)Vertical stress distribution cloud map of slope (negative value represents pressure) (b) Cloud map of pore water pressure distribution in slope body Fig.14 Weak layer was compressed results in the accumulation of pore water pressure(11.5 s) (a)Vertical stress distribution cloud map of slope (negative value represents pressure) (b) Cloud map of pore water pressure distribution in slope body Fig.15 Weak layer was tensile results in the dissipation of pore water pressure(4.5 s)

Conclusion
In this paper, the DGB landslide is taken as the geological model and dynamic response value of saturated the weak layer is simulated by FLAC 3D and the theory and method of fluid-solid coupling dynamic analysis. It is found that the incongruent deformation is caused by the inconsistency in the dynamic response "pace" of upper and lower hard layers, which leads to the dynamic behavior of the vibrating stamping-tension of the weak layer. The vibrating stamping is the cause to the transient amplification of stress in the weak layer, which can trigger the excess pore water pressure of weak layer transiently. While incongruent deformation causing damage is the cause to the accumulation of excess pore water pressure in the weak layer. During strong earthquake, the pore water pressure in mid-reear segment of the bedding fault increases cumulatively. Under the response of excess pore pressure of bedding fault, the gravitational component of landslide mass transfers from the bedding fault to the locked segment at the front part of the slope. It gives rise to the consequence that the