Three-Dimensional Hydromechanical Modeling during Shearing by Nonuniform Crust Movement

Hydromechanical modeling of a geological formation under shearing by the nonuniform crust movement during 10000 years was carried out to investigate the solid stress and pore pressure coupling processes of the formation from the intact to the fractured or faulted. Two three-dimensional numerical models were built and velocities in opposite directions were applied on the boundaries to produce the shearing due to the nonuniform crust movement. The results show that the stress and pore pressure became more and more concentrated in and around the middle of the formation as time progresses. InModel I with no fault, stress and pore pressure are concentrated in the middle of the model during shearing; however, in Model II with a fault zone of weakened mechanical properties, they aremore complex and concentrated along the sides of the fault zone and themagnitudes decreased.Thedistribution of stress determines pore pressure which in turn controls fluid flow. Fluid flow occurs in the middle in Model I but along the sides of the fault zone in Model II. The results of this study improve our understanding of the rock-fluid interaction processes affected by crustal movement and may guide practical investigations in geological formations.


Introduction
It is commonly recognized that fluid flow in bedrock regions is mainly controlled by faults and fractures [1][2][3][4] and that faults and fractures were usually generated by crust movements resulting from the mantle convection based on the theory of plate tectonics [5,6].Existence of faults or fractures significantly affects the flow and transport in geological formation [5,7].It is critically important to understand how the solid stress and fluid pressure evolved in a formation due to crust movement and what is the effect of an existing fault on the evolution and distribution of the stress and pressure in various research areas and engineering applications, such as earthquake prediction and construction of repositories for radioactive waste, dam foundations, underground tunnels, and oil and gas facilities.
The speed of the earth crust movement in a formation is generally nonuniform.Based on the GPS data, for example, Niu et al. [6] pointed out that "the South China block and North China block in eastern China move EES-ward relative to the stable Eurasia.The velocities increase toward the south . ..." The average velocity of the crustal movement to the SEE direction from 1999 to 2011 is about 2.50 mm/year in the Northeast China and about 2.90 mm/year in the same direction in South China [8].The speed difference of 0.40 mm/year seems small but the difference in the accumulative moving distances of the North and South China is huge in geological time.Such a difference in the crustal movement may generate shear stress which in turn creates the strike of faults and fractures along the NWW-SEE direction in eastern China.In fact, Xiao devoted his whole life from 1950s to 2000s to identify the new and active faults and fractures in eastern China for the purpose of locating water-supply wells which commonly have the large capacity of water-supply located in the active faults and fractures [5].He found that most active faults and fractures in eastern China are oriented in the SEE-NWW direction.The objective of this study is to improve our understanding of the coupled processes of the solid stress and fluid pressure in the formations by carrying out a hydromechanical modeling of shallow geological formations under shearing induced by the nonuniform crust movement.
The coupling processes of bedrock deformations and pore pressure changes can be described by the theory of poroelasticity [9][10][11][12][13].Garven investigated the mechanics of flow near a fault on the fold and thrust margin of a basin, especially during active faulting [14].Forster and Evans conducted the experiments and numerical simulations to characterize the permeability structure of thrust zones and found that the regional groundwater flow patterns were controlled by thrust faults and topographies [15].Bredehoeft et al. found that the fluid flow plays an important role in geologic processes and the overpressure affects the development of geological structures [16].Ge and Garven built a two-dimensional numerical model to investigate the effects of the instantaneous compressional tectonics on the regional groundwater flow [17].Saffer [18], Ge and Screaton [19], and Zhou and Hou [20] considered the impacts of the effective stress on fault deformations or other structural developments in the subsurface geological evolution.All above studies considered the case of instantaneous loading rather than continuous loading.On the other hand, most of previous studies on the impacts of faults or fractures on groundwater were either based on the faults or fractures that already exist [21][22][23][24][25] or based on the condition that rock containing microcracking and microvoiding which are the main reasons for deformation and failure of brittle materials like rocks, concrete, and ceramics [26,27].
In this study, we carried out hydromechanical simulations to investigate the coupling processes of solid stress and fluid pressure in a shallow geological formation under a continuous loading caused by the nonuniform crustal movement during 10000 years.Two three-dimensional solidfluid coupling models were built: Model I without a fault zone and Model II with a fault zone.Based on the simulation results, we analyzed the evolution of the solid stress and fluid pressure in the formation during this time period.The results of this study help to improve our understanding of these complex processes.In the following, we will first describe the methodology, then present the results and discussion, and finally draw some conclusions.

