# Three-Dimensional Hydromechanical Modeling during Shearing by Nonuniform Crust Movement

Authors: Yuqing Zhao; You-Kuan Zhang; Xiuyu Liang
Journal: Geofluids
 (2017)

---

## Abstract

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. In Model 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 are more complex and concentrated along the sides of the fault zone and the magnitudes decreased. The distribution 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.


---

## Body

## 1. Introduction

It is commonly recognized that fluid flow in bedrock regions is mainly controlled by faults and fractures [1–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–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–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 solid-fluid 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.


## 2. 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, Figure1(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 x-, y-, and z-axis are Lx = 1000 m, Ly = 500 m, and Lz = 300 m, respectively. The fault zone is located in the middle of the formation and parallel to y-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 XZ-plane in opposite direction were imposed at the boundaries, 0≤x≤Lx/2 and y=Ly, and Lx/2≤x≤Lx and y=0 (the shaded areas in Figure 1), to simulate the shearing process generated by the nonuniform crust movement [28].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,ux, uy, and uz, is meter and that of the velocity u˙y is mm per year.
(a)
(b)


### 2.1. 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.


#### 2.1.1. Solid Phase

Based on the elastic theory and movement equation, the solid stress is governed by the kinetic equation:(1)∂σij∂xj+ρgi=ρ∂u˙i∂t,where σij is stress tensor [N/L2], xi is coordinate vector [L], ρ=1-nρs+nsρw is the bulk density [M/L3], ρs and ρw are densities of the solid and fluid phase, respectively [M/L3],n is porosity [-], s is saturation [-], gi is gravitational acceleration [L/T2], u˙i is velocity vector [L/T], ui is displacement vector [L], and i=x, y, and z 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:(2)ε˙ij=12∂u˙i∂xj+∂u˙j∂xi,where ε˙ij is strain rate tensor [1/T] and εij is strain tensor [-].Based on the linear poroelasticity theory, the strain-stress relation for small deformation of the porous geological media is commonly described by the following constitutive equation:(3)εij=1E1+vσij-vσkkδij,where E is elasticity modulus [N/L2]. v is Poisson’s ratio [-], σkk=σx+σy+σz, δij is Kronecker delta.


#### 2.1.2. Liquid Phase

The fluid flow in the bedrock is described by Darcy’s law:(4)qi=-kijμk^s∂∂xjP-ρwxjgj,where qi is the specific discharge vector [L/T]; kij is permeability tensor of medium [L2]; k^s is the relative permeability which is a function of saturation s [-]; P is pore pressure [N/L2]; μ is fluid viscosity [M/(L × T)].Fluid mass can be expressed with continuity equation as follows:(5)-∂qi∂xi+qv=∂θ∂t,where qv 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] [-].


#### 2.1.3. Coupled Equation

The change in the fluid content is related to the change in pore pressure (P), saturation (s), and mechanical volumetric strains (ε); the response equation for the pore fluid is formulated as(6)1M∂P∂t+ns∂s∂t=1s∂θ∂t-α∂ε∂t,where M is Biot modulus [N/L2], α is Biot coefficient [-], and ε is volumetric strain [-]. In the saturated case (s=1), if the compressibility of the medium’s particles is negligible (α=1), (6) can be reduced to (7)∂P∂t=M∂θ∂t-∂ε∂t.


