Coupled effects of fracture spatial distribution and rock matrix permeability on tunnel water inflow: An embedded discrete fracture model based investigation

Abstract

Groundwater inflow constitutes a critical challenge in rock tunnel engineering. This study systematically investigates the coupled effects of fracture spatial distribution and rock matrix permeability on tunnel water inflow using a novel embedded discrete fracture model based method. A set of quadratic regression models is established to delineate the relationship between inflow rate and fracture distribution parameters over a wide range of fracture-to-matrix permeability ratios (k_{\mathrm{f}}/k_{\mathrm{m}}). Results demonstrate that fracture aperture, spacing, and their interaction dominate the inflow across all permeability ratios. Analysis of variance further reveals a threshold-dependent behavior: coupled effects are significant below a critical k_{\mathrm{f}}/k_{\mathrm{m}} value but decay markedly above it. This threshold decreases with larger aperture and increases with wider spacing, yet remains nearly independent of fracture dip angle. Moreover, when k_{\mathrm{f}}/k_{\mathrm{m}} is below the threshold, aperture and spacing exert greater influence on tunnel inflow at lower permeability ratios, while k_{\mathrm{f}}/k_{\mathrm{m}} gains influence under larger apertures and smaller spacings. Finally, a case study of Nanwan Tunnel shows that matrix permeability plays a dual role—increasing the mean inflow rate while reducing uncertainty from stochastic fracture distribution.

Graphical abstract

Keywords

tunnel water inflow / fractured rock tunnel / fracture spatial distribution / rock matrix permeability / embedded discrete fracture model / fracture-to-matrix permeability ratio

1 Introduction

Groundwater inflow constitutes a principal engineering challenge in tunnelling through fractured rock masses [1]. Statistical analysis of major tunnel construction accidents in China (2010–2020) indicates that groundwater inflow-related hazards caused 38 fatalities, ranking as the third most severe hazard following collapse and explosion [2]. Groundwater inflow compromises tunnel stability through increased hydraulic pressure on linings [3,4], while simultaneously inducing environmental ramifications such as groundwater table drawdown [5], vegetation deterioration [6], and ground settlement [7]. Accurate estimation of tunnel inflow rate is therefore critical for risk mitigation [8]. Overestimation leads to overdesigned waterproofing and drainage systems [9], whereas underestimation elevates disaster risks [10]. Reliable estimation remains challenging due to the complex interplay of controlling factors such as rock mass permeability, groundwater level, excavation geometry, and boundary conditions [11].

Fractures serve as the primary flow conduits from surrounding rock masses to tunnels, with their spatial distribution significantly influencing inflow rate [12]. Previous research has extensively examined the influence of uniformly distributed fractures on tunnel inflow, where fracture spatial distribution is parameterized by fracture spacing, orientation, and hydraulic aperture [11]. These studies demonstrate that under such conditions, tunnel inflow rate is predominantly governed by hydraulic aperture and fracture spacing, with aperture exerting the most pronounced influence [13–15]. In practical situations, fracture spatial distribution commonly exhibits substantial heterogeneity and anisotropy [16], markedly diverging from the uniform distribution. Consequently, the overall distribution characteristics of fracture networks, beyond just fracture aperture, also significantly influence tunnel water inflow. For instance, Cesano et al. [16] reported that increased fracture network heterogeneity leads to greater variability in inflow magnitude. Bang et al. [17] revealed that fracture length-aperture correlation appears to be effective in influencing tunnel inflow.

Most above-mentioned studies assume an impermeable rock matrix, considering only fracture flow [18]. While valid when fracture permeability exceeds matrix permeability by several orders of magnitude, this assumption becomes inadequate under specific geological conditions [19], including: 1) substantial fracture aperture reduction from geochemical (e.g., mineral precipitation) or geomechanical effects (e.g., compaction) [19,20]; 2) sparsely fractured rock masses with non-negligible matrix permeability [21]; 3) lithologies with significant matrix porosity such as porous carbonates [22]. Routinely neglecting matrix permeability without rigorous justification may compromise the generalizability of existing findings, since even moderate matrix permeability can alter fracture flow and, consequently, tunnel inflow. Limited studies have addressed this knowledge gap. Huang et al. [23] explicitly modeled artery fractures while treating ramification fractures and rock matrix as equivalent porous media [24], identifying fracture aperture as the dominant control on tunnel inflow. Similarly, Zhang et al. [25] developed a dual-permeability model for stratified jointed rock masses, employing discrete fracture techniques for major stratification planes and treating secondary fractures and rock matrix as porous media. Their results indicate that the spacing of primary fractures substantially affects water inflow, whereas fracture orientation exhibits minimal impact. Both studies suffer from limited applicability due to their reliance on predefined fracture hierarchies. More importantly, there is a lack of systematic investigation into the coupled effects of fracture spatial distribution and matrix permeability on tunnel water inflow. Specifically, the influence of fracture spatial distribution across varying fracture-to-matrix permeability ratios remains unclear, and the dominant controls on tunnel inflow require further elucidation.

A major challenge in studying the coupled effects lies in establishing a tunnel inflow model that integrates fracture and matrix flow. As mentioned above, the models proposed by Huang et al. [23] and Zhang et al. [25] are limited to stratified jointed rock masses with dominant fracture sets and cannot be readily extended to non-stratified fracture network systems. An alternative model involves explicitly modeling all fractures and utilizing conforming mesh techniques to simulate fracture-matrix hydraulic interactions [26], as implemented in commercial software such as COMSOL [27] and 3DEC [28]. Hokr et al. [29] successfully applied this model to quantify water inflow into a tunnel in northern Czechia. However, this model is limited by its computationally intensive mesh generation process for complex fracture networks, which can become prohibitively expensive or even numerically unstable. Some researchers have utilized the finite-discrete element method (FDEM) to simulate fluid flow in rock masses [30,31]. However, FDEM is generally more suitable for hydro-mechanical coupling problems. For pure flow simulation, FDEM tends to be computationally intensive and bulky, as it explicitly resolves complex fracture initiation, propagation, and contact interactions. The unified pipe network method (UPM) [32–34] simplifies flow simulation by treating the fractured rock mass as a pipe network system. Nevertheless, this strategy may compromise physical fidelity, especially in three-dimensional settings with complex fracture geometries, due to the simplified representation of the flow domain. The embedded discrete fracture model (EDFM) is a recently developed methodology for simulating fluid flow in fractured porous media [35,36]. A major advantage of EDFM over traditional conforming mesh methods lies in its ability to decouple the discretization of fracture and matrix. In this approach, fractures are embedded as lower-dimensional entities into a simple background matrix mesh. Hydraulic interaction between fracture and matrix is achieved through additional flow exchange terms rather than mesh conformity. This non-conforming strategy circumvents the geometric complexity caused by mesh conformity, reduces numerical instability, and improves computational efficiency. On the other hand, EDFM may offer higher physical fidelity than UPM under some complex conditions because it resolves fracture and matrix explicitly through meshes instead of pipes. These attributes make EDFM particularly well-suited for efficient flow simulation in complex fractured porous media. EDFM was originally developed and popularized within the reservoir simulation community for simulating oil and gas flow [37–39]. By contrast, the groundwater modeling community has historically favored methods such as COMSOL and UPM. As a result, the potential of EDFM for groundwater systems remains largely unexplored. This study employs an EDFM-based tunnel inflow model, formulated preliminarily in our prior work [40] and detailed in the present paper, to numerically investigate the coupled effects of fracture spatial distribution and matrix permeability on tunnel water inflow.

In summary, while existing theories adequately capture the influence of fracture spatial distribution on tunnel water inflow under impermeable matrix conditions, the coupled effects of fracture spatial distribution and matrix permeability are not well understood. To address this gap, we develop an EDFM-based tunnel inflow model incorporating both fracture spatial distribution and matrix permeability to obtain tunnel inflow under diverse geological conditions. Using the simulation results, we establish a set of quadratic regression models to delineate the effects of fracture spatial distribution across a wide range of fracture-to-matrix permeability ratios. Analysis of variance (ANOVA) [41] is employed to quantify the significance of the coupled effects. The Nanwan Tunnel project is selected as a case to evaluate the coupled effects of fracture network properties and matrix permeability on the total tunnel inflow rate.

2 Numerical methodology

2.1 Physical model

2.1.1 Tunnel model

Figure 1(a) illustrates the spatial configuration of the tunnel model. To evaluate the impact of fracture spatial distribution and matrix permeability exclusively, the following simplifications are adopted [11]: 1) the tunnel cross-section is simplified as circular; 2) waterproofing and drainage systems are neglected; 3) the computational domain is defined as a square.

Figure 1(b) presents the geometry and boundary conditions of the computational domain. As reported by Farhadian et al. [11], an insufficient model extent (ME) results in an overestimation of the inflow rate, while an excessively large ME incurs substantial and unnecessary computational costs. In this study, the optimum ME is determined based on the empirical chart proposed by Farhadian et al. [11]. The boundary conditions are defined as follows: 1) the top boundary is subjected to a constant water pressure [42]; 2) the left and right boundaries are assigned a linearly varying water pressure [11]; 3) the bottom boundary is modeled either as impermeable (no-flow) or under a constant water pressure [42,43]; 4) the tunnel boundary is set to zero water pressure [42].

2.1.2 Fractured rock mass model

As introduced in the Introduction, this study employs EDFM to characterize fractured rock masses. The EDFM adopts a hybrid-dimensional modeling approach, in which fractures are represented as (n–1)-dimensional entities embedded within an n-dimensional porous matrix [26]. As illustrated conceptually in Fig. 2, the EDFM integrates two core representations: 1) discrete fracture networks (DFNs) for explicit fracture characterization; 2) equivalent porous media for matrix representation. Hydraulic coupling between the fracture and matrix domains is achieved through algebraically formulated flow exchange terms (detailed in subsequent sections).

2.1.3 Groundwater flow model

