Introduction

During the construction of a tunnel, water inrush is one of the most common and complex geological disasters and has a large impact on the construction schedule and on safety (Coli et al. 2008); furthermore, when serious water inrushes occur in tunnel construction, huge economic losses and casualties can occur. Because water inrush causes great harm to a tunnel, the prediction of the groundwater inflow into a tunnel before and during the excavation process is an important task to ensure the safety and schedule during the underground construction process. Additionally, estimated pumping requirements and groundwater control measures are commonly based on predicted groundwater inflows to tunnel.

Inflows can be calculated by means of numerical or analytical solutions. In the last several decades, effort has been made on developing more sophisticated analytical solutions to estimate the water inflow in tunneling (e.g. Lei 1999; Cesano et al. 2000; Celestino et al. 2001; El Tani 2003; Kolymbas and Wagner 2007), but various simplifying assumptions were made in the analyses such as a horizontal tunnel, homogeneous and isotropic aquifer, and even steady-state flow conditions. Lei (1999) investigated the water inflow to a tunnel using analytical solutions while considering steady-state flow to the tunnel and reported some relationships for prediction of this phenomenon, whereas Cesano et al. (2000) investigated parameters regulating groundwater inflow to hard-rock tunnels with a statistical study of Bolmen Tunnel in southern Sweden. All of the methods used in these studies are limited to two-dimensional (2D) cross sections; thus, while useful for some purposes, these methods are not fully applicable to complicated real-world problems due to their limiting assumptions. In reality, a tunnel is drilled progressively rather than instantaneously; therefore, groundwater flow systems, accordingly, should evolve with the progressive drilling process. Recently, there has been increasing interest in estimating the groundwater flow to the tunnel as the drilling progresses. Maréchal et al. (2014) studied water inflow to mechanized tunnels, excavated with tunnel bore machines (TBMs), using analytical solutions during excavation in a heterogeneous unconfined aquifer. It is shown that most of existing analytical solutions overestimate the initial discharge due to the assumption that drilling was instantaneous over the entire tunnel length; however, analytical solutions have been developed to evaluate the transient-drilling speed-dependent discharge rates to a tunnel in homogeneous and heterogeneous formations (Perrochet 2005; Perrochet and Dematteis 2007; Maréchal et al. 2014). Perrochet (2005) investigated the amount of water inflow to a tunnel during excavation using an analytical solution that is based on the seepage of water to tunnels depending on the rate of excavation in an infinite homogeneous aquifer. Perrochet and Dematteis (2007) investigated an analytical model developed to predict transient discharge flow into a tunnel drilled at various speeds through a heterogeneous formation. This model relies on simplifying assumptions commonly enforced in hydrogeology engineering, and combines the convolution and superposition principles to account for composite sections with arbitrary parametric contrasts. An application of the data, measured during the exploratory drilling of an Alpine tunnel, confirms the validity of the approach. Yang and Yeh (2007) developed a mathematical model to describe the groundwater inflow to a tunnel in a multi-layer aquifer system. All these analytical solutions of groundwater discharge into a tunnel are based on the assumption that the drawdown of the water table is constant; however, Anagnostou (1995) found that the drawdown of the water table was considerably affected by the rate of excavation advance in some circumstances. Hwang and Lu (2007) proposed a semi-analytical method to analyze the tunnel water inflow with a constant flow and variable drawdown of the water level, which is a function of the inflow rate. In analytical methods, the information on geological structure and rock mass is usually simplified to a few simple input data, which are not used directly to predict the tunnel inflows from complex aquifers.

