# Spatiotemporal Evolution Characteristics of Fluid Flow through Large-Scale 3D Rock Mass Containing Filling Joints: An Experimental and Numerical Study

Authors: Qian Yin; Hongwen Jing; Tantan Zhu; Lizhou Wu; Haijian Su; Liyuan Yu
Journal: Geofluids
 (2021)

---

## Abstract

This study analyzes the spatiotemporal evolution characteristics of seepage through a large-scale rock mass containing a filling joint. Firstly, a conceptual model was established to characterize the geomechanical occurrence of a typical water-resistant slab adjacent to a water-bearing structure. Then, a special apparatus was developed to conduct a hydromechanical test of a 3D large-scale rock mass. For a certain boundary stress and inlet water pressure, the pore water pressure in the joint first experiences a dramatic increase before approaching a constant value, and the steady pore water pressure presents a linear decrease along the joint length. A water inrush phenomenon happens as a result of connected flowing channels induced by migration of fillings. Using the finite element of COMSOL multiphysics, the influences of filling joint permeability, matrix permeability, and joint thickness as well as the inlet water pressure on seepage evolution in the jointed rock mass were, respectively, investigated. The pore water pressure increases with all these factors, and the stable pressure values increase with the inlet water pressure but decrease along the joint length. The flow velocity undergoes an increase with both the joint permeability and inlet water pressure but presents constant values independent on the matrix permeability or joint thickness. The water pressure contour planes distributed along the flowing path generally transfer from a “long funnel” shape to a “short funnel” shape before reaching a series of parallel pressure planes perpendicular to the joint direction. By using the genetic algorithm, the coupling influences of these factors on the pore water pressure and flow velocity were investigated, and the decision parameters were optimized. The calculated values show a good agreement with the numerical results, indicating a good prediction of the seepage evolution through the filling joint.


---

## Body

## 1. Introduction

In general, major water conservancy projects under construction and planning in China are concentrated in mountainous canyons at the upper reaches of rivers, and the key point of railway and highway construction is also located in Chinese southwestern regions with unprecedentedly complex geological and hydrological conditions [1–5]. During the construction process of deeply buried long tunnels, geological hazards especially water inrush and mud outburst often occur as a result of large burial depth, high water pressure, high geostress, and strong disturbance [6–8]. Additionally, widely distributed disaster-causing structures such as joints, complicated fracture networks, and underground water-bearing structures in surrounding rocks make it more prone to instability of structural plane and groundwater, which would then produce engineering disasters [9–15]. For example, due to complicated hydrogeological conditions along the diversion tunnel of the Jinping-II Hydropower Station (with the deepest buried depth of 2500 m and the highest water pressure of 10 MPa) in China, during the excavation process, the maximum groundwater discharge is 7.3 m3/s, and 12 points with an extra-large water inflow rate exceeding 1.0 m3/s are monitored (Figures 1(a) and 1(b)). Besides, a total of 19 large-scale water inrush accidents occurred in the construction process of the Yesanguan Tunnel of Yiwan Railway, China, with the largest water bursting discharge of approximately 1.5×105m3/h, unfortunately leading to 15 deaths, major economic losses, a construction delay of 2 years, and serious destruction of ecological environment (Figures 1(c) and 1(d)). Therefore, there is no doubt that water inrush geohazards pose a serious threat to the performance and safety of underground excavation projects [16–19].Figure 1
Typical geological disasters of water burst in underground engineering.
(a)
(b)
(c)
(d)In recent years, researches focused on the water inrush mechanisms and prediction and early warning of geological hazards, and disaster prevention and control measures of practical engineering have been extensively performed [1, 4, 9, 18, 20–25]. A lot of interesting findings have been identified. However, few studies have so far given insight into the seepage spatiotemporal evolution characteristics in a water-resistant slab. Thus, in order to evaluate the process of seepage evaluation and water inrush more intuitively, a generalized model was established (Figure 2). A water-bearing structure, taking a karst cave as an example, with a certain hydraulic pressure in the joints or rock matrix, is in front of the tunneling direction. Due to irreversible perturbations during tunnel excavation, as well as the influence of water pressure, fractures are initiated, propagated, and coalesced to generate complex fracture networks, which would then form abundant fluid flow channels and entail significant changes in the permeability of the water-resistant slab [11, 26–28]. Additionally, both original and induced fractures/joints in the rock mass are usually filled with weathered debris of rocks, later intrusions of mud, and other materials, playing a vital role in seepage evolution and overall permeability of the rock mass [14, 15, 29, 30], which has rarely been reported.Figure 2
Schematic diagram of water inrush disaster in a typical underground project.Thus, in this study, a simplified sketch model describing the geomechanical state of a water-resistant slab was established, as shown in Figure2. Here, the water-resistant slab is ideally treated as a cuboid model applied with a boundary stress of σy and σz in the y- and z-directions in the waterproof strata, a water pressure P on the boundary neighboring the water-bearing structure, and the water pressure on the boundary adjacent to the free face of the tunnel excavation that is assumed to be zero [9, 18, 26]. Based on the model established above, a large-scale three-dimensional physical simulation testing system (300×300×900mm×2) for water inrush disasters in a rock mass was self-developed, and experimental investigation on quantitative analysis of seepage evolution and water inrush process in a rock mass containing a single filling joint was conducted. Then, a number of numerical simulations on the spatiotemporal evolution of pore water pressure and flow velocity, as well as the flow streamline distribution were conducted by employing the finite element method of COMSOL multiphysics, with respect to various factors of joint permeability Kj, matrix permeability Km, joint thickness b, and the inlet hydraulic pressure P. Finally, a genetic algorithm was utilized to propose a prediction model identifying the coupling influences of those factors on the flowing characteristics of the filling joint.


## 2. 3D Large-Scale Seepage Test of Rock Mass Containing a Filling Joint

### 2.1. Development of the Testing System