Methodology
As mentioned before, the speed of the earth crust movement is uneven and this nonuniform movement may create fractures and/or faults.Groundwater flow in bedrocks is usually controlled by faults and fractures.The coupling processes of the solid stress and pore pressure of bedrocks from the intact to the fractured or faulted under the nonuniform crust movement are very important in controlling the groundwater flow pattern.Therefore, two solid-fluid coupling numerical models were built.One model is for a homogenous formation without a fault zone (Model I, Figure 1(a)) and the other is for a formation with a fault zone (Model II, Figure 1(b)).The dimensions of two models are the same and the lengths along -, -, and -axis are   = 1000 m,   = 500 m, and   = 300 m, respectively.The fault zone is located in the middle of the formation and parallel to -axis with the size of 80 m × 500 m × 300 m (Figure 1(b)).The medium is isotropic and elastoplastic.Two velocity fields normal to the -plane in opposite direction were imposed at the boundaries, 0 ≤  ≤   /2 and  =   , and   /2 ≤  ≤   and  = 0 (the shaded areas in Figure 1), to simulate the shearing process generated by the nonuniform crust movement [28].

Mathematical Model.
The solid phase is governed by the elasticity equations including the kinetic equation, the geometric equation, and the constitutive equation.The liquid phase is governed by a constitutive equation based on Darcy's law.The coupled behavior of solid and fluid is governed by the equilibrium equation between the volumetric strains, pore pressure, and saturation.

Solid Phase.
Based on the elastic theory and movement equation, the solid stress is governed by the kinetic equation: where ,   and   are densities of the solid and fluid phase, respectively , and  = , , and  present three coordinate orientations.Velocities cause changes in strains, and thus the relationship between the velocity and strain can be described with the geometric (compatibility) equation: where ε  is strain rate tensor [1/T] and   is strain tensor [-].Based on the linear poroelasticity theory, the strainstress relation for small deformation of the porous geological media is commonly described by the following constitutive equation: where =   +   +   ,   is Kronecker delta.

Liquid Phase.
The fluid flow in the bedrock is described by Darcy's law: where   is the specific discharge vector  Fluid mass can be expressed with continuity equation as follows: where  V is the volumetric fluid source intensity [1/T] and  is the fluid content (fluid volume per unit volume of porous material) as introduced by Biot [9] [-].

Coupled Equation.
The change in the fluid content is related to the change in pore pressure (), saturation (), and mechanical volumetric strains (); the response equation for the pore fluid is formulated as where and  is volumetric strain [-].In the saturated case ( = 1), if the compressibility of the medium's particles is negligible ( = 1), (6) can be reduced to

Initial Conditions and Boundary Conditions.
The initial and mechanical boundary conditions for deformation of the two model are the same and are set as follows: The boundary condition (8a) means that the bottom of the model was constrained to have no vertical movements (direction) as generally adopted by previous researches [17] but it can move in the -direction, that is, u  ̸ = 0.There is no displacement in -direction (see (8b)) at the boundaries of  = 0 and  =   .The top boundary ( =   ) is the ground surface and is free.In order to study the changes of the stress and pore pressure of the bedrock or the formation of a fault zone under continuous tectonic geostress, the velocity of u  = −5.0mm/year is loaded on the -plane at  =   and 0 ≤  ≤   /2 (the shaded area on the left of Figure 1) and the velocity of u  = 5.0 mm/year at the -plane at  = 0 and   /2 ≤  ≤   (see (8c) and ( 8d)) (the shaded area on the right in Figure 1).The magnitude of the velocity refers to the GPS observed velocities of the crustal movement in the southwest of China [29][30][31].  is initial stress filed that is obtained using the solutions of steady-state kinetic equation ( 1) with the boundary conditions (8a) and (8b).
The initial and boundary conditions for fluid flow of two models are different.For Model I, ℎ (, , , ) = ℎ 0 when  = 0 For Model I, the initial hydraulic head is equal to the static pressure (see (9a)).All boundaries are set as impermeable boundaries in Model I (see (9b)).The solid medium is considered as saturated and unconfined.The pore pressure at the top boundary is equal to the atmospheric condition.For Model II, the boundaries of the fault zone ( = 0 and  =   when 460 m ≤  ≤ 540 m) are set as constant heads in order to investigate the effect of the boundary condition on the groundwater flow.The rest boundaries are all impermeable.
The densities of the formation and groundwater are taken to be 2500 kg/m 3 and 1000 kg/m 3 , respectively.Relative large porosity values of the bedrock are adopted in order to generate optimistic scenarios for the pore pressure redistribution induced by the continuous plate movement.The other material mechanical and physical parameters used in the computation of Models I and II are listed in Table 1 which are based on the previous literatures [32][33][34][35].