Three-dimensional (3D) numerical models can provide more accurate prediction, considering a progressive drilling process; by contrast, numerical models have rarely, if ever, been used for forecasting inflows during tunnel excavation (Wittke et al. 2006). Numerical models also can be used for more complicated groundwater problems such as modeling the effects of tunnel construction in heterogeneous, anisotropic, and bounded aquifers. Molinero et al. (2002) used TRANMEF-3 to conduct an efficient and accurate simulation of the transient hydrogeological conditions at and around a tunnel during the excavation process. Yoo and Asce (2005) developed a series of 3D stress–pore-water-coupled models in ABAQUS to examine the interaction mechanism between tunneling and groundwater. It is shown that the ground and lining responses are significantly influenced by the relative permeability of the lining, and that the circumferential pregrouting is an effective means for minimizing the tunneling and groundwater interaction. Also highlighted is the importance of the stress–pore pressure coupled analysis in the numerical prediction of tunnel behavior. Dassargues (1997) utilized MODFLOW for simulation of groundwater in order to predict the effect of excavating a tunnel with length of more than 500 m on the groundwater level and the behavior of the aquifer in Liege, Belgium. Zhang et al. (2002) used a dual porosity deformable jointed rock model for simulation of water flow in jointed rock using Universal Distinct Element Code (UDEC). Zhang et al. (2007) developed 3D finite-element numerical model with FEFLOW to simulate groundwater inflow during tunnel construction by considering that both material properties and hydraulic heads along the tunnel alignment change significantly and dynamically. With the Zhang et al. (2007) method, the tunnel was simulated as a time-variable head-dependent flux boundary (called ‘transfer boundary condition’ in FEFLOW). Results of the study showed that the proposed method can effectively simulate groundwater flow conditions during tunnel construction. Li et al. (2009) studied water inflow to two different types of tunnel in the path of Changli Expressway, southern China—the Qingshangang twin-tube and Bimaxi double-arch— using FLAC 3D numerical simulations and innovations regarding dewatering systems. Yang et al. (2009) investigated the effect of tunnel excavation on the hydrogeological environment in Tseng-Wen Reservoir Transbasin Diversion Project in Taiwan. In their research, GMS software together with two numerical codes of MODFLOW (for simulating groundwater inflow to the tunnel) and FEMWATER (for evaluating the impact of tunneling construction on the adjacent hydrogeological environment) and a geographic information system (GIS) software package were used for the analyses. In the simulation, location and properties of rock layers and different geological units, topography, permeability of different units, precipitation, groundwater recharge from stream-flow record, and groundwater level were considered. Font-Capó et al. (2011) tried to predict water inflow to a mechanized tunnel with a diameter of 12 m in granitic and metamorphic rocks in L9 subway line in Barcelona (Spain) using numerical simulations. In their research, influence of the groundwater system on the water inflow to the tunnel was simulated using Visual TRANSIN in three-dimensions. Regarding different geological and geophysical structures in the tunnel route, it was shown that the highest water inflow to the tunnel is in the regions with high frequency of fractures and underground barriers as dikes. Chiu and Chia (2012) simulated groundwater inflow in the long term during excavation of the Hsueh-Shan tunnel in relation to a water reservoir in northern Taiwan. In their work, the long-term effect of excavating a 40-km tunnel on the drainage of the water reservoir in dry and wet seasons, was examined using MODFLOW. Butscher (2012) suggests that numerical models provide estimates of tunnel inflows with sufficient accuracy for practical purposes if the tunnel is lined and has no drainage layer surrounding the lining, and if the hydraulic conductivity of the lining is several orders of magnitude lower than the hydraulic conductivity of the aquifer, or if the lining is thick. Preisig et al. (2014) studied the drainage of groundwater and surface subsidence due to excavation of mechanized tunnels. They predicted the water inflow to the La Praz tunnel in France using analytical solutions and, thereafter, considering the geological structures and hydrological condition of the tunnel route, while performing some numerical simulations for different zones of the tunnel. It was shown that numerically obtained results had higher precision with respect to the results obtained from analytical solutions. Farhadian et al. (2016) investigated the boundary effect and dimensions of a numerical model on the simulation accuracy of water inflow to tunnels in jointed rock masses. In their research, performed with an UDEC commercial software package, the input variables were tunnel radius, joint spacing, and horizontal and vertical extent of the model. Xia et al. (2017) studied the impact of permeability heterogeneity on groundwater inflow by using MODFLOW through a comparison between a homogeneous hydraulic conductivity case and a synthetic heterogeneous one. They used a dynamic modeling approach to simulate the change of groundwater flow step by step in accord with tunnel excavation.

The presented method in this paper was used for modeling groundwater inflow to a mechanized tunnel during the TBM advance, using a 3D finite-different numerical model developed with MODFLOW code of GMS. This method, in comparison with other available methods (e.g. Zhang et al. 2007; Font-Capó et al. 2011; Xia et al. 2017), is different in that all three tunnel seepage conditions (see Lai et al. 2017), material properties and hydraulic heads, change with the drilling advance in the numerical model. The method involves defining the tunnel boundary condition using the Drain package (DRN) in MODFLOW code.

Methodology for the simulation of the groundwater flow

Selection of modeling approach in rock mass

Development of an appropriate numerical model for simulation of groundwater flow in rock mass depends on joint conditions of rock mass and groundwater-modeling study scale. Groundwater flow in rock mass is simulated using three different approaches: equivalent porous media (EPM), dual porosity (DP), and discrete fracture network (DFN)—Shapiro (1987); Schwartz et al. (1990); Gupta and Singhal (1999); Li et al. (2015). Whenever possible, for the purpose of modeling groundwater in rock mass at large scale and when no joint sets have a significant effect on the flow behavior, the equivalent porous medium can be used (Yang et al. 2009). This common and established modeling approach for fractured sites, assumes that the fractured system behavior is equivalent to porous media behavior and can be represented by an equivalent porous medium, with equivalent hydraulic conductivity in a certain area (ITRC 2017). In addition, any fracture zones or fault zones can be modeled using a layer with high permeability, and may require a finer model grid design (Faille et al. 2015). The DP and DFN models are available through commercial, government, and academic sources, but are not used as commonly as equivalent porous media type models, even though these models offer advantages for modeling fractured rock systems (ITRC 2017).

Selection of appropriate modeling code

One of the most widely used programs for simulating groundwater flow in porous media is MODFLOW. The modular finite-difference groundwater flow model (MODFLOW) developed by the US Geological Survey (USGS) is a computer program for simulating common features in groundwater systems (McDonald and Harbaugh 1988; Harbaugh and McDonald 1996). The program was constructed in the early 1980s and has continually evolved since then with development of many new packages and related programs for groundwater studies. The main advantages of MODFLOW are that the code is robust, easy to use, and versatile. Also, it enables the modeler to extract detailed water balance information from the model (using the Flow-Budget subroutine), which greatly assists with model interpretation and trouble-shooting (Wels et al. 2012). There are several actively developed commercial and noncommercial-graphical-user interfaces for MODFLOW, which often include the compiled MODFLOW code with modification. Currently, the most comprehensive graphical user environment for performing groundwater simulations is GMS (Groundwater Modeling System). GMS software supports both finite-difference and finite-element models in 2D and 3D including MODFLOW, MODPATH, FEMWATER, SEEP 2D, MT3DMS/RT3D, SEAM3D, ARTED and UTCHEM (Environmental Modeling Research Laboratory 2005).

