Modelling transient discharge into deep-buried tunnels in karst area based on a coupled discrete-continuum model

Karst aquifers are characteristic of complex hydrogeological conditions and highly developed void structures, often causing serious water and mud inrush, karst collapse and other geological disasters during construction of deep-buried tunnels. Accurate prediction of groundwater discharge into tunnels has become a key issue for tunnel engineering in karst formations. Motivated by safety assessment of tunnel construction in the Jiayan Water Diversion Project located in Guizhou Province, China, we perform groundwater flow simulations to predict transient discharge into the long, deep-buried tunnels. Two numerical models, an equivalent porous medium model and a coupled discrete-continuum model, are used. In the coupled discrete-continuum model, both the tunnels and karst conduits are treated as discrete channels, which are coupled with the surrounding matrix through water exchange. Based on hydrogeological data and site characterization, the model parameters are calibrated by using the PEST inversion algorithm. It is shown that numerical predictions of tunnel discharge during construction agree well with the measurements, indicating reliability of the numerical models. Comparison between the two models shows that the coupled discrete-continuum model can more accurately reflect the changes in the groundwater level and can better estimate the sources of groundwater discharge into the tunnels during the construction. This work suggests that the coupled discrete-continuum model is advantageous over the equivalent porous medium model as the former is a better representation of the karst system. It is demonstrated that the simulation methodology provides a reliable method for assessment of water discharge into tunnels and the optimization of construction scheme. The results are of practical significance for safe construction of tunnels in karst areas.


Introduction
Widely distributed in Southwest China, karst aquifers have strong heterogeneity and highly developed void structures. Groundwater flow processes in karst aquifers are extremely complex, especially under the influence of underground engineering activities. During the construction of deep-buried tunnels in karst areas, excavation would form a new drainage base level, causing profound changes in groundwater runoff and easily inducing sudden water inrush [1]. Common methods for prediction of groundwater inflow into the tunnel can be classified into four categories: empirical methods [2][3], analytical solutions [4][5], numerical solutions based on the equivalent continuum porous media [6][7], and coupled discrete-continuum modelling methods [8][9][10]. In many geological structures, the bulk of groundwater flows along preferred paths or channels within the conduits formed by fracture network [11]. In such situations, discontinuum representation of conduits in coupled discrete-continuum model can give a more realistic representation of the flow system and more accurate prediction of transient discharge into the tunnels.
The coupled discrete-continuum model treats conduits as discrete elements. The flow in the discrete conduits is coupled to the flow in the continuum matrix through usage of water exchange relations. This model representation can accurately reflect the spatial distribution of discrete conduits and simulate the groundwater flow and spring discharge in complex karst aquifers [12]. The objective of this study is to evaluate two numerical modelling approaches, the equivalent porous medium model and the coupled discrete-continuum model, in terms of predicting transient discharge into long, deepburied tunnels in a typical karst area. First, model calibration is performed to obtain uncertain model parameters on the basis of site characterization and water level observations. Then, the transient discharge into the long, deep-buried tunnels is modeled and compared with measurements. Finally, the source of groundwater inflow to the tunnels is estimated during the different phases of construction of the tunnels.

Site description and geology
Located in Bijie and Zunyi in Guizhou Province, China, the Jiayan Water Conservancy Project is a large-scale project with the purpose of transferring water to the areas of high water demand through hydraulic tunnels. The project has a total length of water supply tunnels of about 648 km. In this study, we focus on two sections of the large-scale project: the Shuidaqiao tunnel and the Maochang tunnel (figure 1). The length of the Maochang tunnel and the Shuidaqiao tunnel is 15.696 km and 20.365 km, respectively, with a slope of 1/2800 and 1/3300.
The terrain of the study site is higher in the west and lower in the east, and the altitude is between 1350 and 2200 m. Karst formation is highly developed in this site, and the exposed area of carbonate rocks is more than 60%, mainly in the formations of P2m (limestone) near the Maochang tunnel and of T 1 yn (limestone) near the Shuidaqiao tunnel. The T 1 f siltstone and P 3 l mudstone are of low permeability, and thus divide the study site into two relatively independent parts. The study area has no regional active fault. The Machang fault zone in the south is well cemented and impermeable and other faults in site are permeable. According to the lithology and hydrogeological survey, the hydraulic conductivity of the permeable faults in this study is approximately 5×10 -5 m/s.