#### 2.1.4. 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:(8a)uz=0atz=0m(8b)ux=0atx=0m,x=Lx(8c)u˙y=5mm/yearaty=0,Lx2≤x≤Lx(8d)u˙y=-5mm/yearaty=Ly,0≤x≤Lx2(8e)σij=σowhent=0.The boundary condition (8a) means that the bottom of the model was constrained to have no vertical movements (z-direction) as generally adopted by previous researches [17] but it can move in the y-direction, that is, u˙y≠0. There is no displacement in x-direction (see (8b)) at the boundaries of x=0 and x=Lx. The top boundary (z=Lz) 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˙y = −5.0 mm/year is loaded on the XZ-plane at y = Ly and 0≤x≤Lx/2 (the shaded area on the left of Figure 1) and the velocity of u˙y = 5.0 mm/year at the XZ-plane at y=0 and Lx/2≤x≤Lx (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–31]. σo 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,(9a)hx,y,z,t=h0whent=0(9b)∂hx,y,z,t∂n=0at all boundaries(9c)hx,y,z=zatz=Lzand for Model II(10a)hx,y,z,t=h0whent=0(10b)hx,y,z=Lzaty=0,460m≤x≤540m(10c)hx,y,z=Lzaty=Ly,460m≤x≤540m(10d)∂hx,y,z∂n=0at the other boundaries.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 (y=0 and y=Ly when 460m≤x≤540m) 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/m3 and 1000 kg/m3, 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–35].Table 1
Parameters for Models I and II.
Bulk modulus
Cohesion
Shear modulus
Tension
Friction
Dilation
Hydraulic conductivity
Porosity
(Pa)
(Pa)
(Pa)
(Pa)
(°)
(°)
(m/yr)
Bedrock
7.0
E
+
9
1.0
E
+
7
1.0
E
+
9
1.0
E
+
7
30
0
1.0
E
-
5
0.5
Fault
7.0
E
+
8
1.0
E
+
7
1.0
E
+
8
1.0
E
+
7
15
0
1
0.6


### 2.2. 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 three-dimensional 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 1E-20 and 0.001 yr, respectively.


## 2.1. 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.


### 2.1.1. Solid Phase

Based on the elastic theory and movement equation, the solid stress is governed by the kinetic equation:(1)∂σij∂xj+ρgi=ρ∂u˙i∂t,where σij is stress tensor [N/L2], xi is coordinate vector [L], ρ=1-nρs+nsρw is the bulk density [M/L3], ρs and ρw are densities of the solid and fluid phase, respectively [M/L3],n is porosity [-], s is saturation [-], gi is gravitational acceleration [L/T2], u˙i is velocity vector [L/T], ui is displacement vector [L], and i=x, y, and z 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:(2)ε˙ij=12∂u˙i∂xj+∂u˙j∂xi,where ε˙ij is strain rate tensor [1/T] and εij is strain tensor [-].Based on the linear poroelasticity theory, the strain-stress relation for small deformation of the porous geological media is commonly described by the following constitutive equation:(3)εij=1E1+vσij-vσkkδij,where E is elasticity modulus [N/L2]. v is Poisson’s ratio [-], σkk=σx+σy+σz, δij is Kronecker delta.


### 2.1.2. Liquid Phase

The fluid flow in the bedrock is described by Darcy’s law:(4)qi=-kijμk^s∂∂xjP-ρwxjgj,where qi is the specific discharge vector [L/T]; kij is permeability tensor of medium [L2]; k^s is the relative permeability which is a function of saturation s [-]; P is pore pressure [N/L2]; μ is fluid viscosity [M/(L × T)].Fluid mass can be expressed with continuity equation as follows:(5)-∂qi∂xi+qv=∂θ∂t,where qv 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] [-].


### 2.1.3. Coupled Equation

The change in the fluid content is related to the change in pore pressure (P), saturation (s), and mechanical volumetric strains (ε); the response equation for the pore fluid is formulated as(6)1M∂P∂t+ns∂s∂t=1s∂θ∂t-α∂ε∂t,where M is Biot modulus [N/L2], α is Biot coefficient [-], and ε is volumetric strain [-]. In the saturated case (s=1), if the compressibility of the medium’s particles is negligible (α=1), (6) can be reduced to (7)∂P∂t=M∂θ∂t-∂ε∂t.