Numerical modeling workflow

With rapid increases in computation power and the wide availability of computers and model software, groundwater modeling has become a standard tool for professional hydrogeologists to effectively perform most tasks. A groundwater model is a replica of some real-world groundwater system (US Army Corps of Engineers 1999). Construction of the groundwater model involves several steps and requires inputs and participation from different professions (Refsgaard and Henriksen 2004). Workflow for construction of a groundwater flow model for application in forecasting is presented in Fig. 1.

Fig. 1
figure 1

Workflow for construction of the groundwater flow model (Anderson and Woessner 1992)

Full size image

Data collection procedure

Data sources for groundwater flow modeling should include geomorphology, geology, geophysics, climate, vegetation, soils, hydrology and hydrogeology information and data (Kolm 1996). The procedure of data collection for groundwater modeling studies is shown in Fig. 2. In this research, ArcGIS software was used for data storage, retrieval and analysis for conceptualizing the groundwater system and converting the data to the appropriate formats for the numerical model (Albertson and Hennington 1996).

Fig. 2
figure 2

Procedure of data collection for groundwater modeling studies (Qiu et al. 2015)

Full size image

Definition of the tunnel boundary condition

Simulation of groundwater inflow during tunnel advance is very challenging, difficult or even impossible. Both hydraulic properties and hydraulic heads along the tunnel route change dynamically as tunnel excavation advances. To resolve this problem, the tunnel should be specified as a time-variable head-dependent flux boundary to predict groundwater inflow rate and change over time as the tunnel excavation advances (Zhang et al. 2007). In the finite-difference 3D model of MODFLOW, DRN is used to simulate head-dependent flux boundaries (Harbaugh et al. 2000). According to the presented method, after the numerical model is properly constructed, the tunnel boundary condition is simulated using the DRN in transient condition. That is, groundwater inflow to the tunnel is determined by:

$$ Q=C\ \left(H\hbox{--} {H}^R\right) $$
(1)

where Q is the groundwater inflow to the tunnel [L3/T], C is the conductance of the tunnel [L2/T], which represents the resistance of flow into the tunnel and is determined by the hydraulic conductivity in the immediate vicinity of the tunnel, H is the simulated hydraulic head at the cells adjacent to the tunnel cell [L], and HR is the reference hydraulic head specified at the tunnel cell [L], which is equal to the tunnel elevation. In this step, according to TBM type and tunnel seepage condition, the conductance of the tunnel is set to change during TBM advance.

TBM types and tunnel seepage conditions

Tunnel boring machines (TBM) can be classified by the method for excavation (full face or partial face), the type of cutter head (rotation or non-rotation), and by the method of securing reaction force (from gripper or segment). Several types of TBMs are illustrated in Fig. 3 (ITA WG 14 2001; Wittke et al. 2006); the type of TBM is of paramount importance regarding the drilling method, lining, and waterproofing, with the most important TBMs generally being the open shield, single shield and double shield.

Fig. 3
figure 3

Classification of tunnel boring machines (ITA WG 14 2001; Wittke et al. 2006)

Full size image

An open-shield TBM (Fig. 4) is suitable for application in a rock mass in which support of the excavated cross-section in the area of the temporary face and of the machine is not required or may be achieved with minor additional engineering. Utilizing this type of TBM causes the groundwater to seep not only into the tunnel face, but also into the previously excavated sections (ITA WG 14 2001; Wittke et al. 2006).

Fig. 4
figure 4

A schematic view of the open-shield TBM and location of groundwater inflow (ITA WG 14 2001; Wittke et al. 2006)

Full size image

Single-shielded TBMs (Fig. 5) are utilized if the rock mass, because of its too low strength, is not able to carry the bracing forces of a gripper TBM, which would be necessary to transmit the required thrust forces. A single-shielded TBM without facilities for a face support can also be applied, if the excavation contour is not stable and if rock collapse may occur. The shield skin, which covers the entire machine, serves as a temporary support; however, for a final support, usually pre-cast lining segments of reinforced concrete are used. The lining segments are installed under the protection of the rear part of the shield, the so-called tail-skin, and the space between the shield skin and the excavation contour is referred to as steering gap. In unstable rock mass or soil, the steering gap may be close. The interface between the segmental lining and the excavation contour, the so-called annular gap, is normally grouted with mortar using injection lines which are integrated in the tail-skin. In Fig. 6, a double-shield TBM is illustrated, which, compared with a single-shield TBM, has the advantage, that for waterproofing the excavated sections, the segmental lining can be erected simultaneously with boring. When the segmental linings in the tunnel have been installed with gaps and offsets (between contiguous segments and between neighboring rings), and also the annular gap (the space between the segmental lining and the excavated ground) has been filled with pea gravel and grouting has been carried out in a later stage, the groundwater seeps not only into the tunnel face, but also into the previously excavated sections (ITA WG 14 2001; Wittke et al. 2006).

Fig. 5
figure 5

A schematic view of the single-shield TBM and location of groundwater inflow (ITA WG 14 2001; Wittke et al. 2006)

Full size image
Fig. 6
figure 6

A schematic view of the double-shield TBM and location of groundwater inflow (ITA WG 14 2001; Wittke et al. 2006)

Full size image

Defining the tunnel boundary condition