The special apparatus for a three-dimensional large-scale seepage test of a rock mass is shown in Figure3. This system is mainly equipped with the following four units:
(1)
The main bearing structure of the model. The main body structure consists of a constraint framework and a bearing pedestal. The constraint framework is characterized by a bearing structure with the shape of a square hole inside a circle, and is intensively welded by steel plates with high qualities. Three stiffeners are separately added in each direction in order to improve both the strength and stiffness. The internal bearing structure has a cuboid shape with an allowable model size of 300×300×900mm×2 (Figures 3(a)–3(c))(2)
A boundary stress loading system. The stress loading system consists of a piston-type uniform pressure loader, a servo oil-source system, and a control system. The boundary stress in the y- and z-directions can be individually imposed on a rock mass, with a pressure range of 0–6.0 MPa and a loading accuracy of ±2.0%. The split piston uniform pressure loader can achieve a better uniform load distribution (Figures 3(b) and 3(c))(3)
A gas-liquid composite water supply system. The hydraulic system consists of a group of high-pressure nitrogen tanks and a water tank that can supply a constant water pressure (Figure 3(d)). Nine high-pressure nitrogen cylinders, which can be individually switched on/off and can supply a gas pressure of up to 9.0 MPa for each cylinder, are laid out in parallel. The water tank has a bulk volume of 1.5 m3, a designed maximum hydraulic pressure of 2.0 MPa, and a pressure accuracy of ±5.0%. During the test, water is fed into the model through the water tank that can supply a constant water head at any moment driven by the gas-liquid composite water supply system(4)
A data acquisition system and other auxiliary equipment. The data acquisition system consists of a data acquisition instrument, a sensor power supply module, a circuit board, and a computer (Figure 3(e)). The acquisition instrument is a 32-channel intelligent network distributed data acquisition and processing analyzer, the INV306N-6260, with a maximum sampling rate of 204.8 kHz, which can be utilized to realize the real-time data acquisition, wave display, and spectrum analysisFigure 3
Large-scale three-dimensional fluid flow testing system.


### 2.2. Experimental Materials and Large-Scale Model Preparation

During the test, cement and quartz sand were chosen as the main components to fabricate the matrix with a mass ratio of cement (C32.5 typical Portland cement), quartz sand (20–40 M), and water of 1 : 0.3 : 0.4. The joint was prefabricated and filled with dry and clean quartz sand. First, samples of experimental materials were produced (Figure4(a)) to conduct uniaxial compression (Figure 4(b)), Brazilian splitting, and variable angle shear tests (Figure 4(c)) to determine the basic mechanical properties of the cement mortar. The results indicate that after water curing of 7 days, average unit weight, uniaxial compressive strength, tensile strength, elastic modulus, cohesive force, and internal friction angle of the experimental materials are 19.66 kN/m3, 2.85 MPa, 0.29 MPa, 0.32 GPa, 0.61 MPa, and 31.22°, respectively (Figure 4(d)). Additionally, with an increasing confining pressure from 1.0 to 4.0 MPa, permeability of the samples decreases from 2.406×10−10 to 1.210×10−10m2.Figure 4
Experimental materials and basic mechanical tests.
(a)
Sample preparation
(b)
Uniaxial compression test
(c)
Brazilian splitting and variable angle shear tests
(d)
Stress-strain curves in the uniaxial compression, Brazilian splitting, and angle shear testsThe large-scale model containing a filling joint was poured by utilizing the method of layer filling and compaction, as follows:(1)
Assemble the pouring mould (Figure5(a)), clean the inner surface, and evenly apply a layer of lubricating oil. Prepare the water inlet device (Figure 5(b)). After accurate weighing (Figure 5(c)) and stirring (Figure 5(d)), the experimental materials were poured into the mould laying in three times, respectively, with uniform vibration (Figure 5(e)). The height of each filled layer was 10 mm(2)
During the pouring process, the joint with a length of 900 mm, width of 100 mm, and thickness of 10 mm was prefabricated. Thus, the water inlet device was embedded in the middle position of the model in the second layer (Figure5(f)). After the joint was produced (Figure 5(g)), clean quartz sand was filled (Figure 5(h)), and five high-precision WMS-51 micro water pressure sensors with a pressure range of 0–2.5 MPa and an accuracy of 0.5% were embedded in the joint along a length-wise direction with a spacing of 160 mm (Figure 5(i))Figure 5
A typical process of model pouring, sealing, and installation.


### 2.3. Sealing and Installation of the Large-Scale Model

After completion of the model pouring, the model was made to stand up in a length-wise direction, demoulded, and cured for 7 days. Then, a layer of epoxy resin was evenly applied on the model surface (Figure5(j)) and a layer of crack-resistant seam belt was timely wrapped to bond together with the epoxy resin (Figure 5(k)). After the epoxy resin was solidified after 24 h, a layer of polyurethane was applied on the crack-resistant seam belt (Figure 5(l)). When the polyurethane was solidified after 48 h, the above steps were repeated once again. Finally, the model surface was cleaned, especially the bottom position at which some waste sealing materials may be piled up (Figure 5(m)), and a knife was utilized to cut through the waterproof boundary at the position of the embedded water inlet device, from which the model can be connected with the hydraulic system (Figure 5(n)).The sealed rock mass containing a single filling joint was then placed on a bearing platform with the connector of the water inlet device upwards, and the model was carefully lifted by pushing and the device was precisely positioned using a hoist (Figure5(o)). Then, the model was pushed into the test system, the end of the model was aligned with loading system boundary (Figure 5(p)), and the water inlet device was connected with the gas-liquid composite water supply system using a rubber pipe (Figure 5(q)). Finally, the lead wires of the water pressure sensors were connected with the data acquisition system (Figure 5(r)), debugged, and the initial values cleared.


### 2.4. Test Procedure, Results, and Discussion