### 2.1.4. 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:(8a)uz=0atz=0m(8b)ux=0atx=0m,x=Lx(8c)u˙y=5mm/yearaty=0,Lx2≤x≤Lx(8d)u˙y=-5mm/yearaty=Ly,0≤x≤Lx2(8e)σij=σowhent=0.The boundary condition (8a) means that the bottom of the model was constrained to have no vertical movements (z-direction) as generally adopted by previous researches [17] but it can move in the y-direction, that is, u˙y≠0. There is no displacement in x-direction (see (8b)) at the boundaries of x=0 and x=Lx. The top boundary (z=Lz) 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˙y = −5.0 mm/year is loaded on the XZ-plane at y = Ly and 0≤x≤Lx/2 (the shaded area on the left of Figure 1) and the velocity of u˙y = 5.0 mm/year at the XZ-plane at y=0 and Lx/2≤x≤Lx (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–31]. σo 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,(9a)hx,y,z,t=h0whent=0(9b)∂hx,y,z,t∂n=0at all boundaries(9c)hx,y,z=zatz=Lzand for Model II(10a)hx,y,z,t=h0whent=0(10b)hx,y,z=Lzaty=0,460m≤x≤540m(10c)hx,y,z=Lzaty=Ly,460m≤x≤540m(10d)∂hx,y,z∂n=0at the other boundaries.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 (y=0 and y=Ly when 460m≤x≤540m) 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/m3 and 1000 kg/m3, 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–35].Table 1
Parameters for Models I and II.
Bulk modulus
Cohesion
Shear modulus
Tension
Friction
Dilation
Hydraulic conductivity
Porosity
(Pa)
(Pa)
(Pa)
(Pa)
(°)
(°)
(m/yr)
Bedrock
7.0
E
+
9
1.0
E
+
7
1.0
E
+
9
1.0
E
+
7
30
0
1.0
E
-
5
0.5
Fault
7.0
E
+
8
1.0
E
+
7
1.0
E
+
8
1.0
E
+
7
15
0
1
0.6


## 2.1.1. Solid Phase

Based on the elastic theory and movement equation, the solid stress is governed by the kinetic equation:(1)∂σij∂xj+ρgi=ρ∂u˙i∂t,where σij is stress tensor [N/L2], xi is coordinate vector [L], ρ=1-nρs+nsρw is the bulk density [M/L3], ρs and ρw are densities of the solid and fluid phase, respectively [M/L3],n is porosity [-], s is saturation [-], gi is gravitational acceleration [L/T2], u˙i is velocity vector [L/T], ui is displacement vector [L], and i=x, y, and z 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:(2)ε˙ij=12∂u˙i∂xj+∂u˙j∂xi,where ε˙ij is strain rate tensor [1/T] and εij is strain tensor [-].Based on the linear poroelasticity theory, the strain-stress relation for small deformation of the porous geological media is commonly described by the following constitutive equation:(3)εij=1E1+vσij-vσkkδij,where E is elasticity modulus [N/L2]. v is Poisson’s ratio [-], σkk=σx+σy+σz, δij is Kronecker delta.


## 2.1.2. Liquid Phase

The fluid flow in the bedrock is described by Darcy’s law:(4)qi=-kijμk^s∂∂xjP-ρwxjgj,where qi is the specific discharge vector [L/T]; kij is permeability tensor of medium [L2]; k^s is the relative permeability which is a function of saturation s [-]; P is pore pressure [N/L2]; μ is fluid viscosity [M/(L × T)].Fluid mass can be expressed with continuity equation as follows:(5)-∂qi∂xi+qv=∂θ∂t,where qv 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] [-].


## 2.1.3. Coupled Equation

The change in the fluid content is related to the change in pore pressure (P), saturation (s), and mechanical volumetric strains (ε); the response equation for the pore fluid is formulated as(6)1M∂P∂t+ns∂s∂t=1s∂θ∂t-α∂ε∂t,where M is Biot modulus [N/L2], α is Biot coefficient [-], and ε is volumetric strain [-]. In the saturated case (s=1), if the compressibility of the medium’s particles is negligible (α=1), (6) can be reduced to (7)∂P∂t=M∂θ∂t-∂ε∂t.


