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

This study analyzes the spatiotemporal evolution characteristics of seepage through a large-scale rock mass containing a ﬁ lling 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 ﬁ rst 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 ﬂ owing channels induced by migration of ﬁ llings. Using the ﬁ nite element of COMSOL multiphysics, the in ﬂ uences of ﬁ lling 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 ﬂ ow 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 ﬂ owing 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 in ﬂ uences of these factors on the pore water pressure and ﬂ ow 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 ﬁ lling joint.


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][2][3][4][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][7][8]. Additionally, widely distributed disaster-causing structures such as joints, complicated fracture networks, and underground waterbearing structures in surrounding rocks make it more prone to instability of structural plane and groundwater, which would then produce engineering disasters [9][10][11][12][13][14][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 m 3 /s, and 12 points with an extra-large water inflow rate exceeding 1.0 m 3 /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 × 10 5 m 3 /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][17][18][19].
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][21][22][23][24][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 waterresistant 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][27][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.
Thus, in this study, a simplified sketch model describing the geomechanical state of a water-resistant slab was established, as shown in Figure 2. Here, the water-resistant slab is ideally treated as a cuboid model applied with a boundary stress of σ y and σ z in the yand 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 × 900 mm × 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 K j , matrix permeability K m , 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.

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 Figure 3. 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 × 900 mm × 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 yand 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 highpressure 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 m 3 , 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 analysis 2 Geofluids

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 (Figure 4(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  The 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 ( Figure 5(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 ( Figure 5(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-  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 ( Figure 5(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  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 ( Figure 5(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, gasliquid composite water supply system, water inlet device, and embedding arrangement of the water pressure sensors is displayed in Figure 6. 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.
During the test, both boundary stresses in the yand 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-2500 s) 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.
Epoxy resin

Crack-resistant seam belt Polyurethane
Surface cleaning

Sensor connection
Installation of hydraulic sources Figure 5: A typical process of model pouring, sealing, and installation. 6 Geofluids 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 = 1000 s), the variation trend of p versus L gradually transfers from nonlinearity to linearity (Figure 7(b)).
As 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 Figure 8. 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.

Finite Element Analysis of Seepage
Evolution in a Filling Joint  Figure 6: Schematic diagram of the water supply system and embedding location of the water pressure sensors.
where ρ (kg/m 3 ), 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. 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 Table 1. 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 (K j = 2 × 10 −10 , 4 × 10 −10 , 6 × 10 −10 , 8 × 10 −10 , and 10 × 10 −10 m 2 ), permeability of matrix (K m = 2 × 10 -12 , 4 × 10 −12 , 6 × 10 −12 , 8 × 10 −12 , and 10 × 10 −12 m 2 ), 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.

Dynamic Evolution Process of the Pore Water Pressure.
Variations in p for a rock mass containing a filling joint with respect to various K j , K m , 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.
Figures 10(a)-10(c), respectively, show the variations in p in terms of various K j , K m , and b values. At the initial loading stage of water pressure, a smaller K j , K m , 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 K j , K m , or b. For t = 3600 s, 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:3 m 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.  (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 = 1000 s, 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 K j , 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 K j , the permeability of matrix K m 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 = 1000 s, the p~L curve exhibits an increasing reduction rate with P. Generally, from Figure 11, attenuation characteristics of p with L can be well described using a negative exponential function, as follows: where A, B, and C are fitting coefficients, which are related with K j , K m , b, and P. Figure 12, respectively, shows the variations in coefficient A with K j , K m , b, and P. As K j increases, A shows a decrease. Taking t = 500 s as an example, A decreases from 3.586 to 1.279, by 64.33% in the K j range of 2 × 10 −10 to 10 × 10 −10 m 2 . Coefficient A first increases then decreases with K m . In the b range of 2-10 mm, A shows a decrease from 3.855 to 1.538, by 60.10% at t = 500 s. However, coefficient A keeps constant in the P range of 0.1-0.9 MPa.       Figure 13 presents the changes in flow rate v in the joint with respect to various K j , K m , 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 K j or P. When fluid flowing attains an equilibrium state, with increasing K j or P, the steady v in the joint shows an increase. Taking L = 0:7 m in Figure 13(a) as an example, in the K j range of 2 × 10 −10 to 10 × 10 −10 m 2 , 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:3 m), 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 K m or b. Figure 14 shows the water pressure contour planes and flow streamlines through the jointed rock mass with various K j , K m , 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 K j than K m . 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  Figure 13: Changes in flow rate through the joint with respect to various K j , K m , b, and P. 16 Geofluids transfer to a series of parallel planes perpendicular to the joint direction. From 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 larger K j compared to K m . 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.

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 pressure p depends on all the factors of K j , K m , 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 K j , K m , b, P, and t can be described using the following equations, respectively: where , and C 5 are fitting coefficients. Figure 15 presents the regression analysis of p as a function of K j , K m , 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   good agreement with the fitting curves, with all residual squared R 2 values larger than 0.95, indicating that equations (3), (4), (5), (6), and (7) fit the relations of p~K j , p~K m , p~b, p~P, and p~t very well. The polynomial interpolation method can be applied to well evaluate the relations between various influencing factors and the pore water pressure p. Here, a sevendimensional space of p = f ðK j , K m , 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: where ξ 0 , ξ 11 , ξ 12 , ξ 13 , ξ 14 (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 R 2 value of 0.9946, indicating a high reliability of the genetic algorithm.
As discussed in Section 3.2, the dynamic evolution process of the flow velocity v is also a function of K j , K m , 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 = 3600 s, variations in v versus K j and P can be analyzed with a linear function, while v keeps generally constant with an increase in K m and b. Therefore, a three-dimensional space of v = f ðK j , PÞ was built to evaluate the coupling effects of K j and P: where η 0 , η 1 , η 2 , η 3 , η 4 , and η 5 are decision parameters.
The indexes including MSE, RMSE, MAE, and MAPE are also listed in Table 3, and Figure 16(b) shows the relationships between numerical v and the calculated v values using equation (9), with a R 2 value of 0.9986, implying that equation (9) gives a good evaluation of the coupling influences.

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 × 900 mm) 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 of K j , K m or, b, decreases with the joint length but increases with the inlet water pressure. The stable flow velocity presents an increase with K j and P, but exhibits constant values independent of K m 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 of K j , K m , 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

Data Availability
If needed, we can supply all the data.

Conflicts of Interest
The authors declare that they have no conflicts of interest.