The steady-state Darcy flow model is employed to simulate groundwater flow in fractured rock masses [44]. Within a two-dimensional domain containing n_{\mathrm{f}} fractures, the mass conservation equations for the matrix and fracture domains are given as follows.

1) Matrix domain ({\mathrm{\Omega }}_{\mathrm{m}}\subset R^2):

(1)\mathrm{\nabla }\cdot v_{\mathrm{m}}=-\sum _{i=1}^{n_{\mathrm{f}}}{q}_{\mathrm{m}i},

where v_{\mathrm{m}} (m/s) is the groundwater velocity in the rock matrix, and q_{\mathrm{m}i} (m²/s/m²) is the groundwater flow rate from the matrix to fracture i per unit matrix area.

2) Fracture domain ({\mathrm{\Omega }}_{{\mathrm{f}}_i}\subset R^1, for i\in \left[1,n_{\mathrm{f}}\right]):

(2)\mathrm{\nabla }\cdot v_i=\left(-q_{i\mathrm{m}}-\sum _{j=1,j\ne i}^{n_{\mathrm{f}}}{q}_{ij}\right)/a_i,

where v_i (m/s) is the groundwater velocity in fracture i, q_{i\mathrm{m}} (m²/s/m) is the groundwater flow rate from fracture i to the matrix per unit fracture length, q_{ij} (m²/s/m) is the groundwater flow rate from fracture i to fracture j per unit length, and a_i (m) is the aperture of fracture i.

According to Darcy’s law [45], groundwater velocity v (m/s) is expressed as:

(3)v=-\frac{k}{\mu }\left(\mathrm{\nabla }p-\rho \mathrm{g}\right),

where k (m²) is the permeability of the medium, \mu =1\times {10}^{-3} Pa·s is the dynamic viscosity of water, p (Pa) is the water pressure, \rho =1000 kg/m³ is the water density, and \mathrm{g}=\left[0,-9.8\right] m/s² is the gravitational acceleration vector.

Fracture permeability is given by the cubic law [46]:

(4)k_{\mathrm{f}}=\frac{a^2}{12},

where k_{\mathrm{f}} (m²) is the fracture permeability, and a (m) is the fracture aperture.

Substituting Eq. (3) into Eqs. (1) and (2) yields the final governing equations:

1) Matrix domain ({\mathrm{\Omega }}_{\mathrm{m}}\subset R^2):

(5)\mathrm{\nabla }\cdot \left[-\frac{k_{\mathrm{m}}}{\mu }\left(\mathrm{\nabla }{p}_{\mathrm{m}}-\rho \mathrm{g}\right)\right]=-\sum _{i=1}^{n_{\mathrm{f}}}{q}_{\mathrm{m}i},

where k_{\mathrm{m}} (m²) is the matrix permeability, and p_{\mathrm{m}} (Pa) is the water pressure in the matrix.

2) Fracture domain ({\mathrm{\Omega }}_{{\mathrm{f}}_i}\subset R^1, for i\in \left[1,n_{\mathrm{f}}\right]):

(6)\mathrm{\nabla }\cdot \left[-\frac{k_i}{\mu }\left(\mathrm{\nabla }{p}_i-\rho \mathrm{g}\right)\right]=\left(-q_{i\mathrm{m}}-\sum _{j=1,j\ne i}^{n_{\mathrm{f}}}{q}_{ij}\right)/a_i,

where k_i (m²) is the permeability of fracture i, and p_i (Pa) is the water pressure in fracture i.

2.2 Numerical method

2.2.1 Mesh generation

The EDFM decouples fracture and matrix discretization via a non-conforming mesh strategy. As illustrated in Fig. 3, we implement the strategy comprising three stages: 1) matrix discretization using a structured Cartesian grid [47] as the background mesh; 2) fracture embedding through geometric superposition of discrete fractures onto the matrix mesh; 3) fracture mesh generation by isolating fracture-matrix intersection zones as independent computational elements.

2.2.2 Two-point flux approximation (TPFA)

TPFA method [48] is used for spatial discretization. In TPFA, cells (i.e., meshes) are interconnected through a one-to-one relationship based on the transmissibility concept [26]. As shown in Fig. 4(a), the flux q_{ij} between two adjacent cells {\mathrm{\Omega }}_i and {\mathrm{\Omega }}_j is approximated as:

(7)q_{ij}\approx T_{ij}\left(u_i-u_j\right),

where q_{ij} (m²/s) is the flux across interface {\mathrm{\Gamma }}_{ij}, T_{ij} (m/s) is the transmissibility of {\mathrm{\Gamma }}_{ij}, and u_i (m) and u_j (m) are the potential values at the cell centers of {\mathrm{\Omega }}_i and {\mathrm{\Omega }}_j, respectively.

Let \overline{u} (m) be the potential at {\mathrm{\Gamma }}_{ij} and assume a linear pressure distribution within each cell. By applying Darcy’s law and enforcing flux continuity across the interface, the following relationship is obtained:

(8)-2\text{Δy}{K}_i\left(\overline{u}-u_i\right)/{\mathrm{\Delta }x}_i=2\mathrm{\Delta }y{K}_j\left(\overline{u}-u_j\right)/{\mathrm{\Delta }x}_j,

where K_i (m/s) and K_j (m/s) represent the hydraulic conductivities of {\mathrm{\Omega }}_i and {\mathrm{\Omega }}_j, respectively.

Define T_i (m/s) and T_j (m/s) as:

(9a)\begin{array}{c}{T}_i=2K_i\mathrm{\Delta }\mathrm{y}/{\mathrm{\Delta }x}_i,\end{array}

(9b)\begin{array}{c}{T}_j=2K_j\mathrm{\Delta }\mathrm{y}/{\mathrm{\Delta }x}_j.\end{array}

\overline{u} can be derived from Eq. (8):

(10)\overline{u}=\left(T_i{u}_i+T_j{u}_j\right)/\left(T_i+T_j\right).

Substituting Eq. (10) into the left-hand side of Eq. (8) yields the interfacial flux q_{ij}:

(11)q_{ij}=\left(u_i-u_j\right)/\left(1/T_i+1/T_j\right).

Comparing Eq. (11) with Eq. (7), the interfacial transmissibility T_{ij} is obtained:

(12)T_{ij}=1/\left(1/T_i+1/T_j\right).

In summary, Eq. (7) describes the general TPFA numerical scheme, with computation of interfacial transmissibility being its key component. This transmissibility depends on the hydraulic conductivity and geometry of adjacent cells.

2.2.3 Discretization of the physical model

To discretize the tunnel inflow model using TPFA, the left-hand sides of Eqs. (5) and (6) must be reformulated as interfacial fluxes. Applying the divergence theorem to Eqs. (5) and (6) gives the integral form.

1) Matrix cell ({\mathrm{\Omega }}_i\subset R^2, for i\in \left[1,N_{\mathrm{m}}\right])

(13){\int }_{\mathrm{\partial }{\mathrm{\Omega }}_i}\left[-\frac{k_i}{\mu }\left(\mathrm{\nabla }{p}_i-\rho \mathrm{g}\right)\right]\cdot n_i\mathrm{d}s=-\sum _{j=1}^{N_{\mathrm{f}}}{\int }_{{\mathrm{\Omega }}_i}{q}_{ij}\mathrm{d}x,

where N_{\mathrm{m}} is the number of matrix cells, N_{\mathrm{f}} is the number of fracture cells, {\mathrm{\Omega }}_i is the i\text{th} matrix cell, k_i (m²) is the permeability of {\mathrm{\Omega }}_i, p_i (Pa) is the water pressure in {\mathrm{\Omega }}_i, n_i is the unit outward normal vector of the boundary of {\mathrm{\Omega }}_i, and q_{ij} (m²/s/m²) is the groundwater flow rate from {\mathrm{\Omega }}_i to fracture cell {\mathrm{\Omega }}_j per unit area of {\mathrm{\Omega }}_i.

The left-hand side of Eq. (13) is the matrix flux term (m²/s), and {\int }_{{\mathrm{\Omega }}_i}{q}_{ij}\mathrm{d}x is the matrix-fracture flow exchange term (m²/s) between {\mathrm{\Omega }}_i and {\mathrm{\Omega }}_j (nonzero only when {\mathrm{\Omega }}_i intersects with {\mathrm{\Omega }}_j).

2) Fracture cell ({\mathrm{\Omega }}_j\subset R^1, for i\in \left[1,N_{\mathrm{f}}\right])

(14)\begin{array}{rl}& {\int }_{\mathrm{\partial }{\mathrm{\Omega }}_j}\left[-\frac{k_j}{\mu }\left(\mathrm{\nabla }{p}_j-\rho \mathrm{g}\right)\right]\cdot n_j\mathrm{d}s\\ =& -\sum _{i=1}^{N_{\mathrm{m}}}{\int }_{{\mathrm{\Omega }}_j}{q}_{ji}\mathrm{d}x-\sum _{k=1,k\ne j}^{N_{\mathrm{f}}}{\int }_{{\mathrm{\Omega }}_j}{q}_{jk}\mathrm{d}x,\end{array}

where {\mathrm{\Omega }}_j is the j\text{th} fracture cell, k_j (m²) is the permeability of {\mathrm{\Omega }}_j, p_j (Pa) is the water pressure in {\mathrm{\Omega }}_j, n_j is the unit outward normal vector of the boundary between {\mathrm{\Omega }}_j and neighboring fracture cells, q_{ji} (m²/s/m) is the groundwater flow rate from {\mathrm{\Omega }}_j to matrix cell {\mathrm{\Omega }}_i per unit length of {\mathrm{\Omega }}_j, and q_{jk} (m²/s/m) is the groundwater flow rate exchanged from {\mathrm{\Omega }}_j to fracture cell {\mathrm{\Omega }}_k per unit length of {\mathrm{\Omega }}_j.

The left-hand side of Eq. (14) is the fracture flux term, {\int }_{{\mathrm{\Omega }}_j}{q}_{ji}\mathrm{d}x (m²/s) is the fracture-matrix flow exchange term (equal to the negative of the corresponding matrix-fracture flow exchange term), and {\int }_{{\mathrm{\Omega }}_j}{q}_{jk}\mathrm{d}x (m²/s) is the fracture-fracture flow exchange term between intersecting fractures (nonzero only when {\mathrm{\Omega }}_j intersects with {\mathrm{\Omega }}_k).