Numerical Models.
The two numerical models were constructed with FLAC3D which is a three-dimensional finite difference numerical software for solving problems in engineering mechanics, that is, static and dynamic mechanics and hydraulic effects, including their coupling process [20].According to the previous literatures, such as Papanastasiou and Thiercelin [36], Papanastasiou [37], and Zhou and Hou [20], the material deformation can be simulated using the Mohr-Coulomb constitutive model which was commonly used for simulating crack propagation when the material reached the yield limits.In this study, we coupled the static and dynamic mechanics and hydraulics with the Mohr-Coulomb constitutive model in FLAC3D.
As mentioned above, two models were built with dimension of 1000 m × 500 m × 300 m, and we generated a threedimensional mesh of 50 × 10 × 10 with a total of 6171 nodes for Model I and a mesh of 50 × 10 × 10 (40 × 10 × 10 in the bedrock and 10 × 10 × 10 in the fault zone) with a total of 6171 nodes for Model II.The dynamic damping type is Rayleigh with two relevant parameters  min and  min where  min is set to be 0.05 which is a normal value of the minimum critical damping ratio in geomaterial according to the FLAC3D manual, and  min is 0.122 which was generated by the natural frequency of vibration of two models under the gravity.We used a desktop computer with Intel(R) Core(TM) i7-6700 CPU @ 2.60 GHz 2.60 GHz, 64 GB ram, and 1 TB external storage to run the simulation.It took about 81 hours to run the simulation of Model I and about 45 hours of Model II.The computing efficiency was mainly controlled by the convergence criterion and the size of the time step.The optimal convergence criterion of the maximum unbalanced force and time step we chose is 1 − 20 and 0.001 yr, respectively.

Results and Discussion
The coupled equation (7) for the solid and liquid phases with the initial and boundary conditions (see (8a)-(10d)) for the two models was solved with FLAC3D and the simulation results for the stress and pore pressure are presented and explained as follows.6

Change of the
Geofluids -direction (the first row).It is positive (tensile) at the top and changed to negative (compressive) at the bottom due to the fixed horizontal displacement in -direction and the fixed vertical displacement at the bottom and free in direction.As time progresses,  max became more and more concentrated in the middle of the cross sections due to the shear stress generated by the opposite velocities applied at the two boundary faces.
At  = 2500th year,  max around the middle of the - cross sections decreases as  increases because the stress at the boundary ( = 0) is most affected by the velocity applied and this effect is weakened away from the boundary.At  = 10000th year,  max in the middle of the sections increases as  increases and became the largest in the cross section of  = 250 m (the bottom-right graph) at which a plastic failure zone may appear.At this time a strong tensile stress is formed around the middle of the section (the red area in the bottom-right graph).It is noticed that the distributions of  max at  = 0 and 125 m (the graphs in the left and middle columns) are asymmetric, being mostly compressive on the right half and tensile on the left half (the first graph on second row) while  max at  = 250 m (the graphs in the right column) is symmetric since it is located in the middle in - section.
The changes of the horizontal distribution of  max with time in three cross sections of  = 0, 150 m, and 300 m are provided in Figure 2(b).The initial  max (the top row) is constant and positive on the top ( = 300 m) and negative at the bottom ( = 0).As time goes,  max becomes concentrated in the middle due to the shear stress produced by the boundary conditions.The value of  max at  = 0 (the bottomleft graph) is largest because the bottom boundary is fixed in the vertical direction and decreases upwards.It is noticed that two small red areas of tensile stress (positive  max ) are formed around the middle green zone of high compressive stress (the two bottom graphs in the middle column).The distribution of  max is antisymmetric against the line of  = 250 m in the horizontal cross sections.
Similar to Figure 2 for Model I, Figure 3 presents the changes of  max with time for Model II in the same cross sections.As mentioned previously, the only difference between these two models is that a fault zone (80 m × 500 m × 300 m) with weaker mechanical properties and different hydraulic parameters is contained in Model II in order to investigate the effect of an existing fault on the distribution of the stress and pore pressure under the continuous tectonic movement.The other boundary conditions and parameters are the same with Model I. Overall, the changes of  max in Model II is similar with those in Model I except around the fault zone that has weaker mechanical properties.The effect of the fault zone is to reduce the range of the values of  max : the maximum value of  max is reduced from 0.7 × 10 8 in Model I to 0.25 × 10 8 in Model II and the minimum is increased from −4.0 × 10 8 to −2.0 × 10 8 and to make the stress field more complex around and inside the fault zone.