## 2.1.4. 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:(8a)uz=0atz=0m(8b)ux=0atx=0m,x=Lx(8c)u˙y=5mm/yearaty=0,Lx2≤x≤Lx(8d)u˙y=-5mm/yearaty=Ly,0≤x≤Lx2(8e)σij=σowhent=0.The boundary condition (8a) means that the bottom of the model was constrained to have no vertical movements (z-direction) as generally adopted by previous researches [17] but it can move in the y-direction, that is, u˙y≠0. There is no displacement in x-direction (see (8b)) at the boundaries of x=0 and x=Lx. The top boundary (z=Lz) 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˙y = −5.0 mm/year is loaded on the XZ-plane at y = Ly and 0≤x≤Lx/2 (the shaded area on the left of Figure 1) and the velocity of u˙y = 5.0 mm/year at the XZ-plane at y=0 and Lx/2≤x≤Lx (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–31]. σo 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,(9a)hx,y,z,t=h0whent=0(9b)∂hx,y,z,t∂n=0at all boundaries(9c)hx,y,z=zatz=Lzand for Model II(10a)hx,y,z,t=h0whent=0(10b)hx,y,z=Lzaty=0,460m≤x≤540m(10c)hx,y,z=Lzaty=Ly,460m≤x≤540m(10d)∂hx,y,z∂n=0at the other boundaries.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 (y=0 and y=Ly when 460m≤x≤540m) 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/m3 and 1000 kg/m3, 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–35].Table 1
Parameters for Models I and II.
Bulk modulus
Cohesion
Shear modulus
Tension
Friction
Dilation
Hydraulic conductivity
Porosity
(Pa)
(Pa)
(Pa)
(Pa)
(°)
(°)
(m/yr)
Bedrock
7.0
E
+
9
1.0
E
+
7
1.0
E
+
9
1.0
E
+
7
30
0
1.0
E
-
5
0.5
Fault
7.0
E
+
8
1.0
E
+
7
1.0
E
+
8
1.0
E
+
7
15
0
1
0.6


## 2.2. 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 three-dimensional 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 1E-20 and 0.001 yr, respectively.


## 3. 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.


### 3.1. Change of the Maximum Principal Stress

The changes of the maximum principal stress (Smax) with time in three X-Z cross sections at Y=0 m, 125 m, and 250 m for Model I were presented in Figure 2(a). The initial (t=0) stress field is mainly controlled by gravity and thus Smax is constant along the horizontal X-direction but varies in the vertical Z-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 X-direction and the fixed vertical displacement at the bottom and free in Y-direction. As time progresses, Smax 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.Figure 2
Contour map of the maximum principle stress (Smax) for Model I at different times. (a) is three vertical cross sections of Y=0, 125 m, and 250 m; (b) is three horizontal cross sections of Z=0, 150 m, and 300 m.
(a)
(b)Att = 2500th year, Smax around the middle of the X-Z cross sections decreases as Y increases because the stress at the boundary (Y=0) is most affected by the velocity applied and this effect is weakened away from the boundary. At t = 10000th year, Smax in the middle of the sections increases as Y increases and became the largest in the cross section of Y=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 Smax at Y=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 Smax at Y=250 m (the graphs in the right column) is symmetric since it is located in the middle in X-Z section.The changes of the horizontal distribution ofSmax with time in three cross sections of Z=0, 150 m, and 300 m are provided in Figure 2(b). The initial Smax (the top row) is constant and positive on the top (Z=300 m) and negative at the bottom (Z=0). As time goes, Smax becomes concentrated in the middle due to the shear stress produced by the boundary conditions. The value of Smax at Z=0 (the bottom-left 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 Smax) are formed around the middle green zone of high compressive stress (the two bottom graphs in the middle column). The distribution of Smax is antisymmetric against the line of Y=250 m in the horizontal cross sections.Similar to Figure2 for Model I, Figure 3 presents the changes of Smax 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 Smax 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 Smax: the maximum value of Smax is reduced from 0.7 × 108 in Model I to 0.25 × 108 in Model II and the minimum is increased from −4.0 × 108 to −2.0 × 108 and to make the stress field more complex around and inside the fault zone.Figure 3
Contour map of the maximum principle stress (Smax) for Model II at different times. (a) is three vertical cross sections of Y=0, 125 m, and 250 m; (b) is three horizontal cross sections of Z=0, 150 m, and 300 m.
(a)
(b)