The schematic diagram of the rock mass containing a filling joint, gas-liquid composite water supply system, water inlet device, and embedding arrangement of the water pressure sensors is displayed in Figure6. Note that in the experiment, we just focused on half of the model, with the joint (the length of 900 mm) located at the exact central position of the matrix. In addition, a gravel filter was arranged between the water inlet device and the joint to achieve a uniform distribution of water pressure and prevent the water inlet device being blocked. By using a hydraulic access, variable but uniform inlet water pressures can be applied to the rock mass.Figure 6
Schematic diagram of the water supply system and embedding location of the water pressure sensors.During the test, both boundary stresses in they- and z-directions were constant at 5 MPa, and the inlet water pressure P was set to be 0.3 MPa. Then, variations in the pore water pressure p with time (t=0‐2500s) and the joint length (L=60, 220, 380, 540, 700 mm) can be analyzed. Figure 7(a) presents the spatiotemporal evolution characteristics of p. It can be seen that the pore water pressure at all the five positions through the filling joint shows an increase with t. At the beginning, p changes significantly with t. However, as t continuously increases, the increasing rate steadily diminishes and p gradually approaches a constant value. The variation of p as a function of t can be well described using an exponential function. In addition, the variation rate of p at various positions along the joint is different, indicating a larger value for the position close to the water inlet. For t=2500s, the constant p values for the measure points E_#1, E_#2, E_#3, E_#4, and E_#5 are 0.2767, 0.2097, 0.1450, 0.0815, and 0.0195 MPa, respectively, showing a decrease of 92.95% for the water pore pressure along the joint length. Besides, at a small t (e.g., t=100, 300 s), the decrease in p with L presents a nonlinear trend. With increasing t (t=1000s), the variation trend of p versus L gradually transfers from nonlinearity to linearity (Figure 7(b)).Figure 7
Experimental and numerical evolution characteristics ofp.
(a)
p
~
t
(b)
p
~
LAs the inlet water pressure is continually applied, the filling materials are dissolved in water and carried away by water with time. Connected fluid flowing channels are gradually produced in the filling joint, and a water inrush phenomenon from the joint occurs at the model boundary, as shown in Figure8. Besides, this is due to the fact that the flow capacity of the matrix is several orders of magnitude smaller than that of the filling joint. The matrix surface of the model is just wet, which indicates slight leakage rather than significant flow streamlines.Figure 8
Water seepage characteristics at the model boundary.


## 2.1. Development of the Testing System

The special apparatus for a three-dimensional large-scale seepage test of a rock mass is shown in Figure3. This system is mainly equipped with the following four units:
(1)
The main bearing structure of the model. The main body structure consists of a constraint framework and a bearing pedestal. The constraint framework is characterized by a bearing structure with the shape of a square hole inside a circle, and is intensively welded by steel plates with high qualities. Three stiffeners are separately added in each direction in order to improve both the strength and stiffness. The internal bearing structure has a cuboid shape with an allowable model size of 300×300×900mm×2 (Figures 3(a)–3(c))(2)
A boundary stress loading system. The stress loading system consists of a piston-type uniform pressure loader, a servo oil-source system, and a control system. The boundary stress in the y- and z-directions can be individually imposed on a rock mass, with a pressure range of 0–6.0 MPa and a loading accuracy of ±2.0%. The split piston uniform pressure loader can achieve a better uniform load distribution (Figures 3(b) and 3(c))(3)
A gas-liquid composite water supply system. The hydraulic system consists of a group of high-pressure nitrogen tanks and a water tank that can supply a constant water pressure (Figure 3(d)). Nine high-pressure nitrogen cylinders, which can be individually switched on/off and can supply a gas pressure of up to 9.0 MPa for each cylinder, are laid out in parallel. The water tank has a bulk volume of 1.5 m3, a designed maximum hydraulic pressure of 2.0 MPa, and a pressure accuracy of ±5.0%. During the test, water is fed into the model through the water tank that can supply a constant water head at any moment driven by the gas-liquid composite water supply system(4)
A data acquisition system and other auxiliary equipment. The data acquisition system consists of a data acquisition instrument, a sensor power supply module, a circuit board, and a computer (Figure 3(e)). The acquisition instrument is a 32-channel intelligent network distributed data acquisition and processing analyzer, the INV306N-6260, with a maximum sampling rate of 204.8 kHz, which can be utilized to realize the real-time data acquisition, wave display, and spectrum analysisFigure 3
Large-scale three-dimensional fluid flow testing system.


## 2.2. Experimental Materials and Large-Scale Model Preparation

During the test, cement and quartz sand were chosen as the main components to fabricate the matrix with a mass ratio of cement (C32.5 typical Portland cement), quartz sand (20–40 M), and water of 1 : 0.3 : 0.4. The joint was prefabricated and filled with dry and clean quartz sand. First, samples of experimental materials were produced (Figure4(a)) to conduct uniaxial compression (Figure 4(b)), Brazilian splitting, and variable angle shear tests (Figure 4(c)) to determine the basic mechanical properties of the cement mortar. The results indicate that after water curing of 7 days, average unit weight, uniaxial compressive strength, tensile strength, elastic modulus, cohesive force, and internal friction angle of the experimental materials are 19.66 kN/m3, 2.85 MPa, 0.29 MPa, 0.32 GPa, 0.61 MPa, and 31.22°, respectively (Figure 4(d)). Additionally, with an increasing confining pressure from 1.0 to 4.0 MPa, permeability of the samples decreases from 2.406×10−10 to 1.210×10−10m2.Figure 4
Experimental materials and basic mechanical tests.
(a)
Sample preparation
(b)
Uniaxial compression test
(c)
Brazilian splitting and variable angle shear tests
(d)
Stress-strain curves in the uniaxial compression, Brazilian splitting, and angle shear testsThe large-scale model containing a filling joint was poured by utilizing the method of layer filling and compaction, as follows:(1)
Assemble the pouring mould (Figure5(a)), clean the inner surface, and evenly apply a layer of lubricating oil. Prepare the water inlet device (Figure 5(b)). After accurate weighing (Figure 5(c)) and stirring (Figure 5(d)), the experimental materials were poured into the mould laying in three times, respectively, with uniform vibration (Figure 5(e)). The height of each filled layer was 10 mm(2)
During the pouring process, the joint with a length of 900 mm, width of 100 mm, and thickness of 10 mm was prefabricated. Thus, the water inlet device was embedded in the middle position of the model in the second layer (Figure5(f)). After the joint was produced (Figure 5(g)), clean quartz sand was filled (Figure 5(h)), and five high-precision WMS-51 micro water pressure sensors with a pressure range of 0–2.5 MPa and an accuracy of 0.5% were embedded in the joint along a length-wise direction with a spacing of 160 mm (Figure 5(i))Figure 5
A typical process of model pouring, sealing, and installation.


## 2.3. Sealing and Installation of the Large-Scale Model