Figures 4(b)–4(e) depicts the TPFA discretization of the terms in Eqs. (13) and (14). Pressure is defined at the centroid of each matrix and fracture cell. Due to the non-conforming nature of matrix and fracture meshes, pressure continuity is not enforced geometrically through shared nodes. Instead, it is achieved through algebraic connectivity between cells, which is established by computing the transmissibility-based flow exchange for matrix-matrix, matrix-fracture, and fracture-fracture interactions. The equilibrium of these flow exchanges naturally approximates a continuous pressure field, thereby bypassing the need for mesh conformity. The discrete forms of the terms in Eqs. (13) and (14) are derived as follows.

1) Matrix flux term

The matrix flux term (Fig. 4(b)) is discretized as:

(15a)\begin{array}{rl}& {\int }_{\mathrm{\partial }{\mathrm{\Omega }}_i}\left[-\frac{k_i}{\mu }\left(\mathrm{\nabla }{p}_i-\rho \mathrm{g}\right)\right]\cdot n_i\mathrm{d}s\approx \sum _{j=1}^k{T}_{ij}\left(u_i-u_j\right),\\ & \text{for}i\in \left[1,N_{\mathrm{m}}\right],{\mathrm{\Omega }}_j\in \text{neighbor}\left({\mathrm{\Omega }}_i\right),\end{array}

(15b)\begin{array}{c}{T}_{ij}=1/\left(1/T_i+1/T_j\right),\end{array}

(15c)\begin{array}{c}{T}_i=K_i{l}_{ij}/d_i,\end{array}

(15d)\begin{array}{c}{T}_j=K_j{l}_{ij}/d_j,\end{array}

where {\mathrm{\Omega }}_j is the adjacent cell sharing an interface with {\mathrm{\Omega }}_i, k is the number of neighboring cells of {\mathrm{\Omega }}_i, u_j and u_j are the potential values at the centers of {\mathrm{\Omega }}_i and {\mathrm{\Omega }}_j, respectively, K_i and K_j represent the hydraulic conductivities of {\mathrm{\Omega }}_i and {\mathrm{\Omega }}_j, respectively.

2) Matrix-fracture flow exchange term

The matrix-fracture flow exchange term (Fig. 4(c)) is discretized as:

(16a)\begin{array}{c}{\int }_{{\mathrm{\Omega }}_i}{q}_{ij}\mathrm{d}x\approx T_{ij}\left(u_i-u_j\right),\text{for}i\in \left[1,N_{\mathrm{m}}\right],{\mathrm{\Omega }}_i\cap {\mathrm{\Omega }}_j\ne \varnothing ,\end{array}

(16b)\begin{array}{c}{T}_{ij}=\{\begin{array}{lr}{K}_{ij}{l}_{ij}/{⟨d⟩}_{ij},& \mathrm{i}\mathrm{f}{\mathrm{\Omega }}_j\mathrm{i}\mathrm{s}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{o}\mathrm{f}{\mathrm{\Omega }}_i,\\ 2K_{ij}{l}_{ij}/{⟨d⟩}_{ij},& \mathrm{i}\mathrm{f}{\mathrm{\Omega }}_j\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{e}{\mathrm{\Omega }}_i,\end{array}\end{array}

(16c)\begin{array}{c}{K}_{ij}=\left(V_{\mathrm{p}i}+V_{\mathrm{p}j}\right)/\left(V_{\mathrm{p}i}/K_i+V_{\mathrm{p}j}/K_j\right),\end{array}

(16d)\begin{array}{c}{⟨d⟩}_{ij}={\int }_{{\mathrm{\Omega }}_i}\left|\left(x-x_j\right)\cdot n_j\right|\mathrm{d}x/{\int }_{{\mathrm{\Omega }}_i}\mathrm{d}x,\end{array}

where V_{\mathrm{p}i} and V_{\mathrm{p}j} are the pore volumes of {\mathrm{\Omega }}_i and {\mathrm{\Omega }}_j, respectively. This study adopts a matrix porosity of 0.3 and a fracture porosity of 1.0. {⟨d⟩}_{ij} (m) is the average normal distance from {\mathrm{\Omega }}_i to {\mathrm{\Omega }}_j. In Eq. (16d), x is the coordinate of any point in {\mathrm{\Omega }}_i, x_j is the coordinate of a fixed reference point in {\mathrm{\Omega }}_j, and n_j is the unit normal vector of {\mathrm{\Omega }}_j.

3) Fracture flux term

The fracture flux term (Fig. 4(d)) is discretized as:

(17a)\begin{array}{rl}& {\int }_{\mathrm{\partial }{\mathrm{\Omega }}_j}\left[-\frac{k_j}{\mu }\left(\mathrm{\nabla }{p}_j-\rho \mathrm{g}\right)\right]\cdot n_j\mathrm{d}s\approx \sum _{i=1}^k{T}_{ji}\left(u_j-u_i\right),\\ & \text{for}j\in \left[1,N_{\mathrm{f}}\right],{\mathrm{\Omega }}_i\in \text{neighbor}\left({\mathrm{\Omega }}_j\right),\end{array}

(17b)\begin{array}{c}{T}_{ji}=1/\left(1/T_i+1/T_j\right),\end{array}

(17c)\begin{array}{c}{T}_i=K_i{a}_i/d_i,\end{array}

(17d)\begin{array}{c}{T}_j=K_j{a}_j/d_j,\end{array}

where a_i and a_j are the apertures of {\mathrm{\Omega }}_i and {\mathrm{\Omega }}_j, respectively.

4) Fracture-matrix flow exchange term

Since the fracture-matrix flow exchange term equals the negative of the corresponding matrix-fracture flow exchange term, its numerical scheme is not reiterated here.

5) Fracture-fracture flow exchange term

The fracture-fracture flow exchange term (Fig. 4(e)) is discretized as:

(18a)\begin{array}{c}{\int }_{{\mathrm{\Omega }}_j}{q}_{jk}\mathrm{d}x\approx T_{jk}\left(u_j-u_k\right),\text{for}j\in \left[1,N_{\mathrm{f}}\right],{\mathrm{\Omega }}_j\cap {\mathrm{\Omega }}_k\ne \varnothing ,\end{array}

(18b)\begin{array}{c}{T}_{jk}=1/\left(1/T_j+1/T_k\right),\end{array}

(18c)\begin{array}{c}{T}_j=K_j{a}_j/d_j,\end{array}

(18d)\begin{array}{c}{T}_k=K_k{a}_k/d_k,\end{array}

where d_j (m) and d_k (m) are the average normal distances from the centroid of every subsegment of {\mathrm{\Omega }}_j and {\mathrm{\Omega }}_k to the intersection point between {\mathrm{\Omega }}_j and {\mathrm{\Omega }}_k (Fig. 4(e)).

Within the numerical scheme, the potential value at each cell center is treated as an unknown, resulting in a total of N_{\mathrm{m}}+N_{\mathrm{f}} unknowns. Applying the governing equation (Eqs. (13) and (14)) to each cell yields an equal number of equations, thus forming a well-posed system. Subsequently, the system is solved for the primary variable of pressure. The velocity field is then constructed using Darcy’s law (Eq. (3)). Finally, the total tunnel inflow rate is obtained by integrating the interfacial fluxes over all cells on the tunnel boundary. All computations are implemented numerically in MATLAB [49].

2.3 Model verification

2.3.1 Analytical verification

The proposed model is first verified against Lombardi’s analytical solution [50] for steady-state water inflow into circular tunnels (see Section S1 in Electronic Supplementary Material for details of Lombardi’s solution).

A benchmark case is established for a tunnel embedded in a fractured rock mass containing orthogonal sets of persistent and equally spaced fractures (Fig. 5(a)). The tunnel has a radius of 3 m and is positioned 30 m below the top boundary. Both fracture sets have an aperture of 3.75 × 10⁻⁴ m and a spacing of 6 m. The ME of the computational domain is 90 m. The verification case assumes: 1) impermeable rock matrix; 2) fixed water table at the top boundary; 3) constant water pressure at the bottom boundary.

Figure 5(b) compares the unit-length tunnel inflow rate (q) between the proposed model and Lombardi’s solution, showing excellent agreement with relative errors below 2% across all cases.

2.3.2 Numerical verification

Additional verification is performed under permeable matrix conditions using COMSOL Multiphysics simulations. The numerical verification case replicates the configuration of the analytical case (Fig. 5(a)), with the sole modification of assigning a matrix hydraulic conductivity of 9.8692 × 10⁻⁶ m/s.

The comparative results (Fig. 5(c)) demonstrate good consistency, with the relative errors around 5% for all cases. In summary, these comprehensive verification studies, covering both impermeable and permeable matrix conditions, confirm the reliability of the EDFM-based model for calculating water inflow into fractured rock tunnels.

3 Design of numerical experiments for coupled effects analysis

3.1 Selection of controlling parameters

As established in the Introduction, this study systematically evaluates the coupled effects of fracture spatial distribution and rock matrix permeability on tunnel water inflow. Consequently, the controlling parameters selected for analysis comprise fracture aperture (a), fracture spacing (s), fracture dip angle (\theta), and fracture-to-matrix permeability ratio (k_{\mathrm{f}}/k_{\mathrm{m}}). While other fracture characteristics such as the joint roughness coefficient may affect tunnel inflow behavior, these are reserved for future study.

The numerical analysis adopts a tunnel embedded in a fractured rock mass with orthogonal sets of persistent, equally spaced fractures (Fig. 5(a)), a configuration widely employed in theoretical studies [11,51].

3.2 Selection of parameter levels

To ensure comprehensive coverage of the parameter space, the parameter levels are determined following International Society for Rock Mechanics (ISRM) guidelines for rock mass discontinuity quantification [52], as listed in Tables 1 and 2.

