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.
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 , 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].
(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.
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 and in the - and -directions in the waterproof strata, a water pressure 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 () 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 , matrix permeability , joint thickness , and the inlet hydraulic pressure . 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 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 (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 - and -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 analysis
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 (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 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 to .
(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 tests
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–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))
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 (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 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 (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, gas-liquid 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 - and -directions were constant at 5 MPa, and the inlet water pressure was set to be 0.3 MPa. Then, variations in the pore water pressure with time () and the joint length (, 220, 380, 540, 700 mm) can be analyzed. Figure 7(a) presents the spatiotemporal evolution characteristics of . It can be seen that the pore water pressure at all the five positions through the filling joint shows an increase with . At the beginning, changes significantly with . However, as continuously increases, the increasing rate steadily diminishes and gradually approaches a constant value. The variation of as a function of can be well described using an exponential function. In addition, the variation rate of at various positions along the joint is different, indicating a larger value for the position close to the water inlet. For , the constant 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 (e.g., , 300 s), the decrease in with presents a nonlinear trend. With increasing (), the variation trend of versus gradually transfers from nonlinearity to linearity (Figure 7(b)).
(a)
(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.
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: where (kg/m3), (m/s), (Pa), (), and (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 curves. Based on the above micro numerical simulation parameters, the influences of permeability of filling joints (, , , , and ), permeability of matrix (Km = 2 × 10-12, , , , and ), joint thickness (, 4, 6, 8, and 10 mm), and inlet water pressure (, 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.
3.2. Dynamic Evolution Process of the Pore Water Pressure
Variations in for a rock mass containing a filling joint with respect to various , , , and values in the range of 0–3600 s are displayed in Figure 10, in which, the measuring points at , 0.3, 0.5, and 0.7 m were chosen, respectively. As increases, 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 declines.
(a) , , , various
(b) , , , various
(c) , , , various
(d) , , , various
Figures 10(a)–10(c), respectively, show the variations in in terms of various , , and values. At the initial loading stage of water pressure, a smaller , , or value will result in a smaller value at the same position in the filling joint. However, when reaches stable values, the pore water pressure for an identical is generally the same, which is independent on , , or . For , the stable values for , 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 are attributed to the same inlet water pressure , 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 on the evolution of pore water pressures, fluid flow simulations on a jointed rock mass with , 0.3, 0.5, 0.7, and 0.9 MPa were conducted (Figure 10(d)). For a certain , as increases, the stable exhibits an increase. Taking as an example, in the 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 certain , the numerically simulated shows an attenuation along the flowing paths (Figure 11). At the initial stage of fluid flow (e.g., , 300, and 500 s), the pore water pressure shows an obvious nonlinear decrease along the flowing path, indicating a concave shape. However, as increases, nonlinear characteristics of the curves weaken while the linearity enhances. At , the variation in with has generally transferred to a linear decrease. From Figure 11(a), nonlinear characteristics of against are more obvious for a smaller , with a larger curvature of the curves, which will take longer for the transition of the curves from nonlinearity to linearity. Compared with , the permeability of matrix has little influence on 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 curves is more obvious for a small value. In the steady flowing state, the value is independent on , indicating a linear decrease with . Figure 11(d) presents the transition of distributions for the jointed rock mass with various . Clearly, at , the curve exhibits an increasing reduction rate with .
(a) , , , various
(b) , , , various
(c) , , , various
(d) , , , various
Generally, from Figure 11, attenuation characteristics of with can be well described using a negative exponential function, as follows: where , , and are fitting coefficients, which are related with , , , and .
Figure 12, respectively, shows the variations in coefficient with , , , and . As increases, shows a decrease. Taking as an example, decreases from 3.586 to 1.279, by 64.33% in the range of to . Coefficient first increases then decreases with . In the range of 2–10 mm, shows a decrease from 3.855 to 1.538, by 60.10% at . However, coefficient keeps constant in the range of 0.1–0.9 MPa.
(a)
(b)
(c)
(d)
3.4. Variations in the Flow Velocity
Figure 13 presents the changes in flow rate in the joint with respect to various , , , and . 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 or . When fluid flowing attains an equilibrium state, with increasing or , the steady in the joint shows an increase. Taking in Figure 13(a) as an example, in the range of to , increases by a factor of 4.42. Additionally, at a small joint length (e.g., and 0.3 m), the flow velocity first increases and then decreases with , but for a large joint length (), experiences a monotonous increase with before it approaches a constant value in the steady flowing state. From Figures 13(b) and 13(c), for and 0.3 m, the flow velocity presents at first an increasing and then a decreasing trend, but for and 0.7 m, shows an increase until the steady states. In the steady state, the flow velocity is a constant value independent of or .
(a) , , , various
(b) , , , various
(c) , , , various
(d) , , , various
Figure 14 shows the water pressure contour planes and flow streamlines through the jointed rock mass with various , , , and . At the initial fluid flow stage, distribution of pore water pressure presents a “long funnel” shape, which is due to the larger than . 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 , 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.
(a) , , ,
(b) Kj = 8 × 10-10 m2, , ,
(c) Kj = 8 × 10-10 m2, , ,
(d) , , ,
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 compared to . Then, as 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 pressure depends on all the factors of , , , , , and . Variations in as a function of can be calculated using equation (2), and variations in in terms of , , , , and can be described using the following equations, respectively: where , , , , , , , , , , , , , and are fitting coefficients.
Figure 15 presents the regression analysis of as a function of , , , , and , respectively, using the above equations (3), (4), (5), (6), and (7). Here, for Figures 15(a)–15(d), , 1200, 1800, 2400, 3000, and 3600 s were, respectively, chosen, and for Figure 15(e), the versus curves were chosen from numerical results for the jointed rock mass with various values. Obviously, the numerical values show a good agreement with the fitting curves, with all residual squared values larger than 0.95, indicating that equations (3), (4), (5), (6), and (7) fit the relations of , , , , and very well.
(a) (, )
(b) (, )
(c) (, )
(d) (, )
(e) (, )
The polynomial interpolation method can be applied to well evaluate the relations between various influencing factors and the pore water pressure . Here, a seven-dimensional space of was established by the coupling effects of six factors. According to the effects of these six factors on described in equations (3), (4), (5), (6), and (7), the following equation (6) can be proposed: where , , , , , , , , , , , , , , , , , , , , , , and are decision parameters.
Equation (8) is a combination of equations (2), (3), (4), (5), (6), and (7), and whether the calculated 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 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 value of 0.9946, indicating a high reliability of the genetic algorithm.
(a) Pore water pressure
(b) Flow velocity
As discussed in Section 3.2, the dynamic evolution process of the flow velocity is also a function of , , , and (Figure 17). However, changes in versus and cannot be well characterized using a certain function (Figure 13). For , variations in versus and can be analyzed with a linear function, while keeps generally constant with an increase in and . Therefore, a three-dimensional space of was built to evaluate the coupling effects of and : where , , , , , and are decision parameters.
(a) ()
(b) ()
(c) ()
(d) ()
The indexes including MSE, RMSE, MAE, and MAPE are also listed in Table 3, and Figure 16(b) shows the relationships between numerical and the calculated values using equation (9), with a 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 () 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 , or, , decreases with the joint length but increases with the inlet water pressure. The stable flow velocity presents an increase with and , but exhibits constant values independent of and . 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 , , , , , and 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.
Acknowledgments
The authors gratefully acknowledge the financial support from the National Natural Science Foundation of China (51904290, 51734009), the Natural Science Foundation of Jiangsu Province, China (BK20180663), and the Opening Fund of State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (Chengdu University of Technology), China (SKLGP2020K021).