After the model calibration and verification phases, the model can be used as a tool for simulating groundwater inflows during tunneling, by conceptualization and definition of the tunnel boundary condition. The boundary condition of the tunnel can be defined using DRN via the following steps:

  1. 1.

    Creating the arc (polyline) of the tunnel route in GMS software. The beginning of the arc should match the tunnel inlet and its end should match the tunnel outlet (spatial coordinates and bottom elevation).

  2. 2.

    Dividing the arc into different intervals. The length of each interval is equal to the TBM excavation length in a given time period.

  3. 3.

    Setting the conductance of the tunnel intervals to zero before the tunnel excavation face reaches the location. This prevents groundwater inflows to the “not-yet-excavated” portion of the tunnel ahead of the TBM face (Zhang et al. 2007).

  4. 4.

    Considering the TBM type and the tunnel seepage condition. The conductance can be set to each of the tunnel intervals using one of three following scenarios (Fig. 7):

    • Scenario 1: In open TBM or where seepage into the tunnel is allowable as design. The conductance of the intervals is set to a specific value once the tunnel excavation face reaches the location, and remains unchanged thereafter.

    • Scenario 2: In single- or double-shield TBM, if seepage cannot be completely stopped by installed segmental linings. The conductance of the intervals is set to a specific value once the tunnel excavation face reaches the location. The conductance value is set to change to less than or equal to the specific value (depending on the severity of the lining leakage) once the tunnel excavation face has passed by, and remains unchanged thereafter.

    • Scenario 3: In single- or double-shield TBM if full waterproofing be applied to segmental linings. The conductance of the intervals is set to a given value once the tunnel excavation face reaches the location. The conductance value is set to change to zero once the tunnel excavation face has passed by, and remains unchanged thereafter.

Fig. 7
figure 7

A schematic plot of the tunnel-boundary-condition scenarios

Full size image

Once the tunnel boundary is defined and assigned to the 3D finite difference grid, the model can be run in transient condition. The computer program ZONEBUDGET (Harbaugh 2005) can be applied to calculate the discharge rates into the tunnel. ZONEBUDGET uses cell-by-cell flow data saved by MODFLOW to calculate the water budgets. The observed values of groundwater inflow to the tunnel may be compared to predicted values. If there is reasonable correspondence between predicted and observed values, confidence in the modeling may be established. If there is significant deviation between predicted and observed values, it may be necessary to redefine the model (Wels et al. 2012). The flowchart of the definition of the tunnel boundary condition is given in Fig. 8.

Fig. 8
figure 8

Flowchart of the definition of the tunnel boundary condition

Full size image

Study region of the Qomroud tunnel

The Qomroud water conveyance tunnel is one of the components of a water management system in Iran and involves a 36-km tunnel from the River Dez to the Golpayegan Dam reservoir (Fig. 9). The tunnel was originally divided into four parts, each about 9 km, and was put out to bid as design/build contracts in 2002. The study region for this paper is the first part of this tunnel beginning with the inlet of the tunnel at coordinates x = 377,338, y = 3,683,810 and z = 2061.20, and ending at coordinates x = 385,552, y = 3,681,920 and z = 2046.97. The Sabir Company was the successful bidder for the first part of this project, which was constructed using a Herrenknecht 4.65-m-diameter earth pressure balance (EBP) shield TBM. The concrete lining was created using 5 + 1 tetragonal reinforced concrete segments in each ring. The tunnel hasbeen excavated at a grade of 0.134% and finished with a concrete segmental lining to a diameter of 3.8 m.

Fig. 9
figure 9

a Geographical location and b satellite image of the first part of the Qomroud tunnel

Full size image

Geological condition of the first part of the Qomroud tunnel

According to the structural-sedimental zoning of Iran, the study region is located in the Sanandaj-Sirjan zone (Alavi 1994; Mohajjel and Fergusson 2000; Ghasemi and Talbot 2006), which consists of a series of asymmetric folding and faults and has experienced mild to high metamorphisms which have caused schistosity and recrystallization of minerals. The lithology of the area consists of a sequence of Precambrian formations comprised of schist and slate (Sh) and hard rocks like limestone and metamorphic dolomites (D1) and sometimes quartzite rocks (Qz) and metavolcanic rocks (An); also, Quaternary alluviums comprise soft alluvial zones including saturated fine-grained soils (Q).

Tectonic activities have caused formation of various faults with thick shear and crushed zones along the tunnel route. Most of these faults are of reversed type. In addition, some normal and strike slip faults have been observed in the area. Based on the geological studies, the Moghanak fault zone was recognized along the tunnel route with width of 20 m at 5 + 270 km. The rock mass is classified as “very poor” by standard systems. The geological plan and profile of the first part of the Qomroud tunnel is shown in Fig. 10 with a scale of 1:2000.

Fig. 10
figure 10

The a geological plan and b profile of the first part of the Qomroud tunnel (scale 1:2000). Elevation in m above sea level (asl)

Full size image

Process of the groundwater flow simulation

In this research, after delineation of the model domain to development of the conceptual model, climatological and hydrological studies were performed to calculate the rate of groundwater recharge from precipitation, evaporation from groundwater level and groundwater discharge rate to the rivers. Geological studies were also conducted for development of a stratigraphic model of the model domain, and hydrogeological studies were done to determine bedrock depth of the alluvium, and to assign the initial and boundary conditions, hydraulic properties of the stratigraphic units, and the impact of applied stresses on the groundwater system (such as extraction wells). Using the abovementioned information and the conceptual model, the numerical model was constructed using MODFLOW code of GMS software. After that, the model was calibrated in steady-state condition and a sensitivity analysis was performed for recognizing the uncertainties of the input parameters to the model; then, after defining the tunnel boundary condition using the DRN, the calibrated model was used for prediction of the groundwater inflow to the tunnel during TBM advance. Finally, in order to evaluate the precision and accuracy of the presented method, the simulated rates were compared with the observed rates of groundwater inflow to the tunnel.