Excessively large fracture apertures may lead to non-Darcy flow behavior [53]. To avoid such conditions, the aperture size is limited to below 2.5 mm, with a uniform value assigned to all fractures. Additionally, rock masses exhibiting small fracture spacing (indicative of crushed zones) are excluded by imposing a minimum spacing threshold of 200 mm in the numerical simulations. Given the orthogonality of the fracture sets, the orientation of the fracture network is fully determined by the dip angle of a single set. In line with established research protocols, k_{\mathrm{f}}/k_{\mathrm{m}} is set to 10³, 10⁴, 10⁵, 10⁶, and 10⁷ [19]. A control group with k_{\mathrm{f}}/k_{\mathrm{m}}={10}^{50} is also included to simulate impermeable matrix conditions. The selected parameter levels are summarized in Table 3.

A full factorial design [41] is adopted to systematically evaluate parametric influences, resulting in 384 unique cases. For all cases, the tunnel has a radius of 3 m and is positioned 20 m below the top boundary. To achieve accuracy-efficiency balance, the ME is optimized to 60 m (spacing level 1) and 100 m (spacing levels 2–4). The simulation assumes: 1) fixed water table at the top boundary; 2) constant water pressure at the bottom boundary.

4 Results and discussions

4.1 Influence of fracture spatial distribution on tunnel inflow across varying permeability ratios

4.1.1 Regression model equation

A set of quadratic polynomial regression models is established to characterize the functional relationship between tunnel inflow rate and fracture parameters across varying permeability ratios. The modeling procedure consists of the following stages.

1) Input parameter normalization

The input parameters of the regression model are fracture aperture, spacing, and dip angle. These input parameters are standardized using Z-score normalization [54]. This normalization step is crucial for stabilizing the regression algorithm and ensuring that all input parameters contribute equally to the model, regardless of their original units and scales. Specifically, Z-score normalization standardizes the data to have a mean of zero and a standard deviation of one, expressed as:

(19)X_{\mathrm{Z}}=\frac{X-{\mu }_X}{{\sigma }_X},

where X_{\mathrm{Z}} is the normalized parameter, X is the original parameter value, {\mu }_X is the sample mean, and {\sigma }_X is the standard deviation.

2) Numerical simulation

The 384 simulation cases (as described in Subsection 3.2) are performed using the EDFM-based tunnel inflow model, covering all combinations of fracture parameters and permeability ratios.

3) Model formulation

The quadratic regression model includes linear terms, quadratic terms, and all two-way interaction terms for the input parameters. It takes the form:

(20)\begin{array}{rl}\lg q=& x_1{a}_{\mathrm{Z}}+x_2{s}_{\mathrm{Z}}+x_3{\theta }_{\mathrm{Z}}+x_4{a}_{\mathrm{Z}}^2+x_5{s}_{\mathrm{Z}}^2+x_6{\theta }_{\mathrm{Z}}^2\\ & +x_7{a}_{\mathrm{Z}}{s}_{\mathrm{Z}}+x_8{a}_{\mathrm{Z}}{\theta }_{\mathrm{Z}}+x_9{s}_{\mathrm{Z}}{\theta }_{\mathrm{Z}}+x_{10},\end{array}

where q (m³/s/m) is the water inflow rate per unit tunnel length, a_{\mathrm{Z}}, s_{\mathrm{Z}}, and {\theta }_{\mathrm{Z}} are the normalized fracture aperture, spacing, and dip angle, respectively, and x_1–x_{10} are regression coefficients.

Six separate regression models are derived, corresponding to the six levels of permeability ratios in Table 3.

4) Parameter estimation

Model coefficients are determined using MATLAB’s “fitnlm” function [55], which implements iterative least-squares minimization to compute the coefficients. Results are summarized in Table 4.

5) Model validation

The goodness of fit of each regression model is assessed in Fig. 6 using two metrics: the relative error between the fitted and the numerically simulated inflow rate, and the adjusted coefficient of determination (R²). The x-axis (numbered 1 to 64) represents the index of the numerical experiments in the parametric study, with each number representing a unique combination of fracture parameters under a given permeability ratio. As outlined in Table 3, the four levels for each of the three fracture parameters result in a total of 4\times 4\times 4=64 numerical experiments per permeability ratio. The results demonstrate that the relative error is below 10% for most cases across all investigated permeability ratios. Furthermore, the adjusted R² exceeds 0.95 for every permeability ratio, indicating that each regression model can explain over 95% of the observed variance in tunnel water inflow. These metrics collectively affirm the robustness of these quadratic regression models.

Analysis of the regression coefficients in Table 4 indicates a convergence toward their respective asymptotic limits under the impermeable-matrix condition (i.e., k_{\mathrm{f}}/k_{\mathrm{m}}={10}^{50}). This trend aligns with physical expectations: as k_{\mathrm{f}}/k_{\mathrm{m}} increases, the contribution of matrix flow to total inflow diminishes and ultimately becomes negligible.

4.1.2 Individual effects of fracture parameters

1) Individual effect of fracture aperture

Across the studied range of permeability ratios, fracture aperture consistently dominates tunnel inflow compared to other fracture parameters, as demonstrated by the largest absolute values of its linear (\left|x_1\right|) and quadratic (\left|x_4\right|) coefficients among all parameter coefficients. As shown in Fig. 7(a) (with full results in Fig. S1 in Electronic Supplementary Material), inflow rate exhibits an exponential growth, spanning several orders of magnitude, with increasing aperture, irrespective of the permeability ratio. This profound influence is fundamentally governed by the cubic law for fracture flow, which establishes that the volumetric flow rate through a single fracture is proportional to the cube of its hydraulic aperture [46]. This nonlinear physical principle explains why even marginal changes in aperture can induce disproportionately large variations in inflow, thereby confirming aperture as the most influential factor in controlling tunnel water inflow.

2) Individual effect of fracture spacing

Fracture spacing exhibits an intermediate influence on water inflow, with the absolute values of its linear (\left|x_2\right|) and quadratic (\left|x_5\right|) coefficients being one order of magnitude smaller than those of aperture but two to three orders of magnitude larger than those of dip angle. As shown in Fig. 7(b) (with full results in Fig. S1 in Electronic Supplementary Material), inflow rate decreases with increasing spacing, following a nonlinear threshold behavior related to the percolation theory of fracture networks [56]. At small spacings, the well-connected fracture network establishes preferential flow paths and fracture dominates flow. Once spacing exceeds a critical threshold (approximately 4 m here), network connectivity deteriorates, which diverts a greater proportion of flow through the less permeable rock matrix, thereby reducing the sensitivity of inflow to further spacing increases.

3) Individual effect of fracture dip angle

Fracture dip angle demonstrates the weakest influence on tunnel inflow, as evidenced by the absolute values of its linear (\left|x_3\right|) and quadratic (\left|x_6\right|) coefficients, which are orders of magnitude smaller than those of aperture and spacing. This minimal influence is corroborated by Fig. 7(c) (with full results in Fig. S1 in Electronic Supplementary Material), which shows nearly constant inflow rates across the entire range of dip angles. The negligible effect of fracture dip angle stems from the configuration of the orthogonal fracture sets. Because of their orthogonality, the overall network permeability remains largely unaffected by the rotation of either set, as the directional sensitivity of one set is compensated by the orthogonally oriented set. Consequently, unlike systems with a single dominant fracture set where flow depends on orientation, this orthogonal fracture network makes aperture and spacing, rather than dip angle, the governing factors.

4.1.3 Coupled effects of fracture parameters on tunnel inflow

1) Aperture-spacing coupled effect

The aperture-spacing interaction exhibits the strongest coupled effect among all parameter pairs, as evidenced by its consistently highest regression coefficient (\left|x_7\right|) across most permeability ratios. This dominant coupled effect is visualized in Fig. 8 (with full results in Fig. S2 in Electronic Supplementary Material), showing: 1) enhanced sensitivity of inflow to aperture under smaller spacings; 2) increased sensitivity of inflow to spacing under larger apertures. This behavior can be attributed to the competition between fracture and matrix flow. 1) At small spacings, a well-connected fracture network provides efficient flow conduits. Here, aperture enlargement markedly increases inflow by enhancing these already-connected pathways. Conversely, at large spacings, poor network connectivity forces a greater reliance on matrix flow, and the low matrix permeability suppresses the impact of aperture variation. 2) The sensitivity to spacing is similarly modulated: large apertures enhance fracture dominance, making the system sensitive to spacing, whereas small apertures diminish fracture dominance and thus mitigate the role of spacing.

2) Aperture-dip angle coupled effect

The aperture-dip angle interaction shows consistently weaker coupling than the aperture-spacing coupled effect, evidenced by \left|x_8\right| remaining one to two orders of magnitude smaller than \left|x_7\right| across all permeability ratios. Figure 9 (see complete set in Fig. S3 in Electronic Supplementary Material) visually confirms this weak coupling through two key observations: 1) minimal variation in inflow rate with dip angle at any fixed aperture; 2) essentially parallel relationships between inflow and aperture across different dip angles. This weak coupling stems from the hydraulic compensation inherent to orthogonal fracture networks. The directional hydraulic sensitivity of one fracture set is systematically compensated by the orthogonally oriented set, regardless of aperture variations. Consequently, the overall rock mass permeability remains largely invariant to network rotation, effectively decoupling the influence of fracture aperture from that of dip angle.

3) Spacing-dip angle coupled effect

The spacing-dip angle interaction demonstrates secondary influence, with \left|x_9\right| values being one to two orders of magnitude smaller than \left|x_7\right| across most permeability ratios. Figure 10 (see full set in Fig. S4 in Electronic Supplementary Material) visually confirms this finding through two key observations: 1) negligible inflow variation with dip angle at fixed spacing; 2) nearly identical inflow-spacing relationships across different dip angles. This pattern originates from the same underlying mechanism that governs the weak aperture-dip angle coupling: the orthogonal fracture network configuration renders the system response independent of any single fracture set’s dip angle variations regardless of the spacing.