Change of the Pore Pressure.
Corresponding to  max in Figure 2(a), the changes of the pore pressure () with time for Model I in three - vertical cross sections at  = 0 m, 125 m, and 250 m are shown in Figure 4(a).As expected, the distribution of  is mainly controlled by that of  max .Initially,  is hydrostatic which varies vertically.In earlier times, for example, at  = 500th year, an abnormal positive pressure zone is formed in the middle of the cross sections due to the concentrated compressive stress generated by the boundary conditions. at  = 0 (the first graphs in the second row) is positive and largest around the middle and decreases as  increases.For  = 2500th year,  at  = 0 becomes strongly positive on the right side of the middle line (the yellow-red areas in the first graph on the second row in the left column) and negative on the left side (the blue region) due to the strong compressive stress on the right and tensile stress on the left.As time goes, the positive pressure at  = 0 decreases while the negative pressure becomes more negative.As  increases, the positive pressure decreases while the negative pressure becomes more negative.The distribution of  in these cross sections is asymmetric at  = 0 and 125 m but is symmetric at  = 250 m just like that of  max in Figure 2.
Corresponding to  max in Figure 2(b), the horizontal distribution of  in three - cross sections of  = 0, 150 m, and 300 m at different times is provided in Figure 4(b).The initial  (the top row) is constant horizontally.As time goes,  around the middle either increases to be more positive or decreases to be negative depending on the relative location to the boundary faces.At  = 2500th year, a zone of large positive pressure appears around the middle (yellow-red region) with two small areas of negative pressure (blue region).As time goes, the positive  decreases but the negative  becomes stronger and stronger due to the concentrated characteristics of  max in this area.The positive pressure is the largest at the bottom ( = 0) and decreases upwards at  = 150 m and 300 m while the negative pressure becomes stronger and stronger upwards.Similar to that of  max , the distribution of  on the - planes is antisymmetric against the middle line of  = 250 m in the three horizontal cross sections.
Two typical cross sections of  = 0 and  = 0 (Figure 5) were selected to display the flow direction at  = 500 yr and 5000 yr in order to observe the fluid flow clearly in Model I. Based on the spatial-temporal variations of  in Model I, it is inferred that, initially, the fluid is static since there is no pressure gradient.In early time, for example, at  = 500th year (Figure 5(a)), fluid flows from higher  to lower  in the same height.In later time, fluid flows from positive  (the yellowred area in Figure 5(b) at  = 5000 yr) to negative  (the blue area) as positive pressure is strongly concentrated on the right side of the middle line and negative pressure on the left side.
Corresponding to  max in Figure 3(a), the changes of  for Model II in three - cross sections at  = 0 m, 125 m, and 250 m are shown in Figure 6(a).The boundary conditions for the liquid phase of Model II are the same as those for Model I except for the fact that two constant head boundaries are set for the fault zone between 460 m and 540 m in -axis at  = 0 m and  = 500 m to allow fluid flow to cross the boundaries under the uneven crust movement as the fault zone is an aquifer.It is clearly seen that the initial  at  = 0 (the first row) is constant horizontally and varies vertically.As time progresses, the effect of the fault zone becomes more evident: the positive pressure is concentrated around the fault zone in early times, for example,  = 500 yr (the second row), and, in later times, the higher positive pressure (red region) due to compressive stress is formed on the right side of the fault zone and that of negative pressure (purple region) on the left.As  increases to 125 m, the positive pressure decreases while the negative pressure becomes more negative.At  = 250 m, the  around the fault zone is symmetric and negative (the bottom-right graph in (a)).
Corresponding to  max in Figure 3(b), the effect of the fault zone on the horizontal cross sections of  = 0, 150 m, and 300 m is shown in Figure 6(b).The initial  is constant horizontally.As time progresses,  is antisymmetric against the middle line ( = 250 m) of these cross sections.In earlier times, for example, at  = 500th year, a positive pressure zone is formed in the middle of the cross sections like Model I due to the compressive stress concentrated in this area.In later times, the distribution of  was affected by the fault zone greatly, it is seen that  is mainly concentrated around the sides of the fault zone and is in opposite directions on the two sides of the fault zone: the negative pressure becomes stronger and stronger due to the positive  max concentrated in this area; the positive pressure is stronger concentrated as the negative  max is concentrated here.The  is the largest at the bottom ( = 0) and decreases upwards at  = 150 m and 300 m.
The same cross sections were selected to display the flow direction in Model II (Figure 7).Similar to Model I, the initial fluid is static.As time progresses, the fluid flows from higher to lower pressure at the same height in early time, for example,  = 500 yr.The fluid flows out of the model through the two constant head boundaries (the arrows (a) at  = 500 yr) due to the higher  concentrated in the middle of the model.In later time, for example,  = 5000 yr, fluid flow occurs mainly along the sides of the fault zone and is in opposite directions along the two sides of the fault zone (the arrows (b) at  = 5000 yr): from positive  (the yellow-red area) at  = 0 to negative  (the blue area) at  = 0 m and from positive  at  = 500 m to negative  at  = 500 m; that is, fluid flows from both the matrices and the fault zone to their interfaces.The fluid flows into the fault zone through the constant head boundaries (the arrows at  = 5000 yr of Figure 7(b)).
It is worth noting that the distribution of  max ,  in two models varies significantly under the nonuniform movement of the formation due to the fault zone, and the groundwater flow pattern also varies temporally and spatially.A region with enhanced stress and pore pressure usually occurs around a fault zone, consistent with the theory of mechanics of materials [38].The stress field determines the distribution of the pore pressure which in turn controls the direction of fluid flow.A fault zone may act as a conduit, a barrier, or a combined conduit/barrier system which is essential to the study of faulted flow systems [16,39].Our simulation results show that even with the relatively simple models built in this study a complex stress, pore pressure, and flow fields can be formed in a geological formation under uneven movement of the earth crust.In a formation with no fault, the pore pressure and fluid flow are concentrated in the middle of the formation under different or opposite horizontal velocities.In a formation with a fault zone, however, they are concentrated along the interfaces between the fault zone and the surrounding rock matrices.Geological formations are usually more complicated than the ones simulated and so do the stress, pore pressure, and flow fields.
The distribution of  max in two models varies significantly at different times under the opposite velocity fields due to the fault zone, which results in large changes in .For example, in Model I the initial  max (i.e.,  = 0) is constant horizontally and varies vertically, and thus the initial  has similar changes.As time goes,  max became more and more concentrated in the middle of the cross sections (Figure 2), and  is also concentrated in this area; especially the positive  is concentrated on the right side of the middle line (the yellow-red areas in Figure 4) and the negative  on the left side (the blue areas) due to the strong compressive stress on the right and tensile stress on the left.The changes of  max in Model II (Figure 3) are similar to those in Model I except around the fault zone.The initial  is also constant horizontally and varies vertically.As time progresses,  is more and more concentrated around the sides of the fault zone and the positive  (the yellow-red areas in Figure 6) and the negative  (the blue areas) are in opposite directions on the two sides of the fault zone due to the concentrated  max .