After completion of the model pouring, the model was made to stand up in a length-wise direction, demoulded, and cured for 7 days. Then, a layer of epoxy resin was evenly applied on the model surface (Figure5(j)) and a layer of crack-resistant seam belt was timely wrapped to bond together with the epoxy resin (Figure 5(k)). After the epoxy resin was solidified after 24 h, a layer of polyurethane was applied on the crack-resistant seam belt (Figure 5(l)). When the polyurethane was solidified after 48 h, the above steps were repeated once again. Finally, the model surface was cleaned, especially the bottom position at which some waste sealing materials may be piled up (Figure 5(m)), and a knife was utilized to cut through the waterproof boundary at the position of the embedded water inlet device, from which the model can be connected with the hydraulic system (Figure 5(n)).The sealed rock mass containing a single filling joint was then placed on a bearing platform with the connector of the water inlet device upwards, and the model was carefully lifted by pushing and the device was precisely positioned using a hoist (Figure5(o)). Then, the model was pushed into the test system, the end of the model was aligned with loading system boundary (Figure 5(p)), and the water inlet device was connected with the gas-liquid composite water supply system using a rubber pipe (Figure 5(q)). Finally, the lead wires of the water pressure sensors were connected with the data acquisition system (Figure 5(r)), debugged, and the initial values cleared.


## 2.4. Test Procedure, Results, and Discussion

The schematic diagram of the rock mass containing a filling joint, gas-liquid composite water supply system, water inlet device, and embedding arrangement of the water pressure sensors is displayed in Figure6. Note that in the experiment, we just focused on half of the model, with the joint (the length of 900 mm) located at the exact central position of the matrix. In addition, a gravel filter was arranged between the water inlet device and the joint to achieve a uniform distribution of water pressure and prevent the water inlet device being blocked. By using a hydraulic access, variable but uniform inlet water pressures can be applied to the rock mass.Figure 6
Schematic diagram of the water supply system and embedding location of the water pressure sensors.During the test, both boundary stresses in they- and z-directions were constant at 5 MPa, and the inlet water pressure P was set to be 0.3 MPa. Then, variations in the pore water pressure p with time (t=0‐2500s) and the joint length (L=60, 220, 380, 540, 700 mm) can be analyzed. Figure 7(a) presents the spatiotemporal evolution characteristics of p. It can be seen that the pore water pressure at all the five positions through the filling joint shows an increase with t. At the beginning, p changes significantly with t. However, as t continuously increases, the increasing rate steadily diminishes and p gradually approaches a constant value. The variation of p as a function of t can be well described using an exponential function. In addition, the variation rate of p at various positions along the joint is different, indicating a larger value for the position close to the water inlet. For t=2500s, the constant p values for the measure points E_#1, E_#2, E_#3, E_#4, and E_#5 are 0.2767, 0.2097, 0.1450, 0.0815, and 0.0195 MPa, respectively, showing a decrease of 92.95% for the water pore pressure along the joint length. Besides, at a small t (e.g., t=100, 300 s), the decrease in p with L presents a nonlinear trend. With increasing t (t=1000s), the variation trend of p versus L gradually transfers from nonlinearity to linearity (Figure 7(b)).Figure 7
Experimental and numerical evolution characteristics ofp.
(a)
p
~
t
(b)
p
~
LAs the inlet water pressure is continually applied, the filling materials are dissolved in water and carried away by water with time. Connected fluid flowing channels are gradually produced in the filling joint, and a water inrush phenomenon from the joint occurs at the model boundary, as shown in Figure8. Besides, this is due to the fact that the flow capacity of the matrix is several orders of magnitude smaller than that of the filling joint. The matrix surface of the model is just wet, which indicates slight leakage rather than significant flow streamlines.Figure 8
Water seepage characteristics at the model boundary.


## 3. Finite Element Analysis of Seepage Evolution in a Filling Joint

### 3.1. Numerical Simulation Schemes

As a result of complex operation, heavy workload, and long time consumed by a large-scale physical model test, a large number of experimental studies with multivariables are difficult to achieve. Thus, here, we applied the finite element analysis software COMSOL multiphysics [1, 4, 31–33] to evaluate the seepage evolution in the filling joint. Here, both the numerical model size and boundary conditions are consistent with the large-scale physical model test (Figure 9). The module of fracture flow was chosen, and the governing equations of fluid flowing through the large-scale rock mass with a filling joint are the well-known Navier-Stokes (NS) equations and mass and momentum conservations, as follows:
(1)∇⋅u=0,ρu⋅∇u‐μ∇2u+∇P=0,where ρ (kg/m3), u (m/s), P (Pa), μ (Pa·s), and t (s) denote the density of a fluid, the velocity vector, the pressure, the viscosity coefficient, and time, respectively.Figure 9
Finite element model using COMSOL multiphysics.The numerical model has a total of 105,261 microelements, and the whole calculation process of fluid flow took approximately 1.2 hours by using a personal computer with an i7 core. To determine the micro seepage parameters of the numerical model, a “trial and error” method was utilized by calibrating the porosity and permeability of both the matrix and the filling joint. All the micro input parameters which were finally confirmed for the numerical model are listed in Table1. Note that due to the much larger size of the model compared to the aperture thickness of the joint, the boundary effects of the three-dimensional model were ignored here in this study. Figure 7(a) presents the comparison between the numerically simulated and experimental variations in pore water pressures, which exhibits a good agreement for the p~t curves. Based on the above micro numerical simulation parameters, the influences of permeability of filling joints (Kj=2×10−10, 4×10−10, 6×10−10, 8×10−10, and 10×10−10m2), permeability of matrix (Km = 2 × 10-12, 4×10−12, 6×10−12, 8×10−12, and 10×10−12m2), joint thickness (b=2, 4, 6, 8, and 10 mm), and inlet water pressure (P=0.1, 0.3, 0.5, 0.7, and 0.9 MPa) on evolution characteristics of the pore water pressure, flow velocity, and flow streamline distribution were investigated, respectively. The variable-controlling approach was adopted, and a total of 17 numerical models were established, as listed in Table 2.Table 1
Microparameters of the numerical model in the simulation.
K
j (10-10 m2)
K
m (10-12 m2)
b (mm)
n
m
n
j
P (MPa)
μ (10−3Pa∙s)
8
2
10
0.20
0.45
0.3
1
Note:Kj = permeability of joint; Km = permeability of matrix; b = joint thickness; nm = porosity of matrix; nj = porosity of joint; P = inlet water pressure; μ = dynamic viscosity of water.Table 2
Numerical simulation conditions in this study.
Model
K
j (10-10 m2)
K
m (10-12 m2)
b (mm)
P (MPa)
Note
#1
2
#2
4
#3
6
2
10
0.3
VariousKj
#4
8
#5
10
#6
2
#7
4
#8
8
6
10
0.3
VariousKm
#9
8
#10
10
#11
2
#12
4
#13
8
2
6
0.3
Variousb
#14
8
#15
10
#16
0.1
#17
0.3
#18
8
2
10
0.5
VariousP
#19
0.7
#20
0.9