Hydrogeological investigation

Model domain

The delineation of the model domain is dependent on the selection of suitable external boundary conditions. It is generally preferable to use physical hydrological features which are known to control groundwater flow such as watershed divides, lakes, aquitards, etc. (Wels et al. 2012). In order to determine the boundaries of the numerical model, the aquifers that are influenced by the tunnel excavation should be recognized and separated from the other aquifers. Ridge lines and the river in each drainage basin can be seen as natural boundaries. The ridge lines of the drainage basin are defined as boundaries with no-flow, while rivers are considered to be constant head boundaries (Yang et al. 2009). According to the drainage basins division of the Water Resource Agency of Iran, the first part of the Qomroud tunnel is located in the drainage basins of the Anuj and the Aligudarz rivers (Fig. 11). Therefore, the boundaries of the Anuj and the Aligudarz drainage basins are defined as no-flow boundaries, except for the northwestern and southern borders, which are the rivers flowing out of the drainage basins. These parts of the boundaries are considered to be constant head boundaries (Fig. 12).

Fig. 11
figure 11

Drainage basins in which the first part of the Qomround tunnel is located

Full size image
Fig. 12
figure 12

Domain of the first part of the Qomroud tunnel model. Other watershed boundaries are considered to be no-flow

Full size image

Precipitation infiltration

In the hydrological reports of the Anuj and the Aligudarz drainage basins, using the water budget method, the rate of areal recharge of precipitation, which is one of the most important boundary condition for the numerical model, is estimated. This is based on the mean annual areal recharge in these drainage basins of about 25 mm and its rate is 7 × 10−5 m/d (Table 1)—that is, only 6% of the precipitation is infiltrated to the groundwater system.

Table 1 The water budget results of the Anuj and the Aligudarz drainage basins. MCM  million cubic meters
Full size table

Groundwater discharge to the rivers

Rivers are major surface-water features that are perennial sources or sinks for groundwater. Several methods have been developed to estimate groundwater discharge and recharge from river flow records (Lee et al. 2006). The most commonly used methods are base-flow separation techniques. As shown in Fig. 11, the Anuj River is the location of the beginning of the first part of the tunnel and the Aligudarz River is located close to the end of it; thus, there are two rivers in the model region, but only the Anuj River has a gauge station. The average measured flow rate at the Anuj gauge station is about 0.72 m3/s, based on the base-flow separation method. About 45% of it is runoff directly from precipitation and the rest is groundwater discharge (known as base-flow). According to the experts’ advice, conductance of the riverbeds of both rivers is considered equal to 0.3 m2/d/m for the groundwater modeling. The conductance of the riverbeds was adjusted during calibration until simulated discharge matched measured flow at the Anuj gauge station, and used thereafter (Anderson et al. 2015).

Groundwater evapotranspiration

If the water table is close to the land surface, there may be direct evaporation from the water table and/or phreatophytes (plants whose roots extend into the water table) may extract groundwater through transpiration (Anderson et al. 2015). In groundwater models, the maximum evapotranspiration rate and extinction depth (depth of the dominant plant roots) are usually required.

Because of the shallow depth of the groundwater level in the lands in the margin of the Anuj and the Aligudarz rivers, which is mostly farmland, these lands are considered to be the groundwater evapotranspiration zones. Maximum evapotranspiration rate of these lands is considered equal to 9.0 × 10−4 m/d.; additionally, since the maximum rooting depth of plants in the farmlands is 1.0 m, the evapotranspiration extinction depth is considered equal to 1.0 m (Fig. 13).

Fig. 13
figure 13

Location of farmlands in the study region

Full size image

Stratigraphic units

Considering the complicated stratigraphy of the model domain, and for development of the stratigraphic model in its simplest form, the following were considered:

  1. 1.

    According to the lithology of rocks in the study region, the stratigraphic units are divided into six different groups: schist, slate, metadolomite, marble, quartzite, and andesite-rhyolite.

  2. 2.

    Fault zones are considered to be units with very high permeability (Faille et al. 2015) and a thickness of 200 m, taking into account spatial discretization of the stratigraphic model for the 3D grid.

  3. 3.

    Trends of stratigraphic units along the tunnel route are continued to the model boundaries in order to neglect the stratigraphic unit complications due to minor faults. The extent of alluvial units is extracted from the geological map and their thicknesses are measured from the boreholes along the tunnel route. The stratigraphic units in the study region are presented in Fig. 14.

Fig. 14
figure 14

Stratigraphic units in the study region

Full size image

Permeability of stratigraphic units

Altogether, 26 boreholes were drilled in the first part of the tunnel, the location of which is shown in Fig. 15. The results of Lugeon tests conducted in these boreholes are shown in Table 2. According to Table 2, the values of permeability for schist, slate, quartzite and andesite-rhyolite are equal to 1 Lugeon (9.5 × 10−5 m/d), metadolomite and marble units are equal to 3 Lugeons (2.9 × 10−4 m/d), and permeability of the fault zones are considered equal to 100 Lugeons (9.5 × 10−3 m/d). The permeability of the Moghanak fault zone was determined by Lugeon tests performed in borehole BH-5, close to the fault zone. Lugeons can be related to conventional SI units, as a Lugeon is approximately equal to 1.1 × 10−7 m/s (Banks et al. 1992).