Changes of Pore
Pressure at Selected Points.The changes of the pore pressure with time at seven observation points in the models are depicted in order to better understand the coupling process of rock deformation and fluid flow.Figure 8 shows the locations of the seven observation points (see Table 2 for their coordinators): two along each of the three axes plus the point at the center of the models.It is pointed out that all points except 1 and 2 are located along the plane of  = 500 m, where the largest shear stress is generated by the boundary conditions.Figure 9(a) displays the changes of the pore pressure with time at the seven points for Model I.In general, the pore pressure at all points except 1 and 2 increases with time to reach their respective peak values and then decreases at late time.The pore pressure at 1 (the short-dashed green curve) increases earlier and faster than that at other points because 1 is located on the boundary of  = 0 where the velocity field is applied, and, definitely, the pore pressure at 1 is larger than that at 2 and 0. The  at 1 (the short-dashed blue curve) has the largest peak value because 1 is located at the fixed bottom while 2 is at the free top face of the formation.At  = 10000 years, the pore pressure at 1, 2, and 1 (the two green and the short-dashed blue curves) stays positive while  at 0 and 2 (the black and the long-dashed blue curves) falls below 0, that is, becomes negative, which are due to the tensile stress concentrated area appearing in the middle of the model and becoming larger and larger when time is big enough.The pore pressure at 1 and 2 (the two red lines) which are away from the concentrated shear stress stays almost the same during 10,000 years, indicating that the changes of  max and  mainly occur in the middle of the model.Figure 9(b) shows changes of the pore pressures with time at the same points for Model II.It is seen that the effect of the existing fault zone on the pressure is obvious.First, the magnitudes of the pore pressure at all points are significantly reduced because of the weaker mechanical properties of the fault zone.Secondly, the pore pressure at 2, 1, 2, and 0 reaches their respective peak values earlier and decreases to negative earlier than that at Model I. Little change is observed in the pore pressure at 1 and 2 that are located away from the fault zone.No pressure change at 1 which is located at the boundary where the head is set to be constant.