### 3.2. Dynamic Evolution Process of the Pore Water Pressure

Variations inp for a rock mass containing a filling joint with respect to various Kj, Km, b, and P values in the t range of 0–3600 s are displayed in Figure 10, in which, the measuring points at L=0.1, 0.3, 0.5, and 0.7 m were chosen, respectively. As t increases, p presents an increasing trend, first indicating a remarkable increase and then reaching a stable value, which can be well characterized using an exponential function. Generally, the increase rate of p declines.Figure 10
Variations inp with respect to various Kj, Km, b, and P.
(a)
K
m
=
2
×
10
−
12
m
2, b=10mm, P=0.3MPa, various Kj
(b)
K
j
=
8
×
10
−
10
m
2, b=10mm, P=0.3MPa, various Km
(c)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, P=0.3MPa, various b
(d)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, b=10mm, various PFigures10(a)–10(c), respectively, show the variations in p in terms of various Kj, Km, and b values. At the initial loading stage of water pressure, a smaller Kj, Km, or b value will result in a smaller p value at the same position in the filling joint. However, when p reaches stable values, the pore water pressure for an identical L is generally the same, which is independent on Kj, Km, or b. For t=3600s, the stable p values for L=0.1, 0.3, 0.5 and 0.7 m are 0.26, 0.17, 0.09, and 0.02 MPa, respectively, indicating a generally linear decrease of 92.31% along the joint length. The above variation trends of p are attributed to the same inlet water pressure P, due to the fact that the pore water pressure in steady states are just affected by the applied inlet water pressure. In order to explore the influences of P on the evolution of pore water pressures, fluid flow simulations on a jointed rock mass with P=0.1, 0.3, 0.5, 0.7, and 0.9 MPa were conducted (Figure 10(d)). For a certain L, as P increases, the stable p exhibits an increase. Taking L=0.3m as an example, in the P range of 0.1–0.9 MPa, the pore water pressure increases from 0.060 to 0.538 MPa, or by a factor of 7.97.


### 3.3. Attenuation Characteristics of the Pore Water Pressure along Joint Length

Similar to the large-scale model test results, for a certaint, the numerically simulated p shows an attenuation along the flowing paths (Figure 11). At the initial stage of fluid flow (e.g., t=100, 300, and 500 s), the pore water pressure shows an obvious nonlinear decrease along the flowing path, indicating a concave shape. However, as t increases, nonlinear characteristics of the p~L curves weaken while the linearity enhances. At t=1000s, the variation in p with L has generally transferred to a linear decrease. From Figure 11(a), nonlinear characteristics of p against L are more obvious for a smaller Kj, with a larger curvature of the p~L curves, which will take longer for the transition of the curves from nonlinearity to linearity. Compared with Kj, the permeability of matrix Km has little influence on p distributions along the joint (Figure 11(b)). From Figure 11(c), the pore water pressure is larger for a larger joint thickness in the early fluid flowing stage, and the nonlinearity of the p~L curves is more obvious for a small b value. In the steady flowing state, the p value is independent on b, indicating a linear decrease with L. Figure 11(d) presents the transition of p distributions for the jointed rock mass with various P. Clearly, at t=1000s, the p~L curve exhibits an increasing reduction rate with P.Figure 11
Numerical results and corresponding fitting curves ofp versus L under different conditions.
(a)
K
m
=
2
×
10
−
12
m
2, b=10mm, P=0.3MPa, various Kj
(b)
K
j
=
8
×
10
−
10
m
2, b=10mm, P=0.3MPa, various Km
(c)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, P=0.3MPa, various b
(d)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, b=10mm, various PGenerally, from Figure11, attenuation characteristics of p with L can be well described using a negative exponential function, as follows:
(2)p=e−AL+B+CeB+C×P,where A, B, and C are fitting coefficients, which are related with Kj, Km, b, and P.Figure12, respectively, shows the variations in coefficient A with Kj, Km, b, and P. As Kj increases, A shows a decrease. Taking t=500s as an example, A decreases from 3.586 to 1.279, by 64.33% in the Kj range of 2×10−10 to 10×10−10m2. Coefficient A first increases then decreases with Km. In the b range of 2–10 mm, A shows a decrease from 3.855 to 1.538, by 60.10% at t=500s. However, coefficient A keeps constant in the P range of 0.1–0.9 MPa.Figure 12
Variations in coefficientA with respect to different Kj, Km, b, and P.
(a)
K
j
~
A
(b)
K
m
~
A
(c)
b
~
A
(d)
P
~
A


### 3.4. Variations in the Flow Velocity

Figure13 presents the changes in flow rate v in the joint with respect to various Kj, Km, b, and P. It can be seen that, from Figures 13(a) and 13(d), after the water pressure was applied, the flow velocity will occupy a larger value for the measure points in the joint with a larger Kj or P. When fluid flowing attains an equilibrium state, with increasing Kj or P, the steady v in the joint shows an increase. Taking L=0.7m in Figure 13(a) as an example, in the Kj range of 2×10−10 to 10×10−10m2, v increases by a factor of 4.42. Additionally, at a small joint length (e.g., L=0.1 and 0.3 m), the flow velocity first increases and then decreases with t, but for a large joint length (L>0.3m), v experiences a monotonous increase with t before it approaches a constant value in the steady flowing state. From Figures 13(b) and 13(c), for L=0.1 and 0.3 m, the flow velocity presents at first an increasing and then a decreasing trend, but for L=0.5 and 0.7 m, v shows an increase until the steady states. In the steady state, the flow velocity is a constant value independent of Km or b.Figure 13
Changes in flow rate through the joint with respect to variousKj, Km, b, and P.
(a)
K
m
=
2
×
10
−
12
m
2, b=10mm, P=0.3MPa, various Kj
(b)
K
j
=
8
×
10
−
10
m
2, b=10mm, P=0.3MPa, various Km
(c)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, P=0.3MPa, various b
(d)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, b=10mm, various PFigure14 shows the water pressure contour planes and flow streamlines through the jointed rock mass with various Kj, Km, b, and P. At the initial fluid flow stage, distribution of pore water pressure presents a “long funnel” shape, which is due to the larger Kj than Km. The diffusion velocity of pore water pressure in the joint is greater than that in the matrix, and fluid flow first reaches the boundary of the model through the joint before it spreads to the whole model end boundary, which is generally consistent with the physical test results in Figure 8. Then, with increasing t, the water pressure contour plane gradually transfers from a “long funnel” to a “short funnel” shape. When the water pressure distribution is steady, the pressure isosurfaces transfer to a series of parallel planes perpendicular to the joint direction.Figure 14
Water pressure contour planes and streamlines of flow in the models with variousKj, Km, b, and P.
(a)
K
m
=
2
×
10
−
12
m
2, b=10mm, P=0.3MPa, Kj=6×10−12m2
(b)
K
j = 8 × 10-10 m2, b=10mm, P=0.3MPa, Km=6×10−12m2
(c)
K
j = 8 × 10-10 m2, Km=2×10−12m2, P=0.3MPa, b=2mm
(d)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, b=10mm, P=0.9MPaFrom the evolution characteristics of flow streamlines through the jointed rock mass, at the initial stage, fluid flowing at the left model boundary generally spreads along the horizontal direction. Later, the flow path changes to the “joint-matrix” direction due to a largerKj compared to Km. Then, as t continuously increases, the flow direction transfers to the horizontal direction along the whole model section, and fluid flowing approaches the steady state.