Fig. 15
figure 15

The location of the boreholes along the first part of the Qomroud tunnel

Full size image
Table 2 Permeability of stratigraphic units according to Lugeon tests
Full size table

The permeability measurements in the alluvial units include Lefranc tests (constant and falling head) in the boreholes and single-well pumping tests in three extraction wells (Fig. 16). According to the results of the Lefranc and pumping tests, the permeability of the alluvial units is considered to be 5 × 10−2 m/d.

Fig. 16
figure 16

a Location and (b–d) the results of single-well pumping tests in three extraction wells close to the tunnel route

Full size image

Specific yield

Specific yield (dimensionless), also known as the drainable porosity (Lohman 1970), is the storage property of the unconfined aquifer and is only required for the transient condition of the model. According to the results of the pumping tests, the specific yield of the alluvial units is equal to 0.2 (20%), and based on the experts’ advice, the specific yield values of the rock units and fault zones are respectively considered to be 0.01 (1%) and 0.26 (26%).

Discharge from groundwater sources

According to the data from periodic monitoring of the groundwater sources in the study region from 2003 to 2011, 30 wells, 13 springs, and 14 aqueducts (the Persian Qanats) have been recognized, the locations of which are shown in Fig. 17. The flow rates of the wells range from 4 to 36 L/s and their depths vary from 14 to 100 m. The mean discharge of the springs was 0.1–13 L/s, and the mean discharge of the aqueducts was 4–19 L/s. For simulation of the aqueducts and springs using the DRN in the MODFLOW code, the conductance was considered equal to 0.05 and 1 m2/d/m, respectively. The conductance of the aqueducts and springs was adjusted during calibration by comparing simulated discharge to measured flow (Anderson et al. 2015).

Fig. 17
figure 17

Location of the groundwater sources

Full size image

Groundwater level

In order to specify the initial head along the tunnel route, groundwater level was measured in the boreholes during the planning and design phases of the project. Because of the lack of groundwater level data within the model domain, first, a regression equation was developed using existing land-surface and water-table elevations for the boreholes (Fig. 18). The regression equation was then used to estimate the spatial distribution of groundwater level in the model domain using the mean land-surface elevation values available from the SRTM 90 m (Shuttle Radar Topography Mission) digital elevation model (DEM; USGS 2007). The potentiometric surface map of the model domain is shown in Fig. 19.

Fig. 18
figure 18

The relationship between groundwater level and ground elevation

Full size image
Fig. 19
figure 19

Groundwater-level elevation in the model domain

Full size image

Development of the conceptual model

To solve any site-specific groundwater problem, the hydrogeologist must assemble and analyze relevant field data and articulate important aspects of the groundwater system. The synthesis of what is known about the site is a conceptual model (Kresic and Mikszewski 2012). In general, the closer the conceptual model approximates the field situation, the more likely the numerical model will give reasonable forecasts. Key components of a conceptual model include boundaries, hydrostratigraphy and estimates of hydrogeological parameters, general directions of groundwater flow and sources and sinks of water, and a field-based groundwater budget (Anderson et al. 2015). The conceptual model is usually presented visually as a series of 2D cross-sections or as a 3D block diagram, with supporting text and data tables that describe and quantify the components and features of the model (Wels et al. 2012). The block diagram of the conceptual model for the study region is shown in Fig. 20.

Fig. 20
figure 20

Block diagram of the conceptual model for the study region. Other domain boundaries are defined as no-flow

Full size image

Design of the numerical model

Development of the 3D stratigraphic model and spatial discretization

The starting point in designing the real computer model is discretization of the model extent, on which the continuous natural system is subdivided to smaller segments (i.e., cells, elements, and blocks) and as a result, the model becomes ready for numerical solutions of the partial differential equations. In the MODFLOW code of GMS software, it is possible to use the 3D stratigraphic model to define the hydraulic conductivity and storage properties of each stratigraphic unit and automatically assign them to the 3D finite difference grid.

The Boreholes module of GMS software can be used to create the boreholes by using the borehole drilling logs and creating user-defined cross sections between existing boreholes. The cross section will show the stratigraphy of the site between two boreholes. As soon as the cross sections between all the boreholes are created, the ‘solid’ module can be used for constructing the 3D stratigraphic model. The top elevation of this model corresponds to the SRTM 90 m DEM (USGS 2007) and the bottom elevation of it is the 1,800-m level. In the next step, the 3D stratigraphic model should be spatially discretized. In order to reduce the numerical dispersions and errors which come from decomposition of the velocity vector of groundwater flow in the direction of the grid coordinate axes, the grid is oriented parallel to the principal direction of groundwater flow. The principal flow direction may be controlled by structural features such as faults and joint, strike of the hydrogeological units, or stratification (Wels et al. 2012). Since the dominant structural and stratigraphic direction of the study region is the Azimuth 140°, a 50° rotation is specific to the ‘grid frame’. In order to create the 3D grid, an optimization should be performed regarding the selection of cell sizes, not to be too large to neglect simulating the steep hydraulic gradients of the model nor too small to unreasonably increase the time and number of calculations and create rounding and cut errors in the model. Therefore, considering the aforementioned points and with the model domain area in mind (703 km2), cell dimensions were adopted as 100 × 100 m2. In this way, the model has 360 rows (i), 320 columns (j), and one layer (k). The 3D stratigraphic model and the 3D finite difference grid of the study region are shown in Fig. 21.