### 3.2. Change of the Pore Pressure

Corresponding toSmax in Figure 2(a), the changes of the pore pressure (P) with time for Model I in three X-Z vertical cross sections at Y=0 m, 125 m, and 250 m are shown in Figure 4(a). As expected, the distribution of P is mainly controlled by that of Smax. Initially, P is hydrostatic which varies vertically. In earlier times, for example, at t = 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. P at Y=0 (the first graphs in the second row) is positive and largest around the middle and decreases as Y increases. For t = 2500th year, P at Y=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 Y=0 decreases while the negative pressure becomes more negative. As Y increases, the positive pressure decreases while the negative pressure becomes more negative. The distribution of P in these cross sections is asymmetric at Y=0 and 125 m but is symmetric at Y=250 m just like that of Smax in Figure 2.Figure 4
Contour map of the pore pressure (P) for Model I at different times. (a) is three vertical cross sections of Y=0, 125 m, and 250 m; (b) is three horizontal cross sections of Z=0, 150 m, and 300 m.
(a)
(b)Corresponding toSmax in Figure 2(b), the horizontal distribution of P in three X-Y cross sections of Z=0, 150 m, and 300 m at different times is provided in Figure 4(b). The initial P (the top row) is constant horizontally. As time goes, P around the middle either increases to be more positive or decreases to be negative depending on the relative location to the boundary faces. At t = 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 P decreases but the negative P becomes stronger and stronger due to the concentrated characteristics of Smax in this area. The positive pressure is the largest at the bottom (Z=0) and decreases upwards at Z=150 m and 300 m while the negative pressure becomes stronger and stronger upwards. Similar to that of Smax, the distribution of P on the X-Y planes is antisymmetric against the middle line of Y=250 m in the three horizontal cross sections.Two typical cross sections ofY=0 and Z=0 (Figure 5) were selected to display the flow direction at t=500 yr and 5000 yr in order to observe the fluid flow clearly in Model I. Based on the spatial-temporal variations of P in Model I, it is inferred that, initially, the fluid is static since there is no pressure gradient. In early time, for example, at t = 500th year (Figure 5(a)), fluid flows from higher P to lower P in the same height. In later time, fluid flows from positive P (the yellow-red area in Figure 5(b) at t=5000 yr) to negative P (the blue area) as positive pressure is strongly concentrated on the right side of the middle line and negative pressure on the left side.Figure 5
The flow direction in two cross sections (Y=0 and Z=0) of Model I at t = 500 yr and 5000 yr.Corresponding toSmax in Figure 3(a), the changes of P for Model II in three X-Z cross sections at Y=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 x-axis at Y=0 m and Y=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 P at t=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, t=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 Y increases to 125 m, the positive pressure decreases while the negative pressure becomes more negative. At Y=250 m, the P around the fault zone is symmetric and negative (the bottom-right graph in (a)).Figure 6
Contour map of the pore pressure (P) for Model II at different times. (a) is three vertical cross sections of Y=0, 125 m, and 250 m; (b) is three horizontal cross sections of Z=0, 150 m, and 300 m.
(a)
(b)Corresponding toSmax in Figure 3(b), the effect of the fault zone on the horizontal cross sections of Z=0, 150 m, and 300 m is shown in Figure 6(b). The initial P is constant horizontally. As time progresses, P is antisymmetric against the middle line (Y=250 m) of these cross sections. In earlier times, for example, at t = 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 P was affected by the fault zone greatly, it is seen that P 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 Smax concentrated in this area; the positive pressure is stronger concentrated as the negative Smax is concentrated here. The P is the largest at the bottom (Z=0) and decreases upwards at Z=150 m and 300 m.The same cross sections were selected to display the flow direction in Model II (Figure7). 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, t=500 yr. The fluid flows out of the model through the two constant head boundaries (the arrows (a) at t=500 yr) due to the higher P concentrated in the middle of the model. In later time, for example, t=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 t=5000 yr): from positive P (the yellow-red area) at Y=0 to negative P (the blue area) at Y=0 m and from positive P at Y=500 m to negative P at Y=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 t=5000 yr of Figure 7(b)).Figure 7
The flow direction in two cross sections (Y=0 and Z=0) of Model II at t = 500 yr and 5000 yr.It is worth noting that the distribution ofSmax, P 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 ofSmax in two models varies significantly at different times under the opposite velocity fields due to the fault zone, which results in large changes in P. For example, in Model I the initial Smax (i.e., t=0) is constant horizontally and varies vertically, and thus the initial P has similar changes. As time goes, Smax became more and more concentrated in the middle of the cross sections (Figure 2), and P is also concentrated in this area; especially the positive P is concentrated on the right side of the middle line (the yellow-red areas in Figure 4) and the negative P 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 Smax in Model II (Figure 3) are similar to those in Model I except around the fault zone. The initial P is also constant horizontally and varies vertically. As time progresses, P is more and more concentrated around the sides of the fault zone and the positive P (the yellow-red areas in Figure 6) and the negative P (the blue areas) are in opposite directions on the two sides of the fault zone due to the concentrated Smax.