4.2 Coupled effects of fracture spatial distribution and matrix permeability on tunnel inflow

4.2.1 Individual effect of matrix permeability

Figure 11 illustrates the relationship between tunnel inflow and the permeability ratio across varying fracture apertures and spacings (see full set in Fig. S5 in Electronic Supplementary Material). The results delineate two distinct flow regimes: 1) for , inflow rate decreases markedly with increasing k_{\mathrm{f}}/k_{\mathrm{m}}, reflecting a substantial influence of matrix permeability on tunnel inflow; 2) for k_{\mathrm{f}}/k_{\mathrm{m}}>{10}^6, inflow rate becomes largely insensitive to further increases in k_{\mathrm{f}}/k_{\mathrm{m}}, indicating a negligible effect of matrix permeability. These observations align with the fundamental principle that as the permeability ratio rises, groundwater flow progressively localizes into fracture-dominated pathways. Consequently, the contribution of matrix flow to the total tunnel inflow diminishes significantly, becoming marginal under high k_{\mathrm{f}}/k_{\mathrm{m}} conditions.

4.2.2 Analysis of variance analysis

ANOVA [41] is utilized to quantify the coupled effects of fracture spatial distribution and matrix permeability on tunnel inflow. As a robust statistical method, ANOVA facilitates the decomposition of total variance into attributable components, enabling identification of factors that contribute significantly to observed variability.

As discussed in Subsubsection 4.2.1, a critical permeability threshold of k_{\mathrm{f}}/k_{\mathrm{m}}={10}^6 demarcates two distinct regimes of matrix permeability influence on tunnel inflow. Accordingly, separate ANOVA analyses are performed for low (k_{\mathrm{f}}/k_{\mathrm{m}}={10}^3,{10}^4,{10}^5) and high (k_{\mathrm{f}}/k_{\mathrm{m}}={10}^6,{10}^7,{10}^{50}) permeability ratio conditions to elucidate these differential mechanisms. The ANOVA analysis results are presented in Table 5.

1) Main effect

The analysis of individual factor effects (main effects) on tunnel inflow, as presented in Table 5, reveals three key findings: First, across all permeability ratios, both fracture aperture and spacing demonstrate statistically significant control (P < 0.01). Values of sum of squares (SS) further indicate that aperture exerts a greater influence compared to spacing. Second, while statistically significant (P < 0.01) at high permeability ratios, the SS value of fracture dip angle remains four orders of magnitude lower than that of aperture and spacing, confirming its relatively minor effect. Third, the influence of permeability ratio exhibits distinct behavior across different ranges: statistically significant at low ratios (P < 0.01) but negligible at high ratios (P ≈ 1). These findings demonstrate strong agreement with the regression analysis results presented in Subsubsections 4.1.2 and 4.2.1, thereby cross-validating the scientific robustness of the conclusions.

2) Coupled effect

ANOVA results in Table 5 elucidate the coupled effects controlling tunnel water inflow, revealing the following observations.

(i) a–s coupled effect

This coupled term exhibits statistical significance (P < 0.01) across all permeability ratios. The substantially higher SS value compared to other coupled terms confirms its dominant influence, consistent with the findings in Subsubsection 4.1.3.

(ii) a–θ and s–θ coupled effects

Although statistically significant at high permeability ratios, these coupled terms exhibit SS values approximately four orders of magnitude lower than those of the aperture-spacing coupled term, indicating secondary practical significance.

(iii) a–k_{\mathrm{f}}/k_{\mathrm{m}} coupled effect

The coupled term exhibits distinct threshold-dependent behavior, being statistically significant (P < 0.01) at low permeability ratios but insignificant (P ≈ 1) at high ratios. Figure 12 (see full set in Fig. S6 in Electronic Supplementary Material) confirms this pattern: at low k_{\mathrm{f}}/k_{\mathrm{m}}, the strongly nonlinear q–\lg \left(k_{\mathrm{f}}/k_{\mathrm{m}}\right) and q–a curves reflect pronounced coupling, whereas at high k_{\mathrm{f}}/k_{\mathrm{m}}, the nearly parallel curves indicate hydraulic decoupling. This threshold behavior reflects a shift in the dominant flow mechanism. At high permeability ratios, fracture flow dominates and the matrix contribution becomes negligible, rendering the system response independent of further changes in k_{\mathrm{f}}/k_{\mathrm{m}}. In contrast, at low permeability ratios, the rock mass behaves as a dual-porosity medium in which both fractures and the matrix contribute to flow. In this regime, decreasing the permeability ratio diverts a greater proportion of flow from the matrix into the fractures, amplifying the sensitivity of inflow to aperture variation at lower permeability ratios. Similarly, a larger aperture intensifies flow exchange during changes in k_{\mathrm{f}}/k_{\mathrm{m}}, enhancing the sensitivity of inflow to the permeability ratio.

(iv) s–k_{\mathrm{f}}/k_{\mathrm{m}} coupled effect

As shown in Fig. 13 (see full set in Fig. S7 in Electronic Supplementary Material), the coupled term exhibits a threshold-dependent behavior analogous to the a–k_{\mathrm{f}}/k_{\mathrm{m}} coupling, stemming from the same mechanism: at high permeability ratios, the matrix’s diminished contribution weakens the spacing-permeability ratio interaction. Sensitivity analysis confirms that: 1) inflow is more sensitive to fracture spacing at lower permeability ratios; 2) more sensitive to the permeability ratio at smaller spacings. This behavior can also be attributed to the dual-porosity nature of the rock mass at low permeability ratios. In this regime, a decrease in k_{\mathrm{f}}/k_{\mathrm{m}} transfers more flow from the matrix to the fractures, amplifying the impact of spacing variation on the rock mass permeability and tunnel inflow. On the other hand, at smaller spacings, the well-connected fracture network facilitates more efficient flow redistribution during changes in k_{\mathrm{f}}/k_{\mathrm{m}}, thereby enhancing the sensitivity of tunnel inflow to the permeability ratio.

(v) θ–k_{\mathrm{f}}/k_{\mathrm{m}} coupled effect

This term is statistically insignificant across all investigated conditions (P ≈ 1), a conclusion visually confirmed by the nearly parallel q–\lg \left(k_{\mathrm{f}}/k_{\mathrm{m}}\right) and q–θ curves in Fig. 14 (see full set in Fig. S8 in Electronic Supplementary Material), which indicate negligible coupled effects. This behavior stems from the same hydraulic compensation mechanism that characterizes other dip angle-related terms: the directional hydraulic sensitivity of one fracture set is systematically offset by the orthogonally oriented set, rendering the system response independent of any single fracture set’s dip angle regardless of the permeability ratio.

4.3 Dominant factors controlling tunnel water inflow

To quantitatively assess the relative contributions of various factors to tunnel water inflow, this study employs Partial Eta-squared ({\eta }_{\mathrm{p}}^2) as the effect size metric. This statistical parameter provides an accurate measure of both main effects and coupled effects on water inflow. The calculation formula of {\eta }_{\mathrm{p}}^2 is:

(21){\eta }_{\mathrm{p}}^2=\frac{{SS}_{\text{effect}}}{{SS}_{\text{effect}}+{SS}_{\text{total}}},

where {SS}_{\text{effect}} is the sum of squares for a given factor, and {SS}_{\text{total}} is the total sum of squares for all factors.

The calculated {\eta }_{\mathrm{p}}^2 values are presented in Fig. 15. The results demonstrate that across all permeability ratio conditions, the effect sizes of fracture aperture, aperture-spacing interaction, and spacing are substantially larger than those of other parameters, identifying them as the three primary determinants of tunnel inflow rate. In contrast, terms related to fracture dip angle show negligible effects. Under low permeability ratio conditions, k_{\mathrm{f}}/k_{\mathrm{m}}, a–k_{\mathrm{f}}/k_{\mathrm{m}}, and s–k_{\mathrm{f}}/k_{\mathrm{m}} display measurable effect size, confirming that matrix permeability significantly influences tunnel inflow in this interval. Conversely, at high permeability ratio conditions, all k_{\mathrm{f}}/k_{\mathrm{m}}-related terms demonstrate minimal effect size, indicating that tunnel inflow becomes predominantly governed by fracture spatial distribution, with matrix permeability exerting negligible influence.

4.4 Analysis of permeability ratio threshold

The coupled effect demonstrates distinct threshold-dependent behavior, making the determination of the permeability ratio threshold ((k_{\mathrm{f}}/k_{\mathrm{m}}{)}_{\text{crit}}) both fundamentally and practically significant. This section introduces a quantitative approach for identifying this threshold and systematically investigates its dependence on fracture parameters.

The methodology is based on the principle that beyond (k_{\mathrm{f}}/k_{\mathrm{m}}{)}_{\text{crit}}, the influence of matrix permeability becomes negligible. Adopting a 10% relative error, a widely accepted engineering tolerance, as the criterion, we define (k_{\mathrm{f}}/k_{\mathrm{m}}{)}_{\text{crit}} as the lowest permeability ratio at which the relative error induced by neglecting matrix permeability falls below the 10% tolerance. Using the case of k_{\mathrm{f}}/k_{\mathrm{m}}={10}^{50} as a benchmark for impermeable matrix, we compute the relative error in tunnel inflow resulting from the omission of matrix permeability across all simulation cases. The results are presented in Fig. 16 (see the full set in Fig. S9 in Electronic Supplementary Material).

1) Influence of fracture aperture on the permeability ratio threshold

As shown in Fig. 16(a), the relative error induced by neglecting matrix permeability decreases monotonically with increasing k_{\mathrm{f}}/k_{\mathrm{m}} across all aperture levels. For k_{\mathrm{f}}/k_{\mathrm{m}}>{10}^6, the relative error remains below 10% in all cases, indicating that matrix permeability can be neglected beyond this value. Moreover, the threshold (k_{\mathrm{f}}/k_{\mathrm{m}}{)}_{\text{crit}} decreases with increasing aperture. This observed trend can be explained by the cubic law [46], which states that an increase in aperture elevates fracture flow rate to the third power, thereby reducing the permeability ratio required for fractures to dominate the overall flow. Consequently, in systems with larger apertures, a higher matrix permeability is needed to significantly influence the rock mass permeability and tunnel inflow.