## 3.1. Numerical Simulation Schemes

As a result of complex operation, heavy workload, and long time consumed by a large-scale physical model test, a large number of experimental studies with multivariables are difficult to achieve. Thus, here, we applied the finite element analysis software COMSOL multiphysics [1, 4, 31–33] to evaluate the seepage evolution in the filling joint. Here, both the numerical model size and boundary conditions are consistent with the large-scale physical model test (Figure 9). The module of fracture flow was chosen, and the governing equations of fluid flowing through the large-scale rock mass with a filling joint are the well-known Navier-Stokes (NS) equations and mass and momentum conservations, as follows:
(1)∇⋅u=0,ρu⋅∇u‐μ∇2u+∇P=0,where ρ (kg/m3), u (m/s), P (Pa), μ (Pa·s), and t (s) denote the density of a fluid, the velocity vector, the pressure, the viscosity coefficient, and time, respectively.Figure 9
Finite element model using COMSOL multiphysics.The numerical model has a total of 105,261 microelements, and the whole calculation process of fluid flow took approximately 1.2 hours by using a personal computer with an i7 core. To determine the micro seepage parameters of the numerical model, a “trial and error” method was utilized by calibrating the porosity and permeability of both the matrix and the filling joint. All the micro input parameters which were finally confirmed for the numerical model are listed in Table1. Note that due to the much larger size of the model compared to the aperture thickness of the joint, the boundary effects of the three-dimensional model were ignored here in this study. Figure 7(a) presents the comparison between the numerically simulated and experimental variations in pore water pressures, which exhibits a good agreement for the p~t curves. Based on the above micro numerical simulation parameters, the influences of permeability of filling joints (Kj=2×10−10, 4×10−10, 6×10−10, 8×10−10, and 10×10−10m2), permeability of matrix (Km = 2 × 10-12, 4×10−12, 6×10−12, 8×10−12, and 10×10−12m2), joint thickness (b=2, 4, 6, 8, and 10 mm), and inlet water pressure (P=0.1, 0.3, 0.5, 0.7, and 0.9 MPa) on evolution characteristics of the pore water pressure, flow velocity, and flow streamline distribution were investigated, respectively. The variable-controlling approach was adopted, and a total of 17 numerical models were established, as listed in Table 2.Table 1
Microparameters of the numerical model in the simulation.
K
j (10-10 m2)
K
m (10-12 m2)
b (mm)
n
m
n
j
P (MPa)
μ (10−3Pa∙s)
8
2
10
0.20
0.45
0.3
1
Note:Kj = permeability of joint; Km = permeability of matrix; b = joint thickness; nm = porosity of matrix; nj = porosity of joint; P = inlet water pressure; μ = dynamic viscosity of water.Table 2
Numerical simulation conditions in this study.
Model
K
j (10-10 m2)
K
m (10-12 m2)
b (mm)
P (MPa)
Note
#1
2
#2
4
#3
6
2
10
0.3
VariousKj
#4
8
#5
10
#6
2
#7
4
#8
8
6
10
0.3
VariousKm
#9
8
#10
10
#11
2
#12
4
#13
8
2
6
0.3
Variousb
#14
8
#15
10
#16
0.1
#17
0.3
#18
8
2
10
0.5
VariousP
#19
0.7
#20
0.9


## 3.2. Dynamic Evolution Process of the Pore Water Pressure

Variations inp for a rock mass containing a filling joint with respect to various Kj, Km, b, and P values in the t range of 0–3600 s are displayed in Figure 10, in which, the measuring points at L=0.1, 0.3, 0.5, and 0.7 m were chosen, respectively. As t increases, p presents an increasing trend, first indicating a remarkable increase and then reaching a stable value, which can be well characterized using an exponential function. Generally, the increase rate of p declines.Figure 10
Variations inp with respect to various Kj, Km, b, and P.
(a)
K
m
=
2
×
10
−
12
m
2, b=10mm, P=0.3MPa, various Kj
(b)
K
j
=
8
×
10
−
10
m
2, b=10mm, P=0.3MPa, various Km
(c)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, P=0.3MPa, various b
(d)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, b=10mm, various PFigures10(a)–10(c), respectively, show the variations in p in terms of various Kj, Km, and b values. At the initial loading stage of water pressure, a smaller Kj, Km, or b value will result in a smaller p value at the same position in the filling joint. However, when p reaches stable values, the pore water pressure for an identical L is generally the same, which is independent on Kj, Km, or b. For t=3600s, the stable p values for L=0.1, 0.3, 0.5 and 0.7 m are 0.26, 0.17, 0.09, and 0.02 MPa, respectively, indicating a generally linear decrease of 92.31% along the joint length. The above variation trends of p are attributed to the same inlet water pressure P, due to the fact that the pore water pressure in steady states are just affected by the applied inlet water pressure. In order to explore the influences of P on the evolution of pore water pressures, fluid flow simulations on a jointed rock mass with P=0.1, 0.3, 0.5, 0.7, and 0.9 MPa were conducted (Figure 10(d)). For a certain L, as P increases, the stable p exhibits an increase. Taking L=0.3m as an example, in the P range of 0.1–0.9 MPa, the pore water pressure increases from 0.060 to 0.538 MPa, or by a factor of 7.97.