### 3.3. 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. Figure8 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 x1 and x2 are located along the plane of X=500 m, where the largest shear stress is generated by the boundary conditions.Table 2
Locations of the seven observation points.
Observation points
0
x
1
x
2
y
1
y
2
z
1
z
2
x
500
0
250
500
500
500
500
y
250
250
250
0
250
250
250
z
150
150
150
150
150
0
300Figure 8
Locations of the seven observation points.Figure9(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 x1 and x2 increases with time to reach their respective peak values and then decreases at late time. The pore pressure at y1 (the short-dashed green curve) increases earlier and faster than that at other points because y1 is located on the boundary of Y=0 where the velocity field is applied, and, definitely, the pore pressure at y1 is larger than that at y2 and 0. The P at z1 (the short-dashed blue curve) has the largest peak value because z1 is located at the fixed bottom while z2 is at the free top face of the formation. At t=10000 years, the pore pressure at y1, y2, and z1 (the two green and the short-dashed blue curves) stays positive while P at 0 and z2 (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 x1 and x2 (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 Smax and P 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 y2, z1, z2, 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 x1 and x2 that are located away from the fault zone. No pressure change at y1 which is located at the boundary where the head is set to be constant.Figure 9
Changes of the pore pressure (P) with time at the seven observation points for (a) Model I and (b) Model II.
(a)
(b)


## 3.1. Change of the Maximum Principal Stress

The changes of the maximum principal stress (Smax) with time in three X-Z cross sections at Y=0 m, 125 m, and 250 m for Model I were presented in Figure 2(a). The initial (t=0) stress field is mainly controlled by gravity and thus Smax is constant along the horizontal X-direction but varies in the vertical Z-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 X-direction and the fixed vertical displacement at the bottom and free in Y-direction. As time progresses, Smax 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.Figure 2
Contour map of the maximum principle stress (Smax) for Model I at different times. (a) is three vertical cross sections of Y=0, 125 m, and 250 m; (b) is three horizontal cross sections of Z=0, 150 m, and 300 m.
(a)
(b)Att = 2500th year, Smax around the middle of the X-Z cross sections decreases as Y increases because the stress at the boundary (Y=0) is most affected by the velocity applied and this effect is weakened away from the boundary. At t = 10000th year, Smax in the middle of the sections increases as Y increases and became the largest in the cross section of Y=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 Smax at Y=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 Smax at Y=250 m (the graphs in the right column) is symmetric since it is located in the middle in X-Z section.The changes of the horizontal distribution ofSmax with time in three cross sections of Z=0, 150 m, and 300 m are provided in Figure 2(b). The initial Smax (the top row) is constant and positive on the top (Z=300 m) and negative at the bottom (Z=0). As time goes, Smax becomes concentrated in the middle due to the shear stress produced by the boundary conditions. The value of Smax at Z=0 (the bottom-left 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 Smax) are formed around the middle green zone of high compressive stress (the two bottom graphs in the middle column). The distribution of Smax is antisymmetric against the line of Y=250 m in the horizontal cross sections.Similar to Figure2 for Model I, Figure 3 presents the changes of Smax 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 Smax 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 Smax: the maximum value of Smax is reduced from 0.7 × 108 in Model I to 0.25 × 108 in Model II and the minimum is increased from −4.0 × 108 to −2.0 × 108 and to make the stress field more complex around and inside the fault zone.Figure 3
Contour map of the maximum principle stress (Smax) for Model II at different times. (a) is three vertical cross sections of Y=0, 125 m, and 250 m; (b) is three horizontal cross sections of Z=0, 150 m, and 300 m.
(a)
(b)


## 3.2. Change of the Pore Pressure

Corresponding toSmax in Figure 2(a), the changes of the pore pressure (P) with time for Model I in three X-Z vertical cross sections at Y=0 m, 125 m, and 250 m are shown in Figure 4(a). As expected, the distribution of P is mainly controlled by that of Smax. Initially, P is hydrostatic which varies vertically. In earlier times, for example, at t = 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. P at Y=0 (the first graphs in the second row) is positive and largest around the middle and decreases as Y increases. For t = 2500th year, P at Y=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 Y=0 decreases while the negative pressure becomes more negative. As Y increases, the positive pressure decreases while the negative pressure becomes more negative. The distribution of P in these cross sections is asymmetric at Y=0 and 125 m but is symmetric at Y=250 m just like that of Smax in Figure 2.Figure 4
Contour map of the pore pressure (P) for Model I at different times. (a) is three vertical cross sections of Y=0, 125 m, and 250 m; (b) is three horizontal cross sections of Z=0, 150 m, and 300 m.
(a)
(b)Corresponding toSmax in Figure 2(b), the horizontal distribution of P in three X-Y cross sections of Z=0, 150 m, and 300 m at different times is provided in Figure 4(b). The initial P (the top row) is constant horizontally. As time goes, P around the middle either increases to be more positive or decreases to be negative depending on the relative location to the boundary faces. At t = 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 P decreases but the negative P becomes stronger and stronger due to the concentrated characteristics of Smax in this area. The positive pressure is the largest at the bottom (Z=0) and decreases upwards at Z=150 m and 300 m while the negative pressure becomes stronger and stronger upwards. Similar to that of Smax, the distribution of P on the X-Y planes is antisymmetric against the middle line of Y=250 m in the three horizontal cross sections.Two typical cross sections ofY=0 and Z=0 (Figure 5) were selected to display the flow direction at t=500 yr and 5000 yr in order to observe the fluid flow clearly in Model I. Based on the spatial-temporal variations of P in Model I, it is inferred that, initially, the fluid is static since there is no pressure gradient. In early time, for example, at t = 500th year (Figure 5(a)), fluid flows from higher P to lower P in the same height. In later time, fluid flows from positive P (the yellow-red area in Figure 5(b) at t=5000 yr) to negative P (the blue area) as positive pressure is strongly concentrated on the right side of the middle line and negative pressure on the left side.Figure 5
The flow direction in two cross sections (Y=0 and Z=0) of Model I at t = 500 yr and 5000 yr.Corresponding toSmax in Figure 3(a), the changes of P for Model II in three X-Z cross sections at Y=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 x-axis at Y=0 m and Y=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 P at t=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, t=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 Y increases to 125 m, the positive pressure decreases while the negative pressure becomes more negative. At Y=250 m, the P around the fault zone is symmetric and negative (the bottom-right graph in (a)).Figure 6
Contour map of the pore pressure (P) for Model II at different times. (a) is three vertical cross sections of Y=0, 125 m, and 250 m; (b) is three horizontal cross sections of Z=0, 150 m, and 300 m.
(a)
(b)Corresponding toSmax in Figure 3(b), the effect of the fault zone on the horizontal cross sections of Z=0, 150 m, and 300 m is shown in Figure 6(b). The initial P is constant horizontally. As time progresses, P is antisymmetric against the middle line (Y=250 m) of these cross sections. In earlier times, for example, at t = 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 P was affected by the fault zone greatly, it is seen that P 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 Smax concentrated in this area; the positive pressure is stronger concentrated as the negative Smax is concentrated here. The P is the largest at the bottom (Z=0) and decreases upwards at Z=150 m and 300 m.The same cross sections were selected to display the flow direction in Model II (Figure7). 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, t=500 yr. The fluid flows out of the model through the two constant head boundaries (the arrows (a) at t=500 yr) due to the higher P concentrated in the middle of the model. In later time, for example, t=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 t=5000 yr): from positive P (the yellow-red area) at Y=0 to negative P (the blue area) at Y=0 m and from positive P at Y=500 m to negative P at Y=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 t=5000 yr of Figure 7(b)).Figure 7
The flow direction in two cross sections (Y=0 and Z=0) of Model II at t = 500 yr and 5000 yr.It is worth noting that the distribution ofSmax, P 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 ofSmax in two models varies significantly at different times under the opposite velocity fields due to the fault zone, which results in large changes in P. For example, in Model I the initial Smax (i.e., t=0) is constant horizontally and varies vertically, and thus the initial P has similar changes. As time goes, Smax became more and more concentrated in the middle of the cross sections (Figure 2), and P is also concentrated in this area; especially the positive P is concentrated on the right side of the middle line (the yellow-red areas in Figure 4) and the negative P 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 Smax in Model II (Figure 3) are similar to those in Model I except around the fault zone. The initial P is also constant horizontally and varies vertically. As time progresses, P is more and more concentrated around the sides of the fault zone and the positive P (the yellow-red areas in Figure 6) and the negative P (the blue areas) are in opposite directions on the two sides of the fault zone due to the concentrated Smax.


## 3.3. 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. Figure8 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 x1 and x2 are located along the plane of X=500 m, where the largest shear stress is generated by the boundary conditions.Table 2
Locations of the seven observation points.
Observation points
0
x
1
x
2
y
1
y
2
z
1
z
2
x
500
0
250
500
500
500
500
y
250
250
250
0
250
250
250
z
150
150
150
150
150
0
300Figure 8
Locations of the seven observation points.Figure9(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 x1 and x2 increases with time to reach their respective peak values and then decreases at late time. The pore pressure at y1 (the short-dashed green curve) increases earlier and faster than that at other points because y1 is located on the boundary of Y=0 where the velocity field is applied, and, definitely, the pore pressure at y1 is larger than that at y2 and 0. The P at z1 (the short-dashed blue curve) has the largest peak value because z1 is located at the fixed bottom while z2 is at the free top face of the formation. At t=10000 years, the pore pressure at y1, y2, and z1 (the two green and the short-dashed blue curves) stays positive while P at 0 and z2 (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 x1 and x2 (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 Smax and P 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 y2, z1, z2, 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 x1 and x2 that are located away from the fault zone. No pressure change at y1 which is located at the boundary where the head is set to be constant.Figure 9
Changes of the pore pressure (P) with time at the seven observation points for (a) Model I and (b) Model II.
(a)
(b)


## 4. 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 (Smax) and pore pressure (P) in two three-dimensional models were simulated with FLAC3D and the results show that distributions of Smax and P 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,Smax 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 (X=500 m) appears in the middle of the formation which is surrounded by a zone of high tensile stress, especially on the plane of Y=250 m (the red areas in the bottom-right graph in Figure 2(a)).(2)
The distribution ofP is mainly controlled by that of Smax. Initially, P is hydrostatic and constant horizontally but varies vertically. As time goes, a zone with abnormal P appears in the middle of the formation due to the elevated stress in this zone. For Model I, P 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 ofSmax 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 Smax 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 ofP 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 P, 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) ofP and to move the peak time earlier.



---

*Source: 9605313-2017-12-17.xml*