Hydrogeological conditions
The annual precipitation in the study site is 1120.8 mm, mainly concentrated in May to October, accounting for more than 80% of the amount in the whole year. Rainfall is the main sources of  [13], the coefficient of the average effective infiltration for non-karst formation is between 0.1 and 0.2 and between 0.35 and 0.5 for karst formations. There are six underground rivers in the study area. The outlet section of the Shuidaqiao tunnel is located between the Yinatianba and Shuchang underground rivers. The excavation of tunnels could encounter karst conduits, which become groundwater leakage outlets, causing sudden water inrush accidents and having serious impact on tunnel construction. The Liuchong river in the south and the Baifu river in the east in study site have relatively stable water levels. The rivers are the lowest drainage base level and become a route for groundwater discharge of the aquifer system.

Coupled discrete-continuum modeling method
In the coupled discrete-continuum approach, the matrix continuum is discretized by using a finite difference model, and coupled flow between the discrete conduit elements and the matrix continuum flow is simulated by the conduit flow process (CFP) model [14]. The coupled discrete-continuum process is simulated using the MODFLOW simulator coupled with the CFPM1 module [15]. Assuming that both rock and water are incompressible, the groundwater flow in continuum matrix is governed by the following equation: where K ij is the hydraulic conductivity tensor (subscripts i, j = 1, 2, 3); h is the hydraulic head; x i (and also x j ) is the coordinate; W is a source/sink term; S S is the specific storage. For each discrete conduit node, the mass balance gives the following equation [15]: where np is the number of conduit segments connected to the node; Q ip is the flux from or to other connected conduits; Q ex is the exchange between the conduit node and the continuum matrix; Q R is the direct recharge to the conduit node; and Q s is the change in storage of the node. The flow Q ip in the conduits can be laminar or turbulent. Different equations are used to describe the laminar flow and the turbulent flow. The flow in the discrete conduits is coupled to the flow in the continuum matrix through water exchange in conduit nodes. The exchange is proportional to head difference between conduit node and continuum matrix [15]: where α i,j,k is the conduit conductance at matrix element i,j,k; h i,j,k is the head at matrix element i,j,k; h n is the head at conduit node n. The conduit conductance α i,j,k depends on the hydraulic conductivity of the conduit wall and conduit geometry. According to [15], the empirical formula is as follows: where K i,j,k is the hydraulic conductivity of matrix element i,j,k; subscript ip means the conduit ip connected to node n; r ip is the radius of conduit ip. Due to space limitation, more details can be found in [15]. (2) the top of the model is recharged by precipitation; (3) the Liuchong River and Baifu River are taken as the constant head boundary, whose prescribed water head is determined by the water level interpolation between the measurement points, (4) ridgeline on the north is set as no-flow boundary, (5) the Machang fault zone is an impermeable fault, which is set as no-flow boundary. The starting and ending nodes of the tunnel are set as water head boundary condition. As a comparison, we also use the equivalent porous medium model to simulate the groundwater flow. For simplicity, in the equivalent porous medium model, the matrix elements containing the tunnel segments are assigned a larger hydraulic conductivity calculated based on the tunnel radius and matrix permeability. A more accurate approach may be to set a potential seepage boundary for the tunnel, which will be tested in the future.

Numerical modeling
The hydraulic conductivity distribution is calibrated by the PEST inversion algorithm [16] based on the field data of groundwater level. Transient flow simulations with monthly precipitation are run to predict the discharge into the Shuidaqiao tunnel. The recharge is calculated as the precipitation multiplied by an effective infiltration coefficient in each parameter zone. The construction schedule of the tunnels is taken into account in the simulations. That is, the length of the tunnels in the model progressively increases according the actual construction schedule (figure 3).

Inverse modeling
Steady state flow models with groundwater level data of 21 observation wells are used to inverse hydraulic conductivities of the karst formations by using the PEST algorithm. The initial groundwater head obtained by the inversed models is essentially in agreement with the observations (figure 4), with root mean square error (RMSE) of 9.26 m for the coupled discrete-continuum model and 15.28 m for the equivalent porous medium model. In coupled discrete-continuum model, the errors between the simulated and the observed groundwater level are all less than 30 m. In comparison, the error for a borehole located in the source of Shuchang underground river is greater than 50 m in equivalent porous medium model. This comparison indicates that the coupled discrete-continuum model can better simulate groundwater flow in karst conduits.
The obtained range of hydraulic conductivities from the inverse modeling is 1×10 -7~2 .26×10 -6 m/s, which is in good agreement with the range of site characterization values estimated from packer tests ( Table 1). The hydraulic conductivity of karst formations is mostly in the order of 10 -7~1 0 -6 m/s in numerous relevant literatures [17][18]. Therefore, the obtained hydraulic conductivities in this paper are reasonable. In both the coupled discrete-continuum model and the equivalent porous medium model, karst formations, such as P 2 m and T 1 yn, have relatively large hydraulic conductivity. It should be noted that the hydraulic conductivity of T 1 yn 2 filled with muddy dolomite is smaller than other karst formations.