## 3.3. Attenuation Characteristics of the Pore Water Pressure along Joint Length

Similar to the large-scale model test results, for a certaint, the numerically simulated p shows an attenuation along the flowing paths (Figure 11). At the initial stage of fluid flow (e.g., t=100, 300, and 500 s), the pore water pressure shows an obvious nonlinear decrease along the flowing path, indicating a concave shape. However, as t increases, nonlinear characteristics of the p~L curves weaken while the linearity enhances. At t=1000s, the variation in p with L has generally transferred to a linear decrease. From Figure 11(a), nonlinear characteristics of p against L are more obvious for a smaller Kj, with a larger curvature of the p~L curves, which will take longer for the transition of the curves from nonlinearity to linearity. Compared with Kj, the permeability of matrix Km has little influence on p distributions along the joint (Figure 11(b)). From Figure 11(c), the pore water pressure is larger for a larger joint thickness in the early fluid flowing stage, and the nonlinearity of the p~L curves is more obvious for a small b value. In the steady flowing state, the p value is independent on b, indicating a linear decrease with L. Figure 11(d) presents the transition of p distributions for the jointed rock mass with various P. Clearly, at t=1000s, the p~L curve exhibits an increasing reduction rate with P.Figure 11
Numerical results and corresponding fitting curves ofp versus L under different conditions.
(a)
K
m
=
2
×
10
−
12
m
2, b=10mm, P=0.3MPa, various Kj
(b)
K
j
=
8
×
10
−
10
m
2, b=10mm, P=0.3MPa, various Km
(c)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, P=0.3MPa, various b
(d)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, b=10mm, various PGenerally, from Figure11, attenuation characteristics of p with L can be well described using a negative exponential function, as follows:
(2)p=e−AL+B+CeB+C×P,where A, B, and C are fitting coefficients, which are related with Kj, Km, b, and P.Figure12, respectively, shows the variations in coefficient A with Kj, Km, b, and P. As Kj increases, A shows a decrease. Taking t=500s as an example, A decreases from 3.586 to 1.279, by 64.33% in the Kj range of 2×10−10 to 10×10−10m2. Coefficient A first increases then decreases with Km. In the b range of 2–10 mm, A shows a decrease from 3.855 to 1.538, by 60.10% at t=500s. However, coefficient A keeps constant in the P range of 0.1–0.9 MPa.Figure 12
Variations in coefficientA with respect to different Kj, Km, b, and P.
(a)
K
j
~
A
(b)
K
m
~
A
(c)
b
~
A
(d)
P
~
A


## 3.4. Variations in the Flow Velocity

Figure13 presents the changes in flow rate v in the joint with respect to various Kj, Km, b, and P. It can be seen that, from Figures 13(a) and 13(d), after the water pressure was applied, the flow velocity will occupy a larger value for the measure points in the joint with a larger Kj or P. When fluid flowing attains an equilibrium state, with increasing Kj or P, the steady v in the joint shows an increase. Taking L=0.7m in Figure 13(a) as an example, in the Kj range of 2×10−10 to 10×10−10m2, v increases by a factor of 4.42. Additionally, at a small joint length (e.g., L=0.1 and 0.3 m), the flow velocity first increases and then decreases with t, but for a large joint length (L>0.3m), v experiences a monotonous increase with t before it approaches a constant value in the steady flowing state. From Figures 13(b) and 13(c), for L=0.1 and 0.3 m, the flow velocity presents at first an increasing and then a decreasing trend, but for L=0.5 and 0.7 m, v shows an increase until the steady states. In the steady state, the flow velocity is a constant value independent of Km or b.Figure 13
Changes in flow rate through the joint with respect to variousKj, Km, b, and P.
(a)
K
m
=
2
×
10
−
12
m
2, b=10mm, P=0.3MPa, various Kj
(b)
K
j
=
8
×
10
−
10
m
2, b=10mm, P=0.3MPa, various Km
(c)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, P=0.3MPa, various b
(d)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, b=10mm, various PFigure14 shows the water pressure contour planes and flow streamlines through the jointed rock mass with various Kj, Km, b, and P. At the initial fluid flow stage, distribution of pore water pressure presents a “long funnel” shape, which is due to the larger Kj than Km. The diffusion velocity of pore water pressure in the joint is greater than that in the matrix, and fluid flow first reaches the boundary of the model through the joint before it spreads to the whole model end boundary, which is generally consistent with the physical test results in Figure 8. Then, with increasing t, the water pressure contour plane gradually transfers from a “long funnel” to a “short funnel” shape. When the water pressure distribution is steady, the pressure isosurfaces transfer to a series of parallel planes perpendicular to the joint direction.Figure 14
Water pressure contour planes and streamlines of flow in the models with variousKj, Km, b, and P.
(a)
K
m
=
2
×
10
−
12
m
2, b=10mm, P=0.3MPa, Kj=6×10−12m2
(b)
K
j = 8 × 10-10 m2, b=10mm, P=0.3MPa, Km=6×10−12m2
(c)
K
j = 8 × 10-10 m2, Km=2×10−12m2, P=0.3MPa, b=2mm
(d)
K
j
=
8
×
10
−
10
m
2, Km=2×10−12m2, b=10mm, P=0.9MPaFrom the evolution characteristics of flow streamlines through the jointed rock mass, at the initial stage, fluid flowing at the left model boundary generally spreads along the horizontal direction. Later, the flow path changes to the “joint-matrix” direction due to a largerKj compared to Km. Then, as t continuously increases, the flow direction transfers to the horizontal direction along the whole model section, and fluid flowing approaches the steady state.


## 4. Theoretical Prediction Model of the Coupling Effects