Summary and Conclusions
In this study we carried out hydromechanical modeling of a geological formation: Model I without a fault zone and Model II with a fault zone, to investigate the coupling processes of the solid strain and pore pressure of the formation under shearing by the nonuniform crust movement during 10000 years.The changes of the maximum principal stress ( max ) and pore pressure () in two three-dimensional models were simulated with FLAC3D and the results show that distributions of  max and  are significantly affected by the existence of a fault zone.Our simulation results show that, even with the relatively simple models built in this study, complex stress, pore pressure, and flow fields can be formed in a geological formation under uneven movement of the earth crust.Specific conclusions drawn in this study are as follows: (1) The initial stress field controlled by gravity is positive (tensile) at the top of the formation and changes to negative (compressive) at the bottom.As time progresses,  max became more and more concentrated in the middle of the formation during shearing generated by continuous loadings on two boundary faces.For Model I, a belt of higher compressive stress (the green areas in Figure 2) along the direction of shearing ( = 500 m) appears in the middle of the formation which is surrounded by a zone of high tensile stress, especially on the plane of  = 250 m (the red areas in the bottom-right graph in Figure 2(a)).
(2) The distribution of  is mainly controlled by that of  max .Initially,  is hydrostatic and constant horizontally but varies vertically.As time goes, a zone with abnormal  appears in the middle of the formation due to the elevated stress in this zone.For Model I,  is strongly positive at the bottom-left corner of the loading face (the red-yellow areas in Figure 4) and strongly negative on the other side and near the surface of the formation.
(3) The changes of  max in Model II is similar to those in Model I except around the fault zone that has weaker mechanical properties.The effect of the fault zone is to reduce the magnitudes of  max and to make the stress field more complex around and inside the fault zone.And, fluid flow in Model II occurs mainly along the sides of the fault zone and is in opposite directions on the two sides of the fault zone.
(4) The distribution of  is significant affected by the fault zone: the abnormal positive and negative pressures (the red area) are mainly distributed along the two interfaces between the fault zone and the surrounding matrices, instead of the middle of the formation in the case of no fault zone.Based on the spatial-temporal variations of , it is inferred that the fluid is initially static.As time goes, fluid starts to move from higher to lower pressure.
(5) The pore pressure in the middle of the formation increases with time to reach their respective peak values and then decreases at late time during shearing.At the rest of the formation, the pore pressure stays relatively constant.The effect of the existing fault zone on the pressure is to reduce the magnitude (including the peak) of  and to move the peak time earlier.

2 Figure 1 :
Figure 1: The two conceptual models: (a) is Model I and (b) is Model II.Two opposite velocity fields were applied to the two boundary faces (the shaded areas) to simulate the shear stress generated from nonuniform crust movement.The unit of the displacement,   ,   , and   , is meter and that of the velocity u  is mm per year.

Figure 2 :
Figure 2: Contour map of the maximum principle stress ( max ) for Model I at different times.(a) is three vertical cross sections of  = 0, 125 m, and 250 m; (b) is three horizontal cross sections of  = 0, 150 m, and 300 m.

Figure 3 :Figure 4 :Figure 5 :
Figure 3: Contour map of the maximum principle stress ( max ) for Model II at different times.(a) is three vertical cross sections of  = 0, 125 m, and 250 m; (b) is three horizontal cross sections of  = 0, 150 m, and 300 m.

Figure 6 :Figure 7 :Figure 8 :
Figure 6: Contour map of the pore pressure () for Model II at different times.(a) is three vertical cross sections of  = 0, 125 m, and 250 m; (b) is three horizontal cross sections of  = 0, 150 m, and 300 m.

Figure 9 :
Figure 9: Changes of the pore pressure () with time at the seven observation points for (a) Model I and (b) Model II.

Table 1 :
Parameters for Models I and II.

Table 2 :
Locations of the seven observation points.