2) Influence of fracture spacing on the permeability ratio threshold

As observed in Fig. 16(b), the relative error declines monotonically with increasing k_{\mathrm{f}}/k_{\mathrm{m}} for all spacing values. When k_{\mathrm{f}}/k_{\mathrm{m}} exceeds 10⁶, the relative error remains below 10% in most configurations, suggesting that matrix permeability can generally be neglected beyond this threshold. Furthermore, (k_{\mathrm{f}}/k_{\mathrm{m}}{)}_{\text{crit}} decreases with reduced fracture spacing. This behavior stems from the fact that denser fracture networks form connected pathways more readily, allowing fractures to dominate the flow regime even at lower permeability ratios.

3) Influence of fracture dip angle on the permeability ratio threshold

As demonstrated in Fig. 16(c), the relative error decreases with rising k_{\mathrm{f}}/k_{\mathrm{m}} across all dip angles. Once k_{\mathrm{f}}/k_{\mathrm{m}}>{10}^6, the relative error falls below 10% universally, confirming the negligible role of matrix permeability beyond this point. Notably, (k_{\mathrm{f}}/k_{\mathrm{m}}{)}_{\text{crit}} exhibits minimal sensitivity to dip angle variation. This invariance aligns with the hydraulic compensation inherent to orthogonal fracture networks, as discussed in Subsections 4.1 and 4.2, which renders the global flow behavior largely independent of the orientation of individual fracture sets.

In summary, when k_{\mathrm{f}}/k_{\mathrm{m}} exceeds {10}^6, neglecting matrix permeability introduces less than 10% relative error in most scenarios, confirming that the influence of matrix permeability becomes marginal beyond this threshold. This finding is consistent with the results in Subsection 4.2, which also demonstrate negligible coupled effects above the same permeability ratio.

4.5 Engineering implications of the couple effects

The analysis results reveal significant coupled effects of fracture spatial distribution and rock matrix permeability on tunnel water inflow, which have direct consequences for tunnel engineering practice. The specific manifestations of these coupled effects are discussed below.

1) The coupled effects indicate that the priority of site investigation parameters is context-dependent. In a fracture-dominated system (e.g., sparsely fractured but aperture-rich rock with low permeable matrix), resources should be allocated to meticulously mapping fracture orientation, aperture, and connectivity. However, in a matrix-dominated system (e.g., intensely fractured but with a permeable matrix), determining the matrix permeability through in situ tests (e.g., packer tests) becomes more important than characterizing every single fracture.

2) The coupled effect between fracture aperture and matrix permeability dictates the flow distribution. When a large fracture aperture coincides with a high permeability ratio, it can lead to localized, high-velocity inflows that pose a risk for tunnel face instability. This necessitates pre-grouting or advanced drainage in these specific zones. In contrast, a more homogeneous, matrix-dominated flow (low k_{\mathrm{f}}/k_{\mathrm{m}}) results in widespread seepage, which may be managed effectively by a systematic drainage layer behind the lining rather than localized measures.

3) The coupled effect between fracture spacing and matrix permeability is meaningful for grouting measures. In rock masses with a low permeability ratio, the matrix is relatively permeable. In this scenario, an increase in fracture spacing may not lead to a proportional decrease in total inflow because the matrix compensates for the reduced fracture flow. This suggests that in such conditions, grouting the matrix may be as important as sealing fractures. Conversely, in rock masses with a high permeability ratio, fractures dominate flow. Here, the inflow is more sensitive to fracture spacing, and grouting efforts should prioritize treating the fracture network.

5 Case study

5.1 Project overview

The Zhuhai Connection Line constitutes a critical engineering component of the Hong Kong–Zhuhai–Macao Bridge infrastructure system, serving as a key transition between bridge and tunnel structures. As the sole extra-long tunnel along the connection line, Nanwan Tunnel represents a particularly challenging project.

Nanwan Tunnel traverses the Jiangjun Hill area between Nanping and Wanzai, a region characterized by tectonically denuded hilly terrain [57]. Peak elevations reach approximately 308 m, while the general elevation of the surrounding area ranges between 10 and 276 m. Natural slopes typically vary from 20° to 45°, with localized steep slopes exceeding 60°. The exposed strata within the tunnel zone are predominantly composed of granodiorite and biotite porphyritic granite. Surface vegetation is relatively well-developed and consists mainly of arboreal species, a factor that heightens the need for groundwater environmental protection. Groundwater in the area occurs primarily as pore water, weathered fracture water, and structural fracture water. Furthermore, the tunnel alignment crosses reservoir areas and fault fracture zones twice, resulting in a generally water-rich geological environment. These combined engineering geological and hydrogeological conditions impose stringent requirements on the waterproofing and drainage system design for Nanwan Tunnel.

This case study presents a scenario analysis of Nanwan Tunnel to demonstrate the practical application of the proposed model under site-specific conditions. It is important to note that due to the absence of comprehensive, continuous inflow monitoring data, this analysis is not intended as a direct model validation but rather as a plausible simulation based on the available geological and hydrogeological data. The primary purpose is to evaluate the coupled effects of stochastic fracture distribution and rock matrix permeability on the tunnel’s water inflow.

5.2 Numerical simulation parameters

The fracture network in the surrounding rock masses of Nanwan Tunnel demonstrates pronounced spatial heterogeneity and anisotropy. To statistically characterize the fracture network, we employ DFNs, which explicitly represent individual fracture and network topology [58]. These DFNs are semi-generic, meaning they are not exact field replicas but incorporate several field-measured statistical constraints [59], thereby enhancing the generalizability and applicability of our inflow analysis results. The geometric characteristics of the fracture network are defined by the following statistical distributions.

1) Fracture position

Fracture positions are determined by their geometric centroids, which are randomly distributed within the computational domain \mathrm{\Omega } following a uniform probability distribution:

(22)f\left(x,y\right)=\{\begin{array}{lr}{\displaystyle \frac{1}{S}},& \left(x,y\right)\in \mathrm{\Omega },\\ 0,& \left(x,y\right)\notin \mathrm{\Omega },\end{array}

where \left(x,y\right) represents Cartesian coordinates of fracture centroids, and S denotes the total area of the domain \mathrm{\Omega }.

2) Fracture length

Fracture length obeys a power-law probability density function [60]:

(23)f\left(l\right)=\frac{\alpha -1}{l_{min}^{1-\alpha }}{l}^{-\alpha }\text{for}l\in \left[l_{min},l_{max}\right],

where l is the fracture length, l_{min} and l_{max} define the length range, and \alpha is the characteristic exponent. According to the field measurements reported by Xiong et al. [61], the values of l_{min} and l_{max} are set as 1 and 60 m, respectively. The characteristic exponent \alpha is assigned a value of 2.5, which lies within the typical range for 2-dimensional rock masses (1.7–2.75) as documented in Sadeghnejad et al. [62].

3) Fracture orientation

Fracture orientation follows the Fisher distribution, which describes angular deviations from the mean orientation [63]:

(24)f\left(\theta \right)=\kappa \frac{\sin \theta e^{\kappa \cos \theta }}{e^{\kappa }-e^{-\kappa }},

where \theta is the angular deviation, and κ is the Fisher constant that controls the orientation concentration (higher κ indicates greater orientation concentration). For 2-dimensional fractures, fracture orientation is characterized by the dip angle. Based on the field measurements from Xiong et al. [61], the mean dip angle is set as 50°. The parameter κ is assigned a range of [0,8].

In summary, the distribution parameters of DFNs are assigned according to Table 6. The values of κ and P20 are varied to systematically examine their influence on tunnel inflow. Table 7 details the parameter combinations. For each case, we generate five independent DFN realizations to account for the inherent aleatory uncertainty in DFN modeling.

The rock matrix permeability is set to 1 × 10⁻⁸ m². The ME of the computational domain is set to 40 m. The tunnel has a radius of 3 m and is centered 40 m below the top boundary [61]. All computational cases assume: 1) fixed water table at the top boundary; 2) constant pressure at the bottom boundary. We utilize COMSOL’s 2-dimensional DFN plugin for DFN generation, while the water inflow is computed using the developed EDFM-based model.

We validate the model by comparing its predictions with COMSOL simulations under identical stochastic fracture network configurations from Table 7 (one DFN realization selected per computational case). As shown in Fig. 17, the results exhibit strong consistency, with relative errors remaining below 5% across all cases. This close agreement corroborates the reliability of the proposed model in handling stochastically distributed fracture networks, thereby enhancing its credibility and generalizability.

5.3 Calculation results

5.3.1 Maximum fracture-induced tunnel inflow

Figure 18 presents the maximum fracture-induced inflow for both individual and averaged DFN realizations across all computational cases. The simulated values range from 1 × 10⁻⁶ to 12 × 10⁻⁶ m³/s, showing reasonable agreement with benchmark data (1 × 10⁻⁶–2 × 10⁻⁶ m³/s) [61]. The systematic overestimation primarily results from unaccounted physical factors in the current model formulation, such as non-Darcy flow effects. Despite these limitations, the order-of-magnitude consistency with benchmark data demonstrates the potential of the EDFM-based model as an efficient preliminary assessment tool for tunnel inflow estimation.

As illustrated in Fig. 18, maximum fracture-induced inflow shows no systematic variation with changing κ or P20, indicating limited control by fracture orientation distribution and density.

5.3.2 Influence of fracture network distribution on tunnel inflow

Figure 19 (see the full set in Fig. S10 in Electronic Supplementary Material) presents the hydraulic pressure fields of DFNs with varying Fisher constants. The results reveal highly consistent pressure distributions across all κ values, indicating that fracture orientation distribution has a negligible influence on the seepage field around Nanwan Tunnel.