From the numerical results discussed above, for a rock mass containing a joint, the evolution process of pore water pressurep depends on all the factors of Kj, Km, b, P, L, and t. Variations in p as a function of L can be calculated using equation (2), and variations in p in terms of Kj, Km, b, P, and t can be described using the following equations, respectively:
(3)p=A1+B1×expC1×Kj,(4)p=A2+B2×expC2×Km,(5)p=A3+B3×expC3×b,(6)p=A4+B4×P,(7)p=A5+B5×expC5×t,where A1, A2, A3, A4, A5, B1, B2, B3, B4, B5, C1, C2, C3, and C5 are fitting coefficients.Figure15 presents the regression analysis of p as a function of Kj, Km, b, P, and t, respectively, using the above equations (3), (4), (5), (6), and (7). Here, for Figures 15(a)–15(d), t=600, 1200, 1800, 2400, 3000, and 3600 s were, respectively, chosen, and for Figure 15(e), the p versus t curves were chosen from numerical results for the jointed rock mass with various b values. Obviously, the numerical values show a good agreement with the fitting curves, with all residual squared R2 values larger than 0.95, indicating that equations (3), (4), (5), (6), and (7) fit the relations of p~Kj, p~Km, p~b, p~P, and p~t very well.Figure 15
Regression analysis ofp as a function of Kj, Km, b, P, and t.
(a)
p
~
K
j (L=0.7m, t=600~3600s)
(b)
p
~
K
m (L=0.3m, t=600~3600s)
(c)
p
~
b (L=0.7m, t=600~3600s)
(d)
p
~
P (L=0.5m, t=600~3600s)
(e)
p
~
t (L=0.3m, b=2~10mm)The polynomial interpolation method can be applied to well evaluate the relations between various influencing factors and the pore water pressurep. Here, a seven-dimensional space of p=fKj,Km,b,P,t,L was established by the coupling effects of six factors. According to the effects of these six factors on p described in equations (3), (4), (5), (6), and (7), the following equation (6) can be proposed:
(8)p=ξ0ξ11expξ12t+ξ13+ξ14ξ21expξ22L+ξ23+ξ24ξ31expξ32Kj+ξ33+ξ34ξ41expξ42Km+ξ43+ξ44ξ51expξ52b+ξ53+ξ54ξ61P+ξ62,where ξ0, ξ11, ξ12, ξ13, ξ14, ξ21, ξ22, ξ23, ξ24, ξ31, ξ32, ξ33, ξ34, ξ41, ξ42, ξ43, ξ44, ξ51, ξ52, ξ53, ξ54, ξ61, and ξ62 are decision parameters.Equation (8) is a combination of equations (2), (3), (4), (5), (6), and (7), and whether the calculated p agrees with the numerical results depends on whether the decision parameters are reasonable. Then, a genetic algorithm was constructed to optimize the decision parameters according to the characterization method proposed by Wang and Kong [34] and Wu et al. [24]. This was done through the process of sample series construction, coding, generation of the initial population, calculation of the p values, solution to fitness, selection operation, crossover operation, mutation operation, and the stop breeding condition setting [35, 36]. The indexes including the mean square error (MSE), root mean square error (RMSE), mean absolute error (MAE), and the mean absolute percentage error (MAPE) were introduced to evaluate the correlation between the optimal calculated results and numerical values, as shown in Table 3. Figure 16(a) displays the correlation between the calculated values obtained from the genetic algorithm and numerical results, with the calculated values of the decision parameters in equation (8). The error between the calculated values and numerical results is quite low, with the correlated coefficient R2 value of 0.9946, indicating a high reliability of the genetic algorithm.Table 3
Regression evaluation of the genetic algorithm for the numerical data.
Index
MSE
RMSE
MAE
MAPE (%)
R
R
2
Note
Value
2.60122
×
10
−
6
0.00161
0.00099
0.00024
0.99729
0.99459
p
5.30998
×
10
−
5
0.00729
0.00382
0.01100
0.99929
0.99859
vFigure 16
Correlation between the calculated values of the genetic algorithm and numerical data.
(a)
Pore water pressurep
(b)
Flow velocityvAs discussed in Section3.2, the dynamic evolution process of the flow velocity v is also a function of Kj, Km, b, and P (Figure 17). However, changes in v versus t and L cannot be well characterized using a certain function (Figure 13). For t=3600s, variations in v versus Kj and P can be analyzed with a linear function, while v keeps generally constant with an increase in Km and b. Therefore, a three-dimensional space of v=fKj,P was built to evaluate the coupling effects of Kj and P:
(9)v=η0η1Kj+η2η3P+η4+η5,where η0, η1, η2, η3, η4, and η5 are decision parameters.Figure 17
Regression analysis ofv as a function of Kj, Km, b, P, and t.
(a)
v
~
K
j (t=3600s)
(b)
v
~
K
m (t=3600s)
(c)
v
~
b (t=3600s)
(d)
v
~
P (t=3600s)The indexes including MSE, RMSE, MAE, and MAPE are also listed in Table3, and Figure 16(b) shows the relationships between numerical v and the calculated v values using equation (9), with a R2 value of 0.9986, implying that equation (9) gives a good evaluation of the coupling influences.


## 5. Conclusions

(1)
A conceptual model characterizing the geomechanical states of a water-resistant slab was established, from which, a high-resolution apparatus for seepage tests of a 3D large-scale rock mass and the corresponding experimental method were developed. Then, hydromechanical tests for a large-scale rock mass (300×300×900mm) containing a filling joint with the boundary stress of 5 MPa and an inlet water pressure of 0.3 MPa were conducted
(2)
Evolution characteristics of pore water pressures along the joint length with time were experimentally analyzed, which first sharply increases and then gradually approaches constant values and can be well described using an exponential function. At the steady state, the pore water pressure shows a linear decrease of 92.95% along the joint length. Due to migration of fillings, connected flowing channels are produced and water inrush phenomenon happens
(3)
The influences of filling joint permeability, matrix permeability, joint thickness, and inlet water pressure on seepage evolution were numerically evaluated. The stable pore water pressure, which is independent ofKj, Km or, b, decreases with the joint length but increases with the inlet water pressure. The stable flow velocity presents an increase with Kj and P, but exhibits constant values independent of Km and b. The water pressure contour planes gradually transfer from a “funnel shape” to parallel planes perpendicular to the joint direction, and the flow paths change from “joint-matrix” to horizontal direction along the whole model section
(4)
The comprehensive relation was established to describe the coupling influences ofKj, Km, b, P, L, and t on the pore water pressure and flow velocity, for which, the decision parameters were, respectively, determined and optimized by the genetic algorithm. The discreteness between the calculated results and numerical values was evaluated. The prediction of the evolution characteristics of pore water pressure and flow velocity was realized



---

*Source: 8883861-2021-01-04.xml*