Fig. 21
figure 21

The 3D a stratigraphic model and b finite difference grid of the study region

Full size image

Definition of the boundary conditions

After constructing the grid, the boundary conditions (external and internal) have to be implemented into a numerical model. The boundary conditions can be defined in two different ways according to the grid and conceptual model approaches for the grid in the MODFLOW model of GMS software. The grid approach involves directly assigning the boundary conditions to the 3D grid on a cell-by-cell basis; however, in the conceptual model approach, all of the boundary conditions are defined in the model using feature objects in the ‘map’ module, and GMS automatically assign the values to the cells of a grid. Except for simple problems, the conceptual approach is typically the most effective.

In this section, MODFLOW packages are first applied for definition of the model boundary conditions (Fig. 22), and finally, they are assigned to the 3D finite difference grid (Fig. 23). All the used boundary condition packages are as follows:

  • Internal boundary conditions (stress packages)

  1. 1.

    Well package: for simulating the extraction wells

  2. 2.

    Areal Recharge package: for simulation of precipitation infiltration

  3. 3.

    River package: for simulating the Anuj and Aligoudarz rivers

  4. 4.

    Evapotranspiration package: for simulating the evaporation from the groundwater surface and the transpiration by plants from the underlying root zone

  5. 5.

    Drain package: for simulating the springs and aqueducts

  • External boundary conditions

  1. 1.

    Specific Head package: for simulation of a specific head for the outer boundary of the model

  2. 2.

    No-flow Boundaries: not defining any boundary conditions for the external boundary of the model in GMS software means having no-flow boundaries without requiring any package or assumptions.

Fig. 22
figure 22

Definition of the boundary conditions for the numerical model using the proper packages in the MODFLOW code

Full size image
Fig. 23
figure 23

Assigned boundary conditions to the 3D finite difference grid

Full size image

Definition of the calibration targets

During model calibration, field observations, including heads and fluxes, are used as calibration targets, which are compared to simulated equivalent values computed by the model (Anderson et al. 2015). Measured groundwater level in the boreholes and measured springs and aqueducts discharge, along with the Anuj River flow rate at the gauge station, are considered to be observation points/flows in the model (Fig. 24).

Fig. 24
figure 24

Definition of the calibration targets for the numerical model using the proper packages in the MODFLOW code

Full size image

Definition of the initial condition

In the steady-state condition, in order to accelerate convergence of the numerical solutions, the initial groundwater levels are defined as the initial conditions for the model, and for the transient condition, the initial conditions should preferably be adopted from the steady-state simulation (US Army Corps of Engineers 1999). For this purpose, the potentiometric surface maps are used for definition of the initial conditions for the numerical model (Fig. 25).

Fig. 25
figure 25

Application of the groundwater level as the initial conditions for steady-state simulation of the MODFLOW code

Full size image

Correcting the errors and running the model

After defining the boundary and initial conditions and hydraulic properties of the stratigraphic units for the grid in MODFLOW, the model should be checked with Model Checker. Because of the significant amount of data required for a simulation for the model, it is often easy to neglect important data or to define inconsistent or incompatible options and parameters. Such errors will either cause the model to crash or to generate an erroneous solution. The purpose of the Model Checker code is to analyze the input data currently defined for a model simulation and report any obvious errors or potential problems. Running Model Checker successfully does not guarantee that a solution will be correct. It simply serves as an initial check on the input data and can save a considerable amount of time that would otherwise be lost tracking down input errors. Finally, after finding and correcting of the errors, the numerical model can be run in the proper way.

Model calibration

Calibration is the process of adjusting model inputs to achieve a desired degree of correspondence between the model simulations and the natural groundwater flow system (US Army Corps of Engineers 1999). In other words, in order for a groundwater model to be used in any type of predictive role, it must be demonstrated that the model can successfully simulate observed aquifer behavior. In the calibration procedure, some of the input parameters of the model such as recharge and hydraulic conductivity are altered in a systematic fashion and the model is repeatedly run until the computed solution matches field-observed values (hydraulic head and flow) within an acceptable level of accuracy.

In GMS software, the calibration procedure of the model can be performed in different ways, as a manual method with trial and error, or completely automatically using inverse modeling with PEST (model-independent Parameter ESTimation) for MODFLOW. PEST is a general purpose parameter estimation utility developed by Doherty (2007) to assist in data interpretation, model calibration and predictive analysis. The calibration of the model can be performed in steady state, transient, or both of the conditions.

Because periodic groundwater level was not measured in the boreholes along the tunnel route, calibration of the model was performed only for steady-state condition. As a result, after defining the calibration targets in the conceptual model, in order to save time in calibration, the PEST utility is used to reach the calculation error and, thus, achieve green color for the calibration target’s bar in the GMS software (Fig. 26; Table 3).

Fig. 26
figure 26

Calibration results of the numerical model

Full size image
Table 3 The final results of the calibrated input parameters of the model
Full size table

Sensitivity analysis of the model

A sensitivity analysis is a quantitative evaluation of the influence on model outputs from variation of model inputs. A sensitivity analysis identifies those parameters most influential in determining the accuracy and precision of model predictions (US Army Corps of Engineers 1999).