Similarly, Fig. 20 (see the full set in Fig. S11 in Electronic Supplementary Material) displays the hydraulic pressure fields of DFNs with varying P20 values. The observed hydraulic pressure fields show minimal variation across the tested P20 range, confirming that fracture density has a limited effect on the seepage field.

Table 8 summarizes the statistical characteristics of tunnel water inflow, including the mean, standard deviation, and coefficient of variation (CV), for different computational cases. The key observations are as follows: 1) the mean inflow rate varies by less than 1% across all cases; 2) the intra-case CV remains below 0.1 for all DFN realizations. Since the differences between computational cases arise solely from variations in fracture network configurations, these results further demonstrate that fracture spatial distribution plays a secondary role in governing the inflow into Nanwan Tunnel.

5.3.3 Influence of matrix permeability on tunnel inflow

The influence of matrix permeability on tunnel inflow is further investigated by conducting numerical simulations under the matrix-permeability-excluded condition. The simulations employ identical parameters to those in the matrix-permeability-included scenario, with the exclusion of matrix permeability to isolate its hydraulic contribution.

Figure 21 presents the statistical characteristics of tunnel water inflow under the matrix-permeability-excluded condition, including the mean value and CV. The results indicate that the mean water inflow rate ranges from 1 × 10⁻³ to 8 × 10⁻³ m³/s, representing a reduction of more than 10% compared to the matrix-permeability-included case. This substantial discrepancy underscores the critical role of matrix permeability in accurately estimating water inflow for the Nanwan Tunnel project. Moreover, Fig. 21 shows that the CV of water inflow rates increases significantly when matrix permeability is neglected, implying that matrix permeability helps stabilize inflow rate by reducing uncertainty induced by stochastic fracture distribution.

In summary, for the Nanwan Tunnel project, matrix permeability plays a predominant role in governing water inflow, with its hydraulic mechanisms operating through two primary pathways: 1) substantially increasing the mean water inflow rate; 2) effectively mitigating prediction uncertainty arising from stochastic fracture distribution. Conversely, the fracture spatial distribution exhibits only secondary influence on tunnel water inflow.

6 Conclusions

6.1 Principal conclusions

This study systematically investigates the coupled effects of fracture spatial distribution and matrix permeability on tunnel water inflow utilizing a newly developed EDFM-based tunnel inflow model. Quadratic regression and ANOVA are employed to quantify the coupled effects. Nanwan Tunnel is selected as a case study. The principal conclusions are summarized as follows.

1) The proposed EDFM-based model demonstrates strong agreement (relative errors generally below 5%) with both Lombardi’s analytical solution and COMSOL simulations, demonstrating its great potential for simulating water inflow into fractured rock tunnels.

2) Individual parameter effects. (i) Inflow rate exhibits a strong positive correlation with fracture aperture. (ii) Inflow rate decreases with increasing fracture spacing. However, the reduction rate diminishes once the spacing exceeds 4 m, indicating a weakened hydraulic influence under sparsely fractured conditions. (iii) Fracture dip angle has a negligible influence on water inflow. (iv) A critical permeability ratio threshold ((k_{\mathrm{f}}/k_{\mathrm{m}}{)}_{\text{crit}}={10}^6) governs the influence of matrix permeability, which is significant below this value and minimal above it. This threshold decreases with larger aperture and increases with wider spacing, yet remains nearly independent of fracture dip angle.

3) Coupled effects. (i) The aperture-spacing interaction dominates pairwise effects. The influence of spacing is amplified at larger apertures, while the effect of aperture is more pronounced at smaller spacings. (ii) In contrast, interactions involving the dip angle (aperture-dip angle, spacing-dip angle, and dip angle-permeability ratio) are statistically insignificant. (iii) Both the aperture-permeability ratio and spacing-permeability ratio interactions exhibit a threshold-dependent behavior, being significant only for . Within this regime, fracture aperture and spacing exert greater influence at lower permeability ratios, whereas permeability ratio gains prominence at larger apertures and smaller spacings, respectively. (iv) Fracture aperture and spacing, along with their interaction, are the dominant factors controlling tunnel inflow.

4) For Nanwan Tunnel, neglecting matrix permeability would result in an underestimation error of more than 10% in the calculated tunnel inflow, as it is a dominant factor. The influence of matrix permeability is twofold, by both elevating the mean inflow rate and reducing the uncertainty induced by stochastic fracture distribution.

6.2 Limitations and future work

The EDFM-based tunnel inflow model incorporates several deliberate simplifications to focus on the core physical mechanisms. The rationale and associated limitations are outlined as follows.

1) Tunnel model simplifications

(i) A circular cross-section is adopted for geometric symmetry, enabling focused analysis of the coupled effects. While absolute tunnel inflow values may vary with tunnel profiles, the identified global trends in coupled mechanisms remain valid, as these are governed primarily by rock mass properties rather than excavation shape.

(ii) Waterproofing and drainage systems are omitted to focus on the hydrogeological response of surrounding rock masses. This approach represents a conservative scenario, providing an upper-bound estimate of tunnel water inflow. The results thus offer a critical baseline for designing waterproofing and drainage systems in practical applications.

2) Fractured rock mass model simplifications

(i) Fractured rock mass models containing orthogonal sets of persistent and equally spaced fractures are used in parametric analysis to elucidate fundamental coupled effects. In rock masses with random fracture networks, groundwater flow is strongly influenced by network topology, a factor not sufficiently addressed in the present study. Nevertheless, the core coupled mechanism identified in this work, the threshold-dependent interaction between fracture spatial distribution and matrix permeability, is expected to remain conceptually valid in stochastic systems. Specifically, the coupled effect is likely significant when k_{\mathrm{f}}/k_{\mathrm{m}} falls below a system-specific threshold and becomes negligible beyond it. The precise value of this critical threshold and the magnitude of factor effects should be further studied in future work.

(ii) In the case study, site-based random fracture networks are established to investigate the coupled effects of stochastic fracture distribution and rock matrix permeability on water inflow into Nanwan tunnel. This analysis remains preliminary and site-specific. Future work would focus on generalizing the coupled behaviors under spatially variable fracture distributions.

3) Groundwater flow model simplifications

(i) A 2-dimensional flow model is employed to ensure computational tractability for parametric analysis. This simplification neglects 3-dimensional effects such as longitudinal geological variation, potentially leading to underestimated tunnel water inflow.

(ii) Steady-state flow is assumed to reflect equilibrium groundwater conditions, relevant for long-term tunnel operation. This approach does not capture transient drawdown during construction, where inflow initially peaks and gradually declines. Thus, the results may not be directly suitable for the initial excavation phase but represent the final, stabilized condition. However, the underlying EDFM framework is inherently extensible to transient flow because the governing equations can be extended to include time-dependent terms (e.g., storage effects). Simulating transient flow is an important direction for our future work.

(iii) Darcy flow is assumed, which is valid under common geological settings. However, in extreme conditions such as large-aperture fractures, non-Darcy flow (e.g., turbulence) may occur, which is beyond the scope of the current model.

Future work will focus on extending the current research toward more realistic conditions and broader applicability.

References

[1] Palmstrom A , Stille H . Ground behavior and rock engineering tools for underground excavations. Tunnelling and Underground Space Technology, 2007, 22(4): 363–376

[2] Zhu Y M , Zhou J J , Zhang B , Wang H C , Huang M Q . Statistical analysis of major tunnel construction accidents in China from 2010 to 2020. Tunnelling and Underground Space Technology, 2022, 124: 104460

[3] Liu J H , Li X J . Analytical solution for estimating groundwater inflow into lined tunnels considering waterproofing and drainage systems. Bulletin of Engineering Geology and the Environment, 2021, 80(9): 6827–6839

[4] Gong C J , Cheng M J , Fan X , Peng Y C , Ding W Q . Hydraulic fracturing-based analytical method for determining seepage characteristics at tunnel-gasketed joints. Journal of Central South University, 2025, 32(4): 1520–1534

[5] Gokdemir C , Li Y D , Rubin Y , Li X J . Stochastic modeling of groundwater drawdown response induced by tunnel drainage. Engineering Geology, 2022, 297: 106529

[6] Xu H , Li X J , Gokdemir C . Modeling and assessing the impact of tunnel drainage on terrestrial vegetation. Tunnelling and Underground Space Technology, 2021, 116: 104097

[7] Pujades E , Jurado A . Groundwater-related aspects during the development of deep excavations below the water table: A short review. Underground Space, 2021, 6(1): 35–45

[8] Shahbazi A , Saeidi A , Chesnaux R , Rouleau A . Effects of fracture-system geometrical parameters on the inflow rate into a tunnel in rock: A numerical modeling experiment. Quarterly Journal of Engineering Geology and Hydrogeology, 2023, 56(2): 618850

[9] Bi J F , Jiang H . Water infiltration in drained circular tunnels: An analytical solution and its simplification. Geotechnique, 2025, 75(2): 180–191

[10] Sun H Y , Le L , Dai Y M , Rui Y , Zhu H H , Li X J , Li X H , Xue G W . Groundwater environmental effects risk evaluation in mountain tunnel construction: A dynamic decision support system based on fuzzy two-dimensional cloud probability model. Tunnelling and Underground Space Technology, 2025, 156: 106276

[11] Farhadian H , Katibeh H , Huggenberger P , Butscher C . Optimum model extent for numerical simulation of tunnel inflow in fractured rock. Tunnelling and Underground Space Technology, 2016, 60: 21–29

[12] Javadi M , Sharifzadeh M , Shahriar K . Uncertainty analysis of groundwater inflow into underground excavations by stochastic discontinuum method: Case study of Siah Bisheh pumped storage project, Iran. Tunnelling and Underground Space Technology, 2016, 51: 424–438

[13] Sharifzadeh M , Karegar S , Ghorbani M . Influence of rock mass properties on tunnel inflow using hydromechanical numerical study. Arabian Journal of Geosciences, 2013, 6(1): 169–175

[14] Shahbazi A , Saeidi A , Chesnaux R , Rouleau A . An index based on fracture length and aperture to predict groundwater inflow rates in tunnels excavated in fractured-rock. Transportation Geotechnics, 2024, 45: 101217