Transient tunnel discharge simulation
Based on initial groundwater head and measured precipitation, transient groundwater flow simulations are performed to predict discharge into the Shuidaqiao tunnel. At the end of the construction period, the simulated discharge into Shuidaqiao tunnel is 42682.97 m 3 /d for the coupled discrete-continuum  81 m 3 /d), the numerical result from the coupled discrete-continuum model is more accurate because these empirical formulas usually underestimate the values, especially in heterogeneous geological conditions [9].
The simulated transient tunnel discharge agrees well with measured monthly average tunnel discharge, indicating reliability of the numerical models ( figure 5). In the early stage (from January 2016 to February 2017), the tunnel discharge increases with the length of tunnel. When the length of tunnel is long enough and all sections are almost connected (After July 2017), transient discharge into Shuidaqiao tunnel is closely related to precipitation. In the later stage (from January 2018 to July 2018), the tunnel discharge simulated by equivalent porous medium model tend to be stable, while the tunnel discharge varies with precipitation in coupled discrete-continuum model. After July 2018, the Shuidaqiao tunnel is connected to the Maochang tunnel, resulting in an increase in tunnel discharge. Comparison between the two models suggests that the coupled discrete-continuum model can better respond to precipitation and more accurately simulate transient tunnel discharge during the construction. In order to analyse the source of groundwater discharge into the tunnels, the surface area with groundwater drawdown larger than 10 m is mapped during the construction (figure 6). At the end of tunnel construction, the predicted area of significant drawdown (>10 m) is different, both in shape and quantity, between the coupled discrete-continuum model and the equivalent porous medium model. This is likely due to the different approaches of treating the tunnels in the models. In coupled discretecontinuum model, the tunnel is simulated by discrete line segments embedded into the corresponding matrix grid. In the equivalent porous medium model, the tunnel is simulated by assigning a larger hydraulic conductivity to the grid where the tunnel is located in. In the equivalent porous medium model, the groundwater discharge into the tunnel is concentrated at the outlet of the Shuidaqiao tunnel, which is related to the tunnel construction process. The outlet of the Shuidaqiao tunnel is constructed at the earliest and connected to a long continuous tunnel, becoming a new drainage outlet of groundwater. Consequently, a large part of the area along the tunnel far from the outlet of the Shuidaqiao tunnel has drawdown less than 10 m, so in figure 6b the colored area appears to be discontinuous. Compared with equivalent porous medium model, the source of groundwater discharge into the tunnel in the coupled discrete-continuum model is concentrated on some discrete nodes in the grids at the beginning, then distributed continuously along the tunnel.  Figure 6. Evolution of surface area with significant drawdown during excavation of the tunnel.

Conclusion
Large-scale tunnel construction in karst aquifers often encounters severe water discharge issues which not only threaten the safety of construction but also exert enormous influence on the hydrogeological conditions. In this study, motivated by safety assessment of tunnel construction as part of the Jiayan Water Diversion Project, two different large-scale finite difference models are established based on site conditions and hydrogeological characterization. Groundwater level observations are used to calibrate hydraulic conductivities. Transient simulations are performed to predict the time-series of tunnel discharges. The result shows that coupled discrete-continuum model can better reflect the initial groundwater head than the equivalent porous medium model. The transient simulation results suggest that the groundwater drainage is mainly controlled by tunnel construction process and precipitation.
Comparison between the two models shows that the coupled discrete-continuum model can more accurately reflect the changes in the groundwater level and can better estimate the sources of groundwater discharge into the tunnels during the construction. This research suggests that the coupled discrete-continuum model is advantageous over the equivalent porous medium model in karst aquifers affected by long, deep-buried tunnel. The simulation method provides a reliable way for evaluating tunnel discharge and optimizing construction schemes. The results have practical significance for safe construction of tunnels in karst areas.