In MODFLOW code of GMS software, a parameter sensitivity plot is used to display the sensitivity of the PEST parameters (hydraulic conductivity and areal recharge). This plot is automatically and easily created in the Plot Wizard by setting the plot type to ‘parameter sensitivity’. The result of the sensitivity analysis of the PEST parameters of the model is shown in Fig. 27. According to the shown graph, the model was more sensitive to the areal recharge than hydraulic conductivity of the stratigraphic units. Recharge from precipitation is often one of the greatest uncertainties in groundwater modeling studies and has to be evaluated further during mathematical model calibration and sensitivity analysis (Wels et al. 2012).

Fig. 27
figure 27

Diagram of sensitivity analysis of the PEST parameters of the model

Full size image

Scenario definition and prediction of the water inflow to the tunnel

After successfully performing the calibration and sensitivity analysis of the model, by defining the tunnel boundary condition, the data can be utilized to predict groundwater inflow to the tunnel, which is the objective of the modeling. The data of the excavation rate of the TBM is presented in Table 4 for every 3 months corresponding to the tunnel as-built. The first part of the Qomroud tunnel was constructed using EBP shield TBM, and seepage occurred throughout the previously excavated sections—because the segmental linings were installed with the gaps and the offsets, and also the annular gap was backfilled using pea gravel (Fig. 28), thus the tunnel boundary condition is defined with the first scenario. Therefore, after creating the arc (polyline) of the tunnel route in the GMS software, the arc was divided into 14 intervals, so that the length of each interval equals the length of the tunnel excavation in 3 months. Thereafter, the conductance is assigned to each interval separately using the XY Series editor (a simple list of time/data pairs) in the DRN as follows:

  1. 1.

    The conductance of the intervals, before the tunnel excavation face reaches the location, was set to zero.

  2. 2.

    The conductance of the intervals was set to 0.1 m2/d/m once the tunnel excavation face reached the location, and remained unchanged thereafter. The conductance of the tunnel excavation face and the excavated portion of the tunnel, based on the experts’ advice, was adopted as 0.1 m2/d/m.

Table 4 The data of the excavation rate of the TBM of the first part of the Qomroud tunnel
Full size table
Fig. 28
figure 28

The groundwater inflow to the tunneling excavation: a at chainage 2 + 015 km, b at chainage 8 + 271 km

Full size image

The definition of the tunnel boundary condition for the model, using the presented method, is shown in Fig. 29. After defining the tunnel boundary, the model was run in the transient condition. The beginning of the simulation is 29 October 2007 and the ending is 27 April 2011; additionally, the length of the time steps is defined as 3 months. The accumulated groundwater inflow to the tunnel in each time step of the simulation is presented in Table 5. The total amount of groundwater inflow to the tunnel has been estimated to be about 33% of the mean annual precipitation infiltration in the model domain.

Fig. 29
figure 29

Definition of the tunnel boundary conditions for the model with the presented method

Full size image
Table 5 Simulated vs. observed groundwater inflows during TBM advance
Full size table

Comparison of predicted and observed data

Comparing the predicted groundwater inflow during TBM advance and the observed rates from the as-built inflow, one can evaluate the precision and validity of the presented method. The observed accumulated groundwater inflow to the first part of the Qomroud tunnel is presented in Table 5. This table shows that the total amount of groundwater inflow to the tunnel is about 70 L/s which is 28% of the mean annual precipitation infiltration. Comparing the simulated and the observed groundwater inflow, it can be deduced that the prediction precision is about 80%, which indicates an acceptable precision of modeling with the presented method. Moreover, Table 5 shows that at 5 + 270 km, which is the location of the damage zone of the Moghanak fault, the observed groundwater inflow to the tunnel is 33 L/s, while the predicted value is 36 L/s. Both the predicted and observed accumulated groundwater inflow during tunnel excavation are compared in Fig. 30. As can be seen, the predicted rates are slightly higher than the observed rates, but with similar trends for the curves. The difference between predicted and observed rates may come from uncertainties in the input parameters and boundary conditions. Uncertainty is typically reduced by utilizing additional information collected during the construction phase and will become available for updates of the model.

Fig. 30
figure 30

Comparing the observed and simulated accumulated groundwater inflow to the first part of the Qomroud tunnel

Full size image

Conclusion

The presented method for simulation of the groundwater inflow into a tunnel at a large scale, implemented using a MODFLOW code, is a procedure which predicts groundwater inflow to the mechanized tunnel during TBM advance. A more important fact, which should be considered in this simulation model by defining the tunnel boundary condition using DRN in MODFLOW code, is that the tunnel seepage condition, material properties and hydraulic heads change during the TBM advance in the numerical model. In this process, two issues were considered: (1) inflows were concentrated at the face and the rock-machine contact because of the theoretical imperviousness of the tunnel due to the lining installation, which restricted the main entry of water to a “moving interval”, and (2) there was a possibility of some residual leakage in the lining. These two considerations led the researchers to adopt a variable leakage boundary condition, in contrast to open tunnels (Molinero et al. 2002), where a transient flow state without restrictive leakage is applied. The results of the simulation show that maximum discharge rate occurred when the high-permeability zone (e.g. fault zone or alluvial unit) was exposed. Also, comparing the simulated and the observed groundwater inflow curves indicates that the method works effectively in transient simulation of groundwater flow conditions for tunnel construction and produces realistic and rational results. In addition, the presented method is more comprehensive compared to the ones performed by previous researchers as Dassargues (1997), Yang et al. (2009), Zhang et al. (2007a), Font-Capó et al. (2011), Chiu and Chia (2012), Preisig et al. (2014) and Xia et al. (2017), and the obtained results have an acceptable precision and accuracy.