[15] Shahbazi A , Chesnaux R , Saeidi A . A new combined analytical-numerical method for evaluating the inflow rate into a tunnel excavated in a fractured rock mass. Engineering Geology, 2021, 283: 106003

[16] Cesano D , Bagtzoglou A C , Olofsson B . Quantifying fractured rock hydraulic heterogeneity and groundwater inflow prediction in underground excavations: The heterogeneity index. Tunnelling and Underground Space Technology, 2003, 18(1): 19–34

[17] Bang S H , Jeon S W , Kwon S K . Modeling the hydraulic characteristics of a fractured rock mass with correlated fracture length and aperture: Application in the underground research tunnel at KAERI. Nuclear Engineering and Technology, 2012, 44(6): 639–652

[18] Rong G , Peng J , Wang X Z , Liu G , Hou D . Permeability tensor and representative elementary volume of fractured rock masses. Hydrogeology Journal, 2013, 21(7): 1655–1671

[19] Wang C C , Liu X L , Wang E Z , Wang M Y , Liu C . Dependence of connectivity dominance on fracture permeability and influence of topological centrality on the flow capacity of fractured porous media. Journal of Hydrology, 2023, 624: 129883

[20] Romano V , Bigi S , Carnevale F , Hyman J D , Karra S , Valocchi A J , Tartarello M C , Battaglia M . Hydraulic characterization of a fault zone from fracture distribution. Journal of Structural Geology, 2020, 135: 104036

[21] He CYao CShao Y LFan HZhou C B. Effective permeability of three-dimensional fractured rock with low fracture densities. Journal of Tsinghua University (Science and Technology), 2021, 61(8): 827–832 (in Chinese)

[22] Regnet J B , David C , Robion P , Menéndez B . Microstructures and physical properties in carbonate rocks: A comprehensive review. Marine and Petroleum Geology, 2019, 103: 366–376

[23] Huang Y , Yu Z B , Zhou Z F . Simulating groundwater inflow in the underground tunnel with a coupled fracture-matrix model. Journal of Hydrologic Engineering, 2013, 18(11): 1557–1561

[24] Ma L , Gao D , Qian J Z , Han D , Xing K , Ma H C , Deng Y P . Multiscale fractures integrated equivalent porous media method for simulating flow and solute transport in fracture-matrix system. Journal of Hydrology, 2023, 617: 128845

[25] Zhang K , Xue Y G , Xu Z H , Su M X , Qiu D H , Li Z Q . Numerical study of water inflow into tunnels in stratified rock masses with a dual permeability model. Environmental Earth Sciences, 2021, 80(7): 260

[26] Flemisch B , Berre I , Boon W , Fumagalli A , Schwenck N , Scotti A , Stefansson I , Tatomir A . Benchmarks for single-phase flow in fractured porous media. Advances in Water Resources, 2018, 111: 239–258

[27] COMSOL Multiphysics. Version 6.3. Stockholm: COMSOL. 2024

[28] 3DEC: Three-Dimensional Distinct Element Code. Version 7.0. Minneapolis, MN: Itasca Consulting Group, Inc. 2020

[29] Hokr M , Škarydová I , Frydrych D . Modeling of tunnel inflow with combination of discrete factures and continuum. Computing and Visualization in Science, 2012, 15(1): 21–28

[30] Yan C Z , Wang Y X , Xie X , Ali S , Sheng Z G . A 2D continuous-discrete mixed seepage model considering the fluid exchange and the pore pressure discontinuity across the fracture for simulating fluid-driven fracturing. Acta Geotechnica, 2023, 18(11): 5791–5810

[31] Zhao L H , Di Y T , Shao L Y , Mao J . An immersed finite-discrete element method (IFDEM) framework for water entry with fracture dynamics. Computer Methods in Applied Mechanics and Engineering, 2025, 442: 118026

[32] Ren F , Ma G W , Wang Y , Li T , Zhu H H . Unified pipe network method for simulation of water flow in fractured porous rock. Journal of Hydrology, 2017, 547: 80–96

[33] Yan X , Yu H T . Numerical simulation of hydraulic fracturing with consideration of the pore pressure distribution based on the unified pipe-interface element model. Engineering Fracture Mechanics, 2022, 275: 108836

[34] Ma G W , Wang H D , Fan L F , Chen Y . A unified pipe-network-based numerical manifold method for simulating immiscible two-phase flow in geological media. Journal of Hydrology, 2019, 568: 119–134

[35] Lee S H , Lough M F , Jensen C L . Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resources Research, 2001, 37(3): 443–455

[36] Hajibeygi H , Karvounis D , Jenny P . A hierarchical fracture model for the iterative multiscale finite volume method. Journal of Computational Physics, 2011, 230(24): 8729–8743

[37] Moinfar A , Varavei A , Sepehrnoori K , Johns R T . Development of an efficient embedded discrete fracture model for 3D compositional reservoir simulation in fractured reservoirs. SPE Journal, 2014, 19(2): 289–303

[38] Fumagalli A , Pasquale L , Zonca S , Micheletti S . An upscaling procedure for fractured reservoirs with embedded grids. Water Resources Research, 2016, 52(8): 6506–6525

[39] Rao X , Cheng L S , Cao R Y , Jia P , Wu Y H , He Y M , Chen Y . An efficient three-dimensional embedded discrete fracture model for production simulation of multi-stage fracture horizontal well. Engineering Analysis with Boundary Elements, 2019, 106: 473–492

[40] Zhang YRui YLi Z FLv Y YLi X J. Computing water inflow into fractured rock tunnels based on embedded discrete fracture model. Journal of Tongji University (Natural Science), 2025, 53(11): 1691–1701 (in Chinese)

[41] Jobson J D. Analysis of variance and experimental design. In: Jobson J D, ed. Applied Multivariate Data Analysis: Regression and Experimental Design. New York, NY: Springer, 1991

[42] Jiang X W , Wan L , Yeh T C J , Wang X S , Xu L . Steady-state discharge into tunnels in formations with random variability and depth-decaying trend of hydraulic conductivity. Journal of Hydrology, 2010, 387(3–4): 320–327

[43] Chen Z Q , He C , Zhang Y H , Xu Z Y , Li Z , Yu B X . The impact of formation heterogeneity on water discharge and groundwater depletion of an excavated tunnel. Journal of Hydrology, 2023, 627: 130403

[44] Nikvar Hassani A , Farhadian H , Katibeh H . A comparative study on evaluation of steady-state groundwater inflow into a circular shallow tunnel. Tunnelling and Underground Space Technology, 2018, 73: 15–25

[45] Atangana A. Fractional Operators with Constant and Variable Order with Application to Geo-Hydrology. Amsterdam: Elsevier Inc., 2018

[46] Snow D T . Anisotropie permeability of fractured media. Water Resources Research, 1969, 5(6): 1273–1289

[47] Blazek J. Principles of grid generation. In: Blazek J, ed. Computational Fluid Dynamics: Principles and Applications. 3rd Ed. Amsterdam: Elsevier, 2015

[48] Møyner O , Lie K A . A multiscale two-point flux-approximation method. Journal of Computational Physics, 2014, 275: 273–293

[49] Wong DDoster FGeiger S. Chapter 9: Embedded discrete fracture models. In: Advanced Modelling with the MATLAB Reservoir Simulation Toolbox. Cambridge: Cambridge University Press, 2021

[50] El Tani M . Circular tunnel in a semi-infinite aquifer. Tunnelling and Underground Space Technology, 2003, 18(1): 49–55

[51] Farhadian H , Katibeh H , Huggenberger P . Empirical model for estimating groundwater flow into tunnel in discontinuous rock masses. Environmental Earth Sciences, 2016, 75(6): 471

[52] ISRM (International Society for Rock Mechanics) . Suggested methods for the quantitative description of discontinuities in rock masses. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1978, 15(6): 319–368

[53] Rong G , Tan J , Zhan H B , He R H , Zhang Z Y . Quantitative evaluation of fracture geometry influence on nonlinear flow in a single rock fracture. Journal of Hydrology, 2020, 589: 125162

[54] Peng L B , Lu Z , Lei T , Jiang P . Dual-structure elements morphological filtering and local Z-score normalization for infrared small target detection against heavy clouds. Remote Sensing, 2024, 16(13): 2343

[55] MATLAB . Version R2023b. Santa Clara, CA: The MathWorks, Inc. 2023

[56] Sanderson D J , Nixon C W . Topology, connectivity and percolation in fracture networks. Journal of Structural Geology, 2018, 115: 167–177

[57] China Railway 14th Bureau Group Co.Ltd. Guanyue Road & Bridge Co.Guangdong Ltd . Monitoring and Measurement Plan for Nanwan Tunnel. 2013 (in Chinese)

[58] Gómez S , Sanchidrián J A , Segarra P , Bernardini M . A non-parametric discrete fracture network model. Rock Mechanics and Rock Engineering, 2023, 56(5): 3255–3278

[59] Hyman J D . Flow channeling in fracture networks: Characterizing the effect of density on preferential flow path formation. Water Resources Research, 2020, 56(9): e2020WR027986

[60] de Dreuzy J R , Méheust Y , Pichot G . Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (DFN). Journal of Geophysical Research: Solid Earth, 2012, 117: B11207

[61] Xiong FJiang Q HChen S YHu X C. Modeling of coupled Darcy–Forchheimer flow in fractured porous media and its engineering applications. Chinese Journal of Geotechnical Engineering, 2021, 43(11): 2043 (in Chinese)

[62] Sadeghnejad S , Masihi M , King P R . Dependency of percolation critical exponents on the exponent of power law size distribution. Physica A, 2013, 392(24): 6189–6197

[63] Huang N , Liu R C , Jiang Y J , Cheng Y F . Development and application of three-dimensional discrete fracture network modeling approach for fluid flow in fractured rock masses. Journal of Natural Gas Science and Engineering, 2021, 91: 103957

RIGHTS & PERMISSIONS

Higher Education Press