Prediction of tunnel squeezing in soft sedimentary rocks by geoelectrical data

Squeezing is time-dependent deformation that can cause technical difficulties and financial consequences in underground structures. This study employs electrical resistivity data to predict squeezing intensity along the Beheshtabad tunnel in the Sanandaj–Sirjan zone in Iran. For comparison analysis, a semi-empirical approach was correlated with numerical modeling to predict tunnel squeezing at the design stage. The squeezing intensity obtained for the Beheshtabad tunnel was then compared with the instability observations along the Golab tunnel excavated in sedimentary rocks of the same zone. We developed a relationship between electrical resistivity and strains and provided a new strain-based squeezing classification system. The calibrated electrical data produced more accurate results for predicting tunnel squeezing than the conventional methods. The results show that rock type, joint properties, and water saturation impact squeezing.

The opening radius of tunnels (m) E c Conglomerate K s,sh

Introduction
Time-dependent rock deformations by disturbance of the primary stress field can impose high pressure on support installations, especially in deep tunnels excavated in soft rocks (Swannell et al. 2016). This type of deformation, called squeezing, can continue during the construction phase or even over a long period. This phenomenon is crucial for mechanized tunneling (Barla and Barla 2008).
Squeezing is more common in relatively deep underground structures excavated in soft rocks. However, this phenomenon can occur (due to topographic or tectonic conditions) for weak rocks at shallow depths (Shrestha and Broch 2008) and in deep tunnels excavated in hard rocks (Malan 1999). Predicting squeezing potential around underground structures is challenging, especially in the design phase, and requires a deep understanding of rock mass behavior.
Squeezing problems have been recorded in many tunnels, e.g., the San Bernardino and Gotthard road tunnels in Switzerland (Thut et al. 2000) and the second part of the watersupply Karaj tunnel in Iran (Khanlari et al. 2012). Some previous studies also included comprehensive databases for various squeezing cases (Feng and Jimenez 2015;Sun et al. 2018;Chen et al. 2020).
Geotechnical equations and exploration approaches can provide a comprehensive perception of rock mass conditions and evaluate deformations of excavated tunnels. Geotechnical (conventional) methods are classified as empirical and semi-empirical equations (Hoek and Marinos 2000;Palmstrom 1996) and physics-based numerical modeling (Hasanpour et al. 2014). These equations are usually obtained based on a few geotechnical parameters that can affect squeezing. For example, the Aydan approach was derived from some tunnels excavated in Japan (Aydan et al. 1996). Applying geotechnical equations can provide erroneous results when geotechnical parameters are uncertain (i.e., in the design phase).
Empirical equations are developed based on rock mass classification indices (Singh et al. 1992;Goel et al. 1995), while semi-empirical equations are provided based on deformations. Empirical equations determine squeezing occurrence based on rock mass indices, e.g., using the Norwegian rock mass classification index (Singh et al. 1992). Thus, their results cannot be compared to semi-empirical approaches, which can predict the squeezing intensity. Semi-empirical equations employ different ratios: e.g., the uniaxial compressive strength of rock mass to vertical in situ stress (Jethwa et al. 1984), the calculated strain to critical strain (Aydan et al. 1996), the rock mass index (RMi) to induced tangential stress (Palmstrom 1996) and the uniaxial compressive strength of rock mass to strains (Hoek and Marinos 2000).
Physics-based numerical modeling can cover a large-scale medium. However, a proper numerical approach is required to simulate coupled processes to predict squeezing potential. Numerical modeling was previously used to evaluate deformations and squeezing, e.g., by incorporating elastoviscoplastic models into numerical tools to simulate timedependent deformations (Barla et al. 2012;Hasanpour et al. 2014). The accuracy of numerical models can be affected by uncertain geotechnical data and modeling simplifications, particularly at the design stage.
The exploration approaches (including laboratory experiments, in situ tests, and geophysical prospecting) can also be used to evaluate squeezing. In situ tests can cover larger areas comparing laboratory experiments, i.e., the tested area is excavated galleries instead of drilled boreholes. However, laboratory experiments and in situ tests can cover the study site partially.
Geophysical data can be a practical approach to determining the rock mass condition at a large scale. This approach has been employed to assess site characteristics of roads, dams, offshore installations, slopes, and ancient sites (Akinrinmade et al. 2013;Whiteley 2006). Geoelectrical data have not been considered as an index to detect underground deformations. Geoelectrical prospecting (i.e., relying on electrical resistivity contrasts of different geomaterials) was frequently used (Ammar and Kamal 2018).
The geotechnical and geophysical parameters (especially ER) are usually momentary and time-independent. However, the main difference between the ER and geotechnical data is the coverage area; i.e., geotechnical properties are based on drilled boreholes and can be uncertain at the design stage.
Complementary geotechnical data at the production phase usually alter the initial geotechnical data.
One challenging task in predicting squeezing was the time-dependent behavior of rock mass. All previous studies used indices such as rock mass classification systems (Jimenez and Recio 2011;Dwivedi et al. 2013) or geotechnical parameters (Aydan et al. 1996;Hoek and Marinos 2000) as an indirect scale of the time-dependent behavior. The driving force behind this indirect application was the possibility to define time-dependent deformation using fractures properties (e.g., weathering, orientation, spacing, filling materials, roughness, and water inflow) and geotechnical characteristics (including rock mass uniaxial compressive strength and overburden). Using joint properties and geotechnical conditions to predict squeezing at the design stage was limited to numerical modeling, e.g., correlating a simple numerical model and a semi-empirical equation to predict the squeezing (Khanlari et al. 2012).
The squeezing intensity is substantially related to rock composition and structure, and rocks with significant clay minerals (e.g., claystone, shale, siltstone, mudstone, marlstone) and a water-bearing porous medium are susceptible to tunnel squeezing (Kovari and Staus 1996;Aceh 2019;Arora et al. 2020).
The relationships between rock strength parameters and squeezing were studied (Aydan et al. 1996;Vrakas et al. 2017;Wood 1972). Under high overburden pressure, the ratio of rock strength to in situ stress decreased, and the squeezing potential increased (Gutierrez and Xia 2008;Jethwa et al. 1984;Hoek and Marinos 2000). Field observations revealed the direct effect of groundwater conditions and stress on squeezing (Arora et al. 2020;Sun et al. 2018).
In the case of fracture effects, a few squeezing approaches were developed based on the frequency and properties of discontinuities (Singh et al. 1992;Goel et al. 1995;Jimenez and Recio 2011;Dwivedi et al. 2013). Saturated weak rocks can deteriorate over time. This weakening process can also be related to weathering. Similar to weathering, the initiation and development of cracks can increase the squeezing intensity (Agan 2015).
Other factors that can affect squeezing, albeit to a lesser extent (Sun et al. 2018), are the installed support system as well as methods and rates of excavation, mainly related to decisions at the production phase (Sun et al. 2018), and are out of the scope of this research.
The prementioned factors can alter ER. The ER value of rocks with high clay compositions or porous rocks (with high water content) is low (Kovari and Staus 1996). Laboratory experiments, field data, and numerical modeling confirmed the relationship between parameters affecting ER and squeezing (Ammar and Kamal 2018) (Fig. 1). For example, the fluid-bearing capacity (Bhatt and Jain 2014), fluid pressure, weathering, fracture characteristics, faults activations (Backstrom 2004;Laszlo 2011), porosity, and lithology of rock matrix (Aceh 2019) can influence ER, similar to tunnel squeezing. For saturated rock specimens, alterations of ER are related to high vertical stress during experimental studies in plastic deformation zone (Brace and Orange 1968;Dayuan et al. 1988). However, lower ER values show a higher squeezing intensity.
Water exploitation by geoelectrical prospecting is a frequently used method (Ammar and Kamal 2018;Rolia and Sutjiningsih 2018). At the deformation zone, ER changes are due to the volume increase under shear deformations and micro-crack development (Chen and Lin 2004). The decreasing trend of ER due to fracture propagation, including the effects of crack frequency, was also detected (Li et al. 2015;Sandler et al. 2009).
The higher the overburden, rock strength, and fracture density are, the lower the ER is. Therefore, the relationship between these parameters and squeezing intensity is direct. Water saturation can also decrease ER. A combination of porosity, water-bearing capacity, structure, and mineralogy affect ER (and rock squeezing).
This study aims to predict the time-dependent behavior of rock mass based on ER, which may provide remarkable insight into the general situation of joints and rock mass. The in-common factors affecting squeezing and ER confirmed the correlation between rock squeezing and ER.
The chainage of 7.6-16 km of the Beheshtabad tunnel in the Sanandaj-Sirjan zone is prone to the squeezing phenomenon. We compared the predicted squeezing with observational data (obtained from chainage of 0-8.8 km of Fig. 1 Factors controlling (before excavation) ER and squeezing are presented. The factors decrease ER while increasing the squeezing intensity the Golab tunnel excavated in the sedimentary rocks of the same geological zone). Then, we developed a relationship between ER and squeezing based on in situ strain measurements of strain during the excavation of the Golab tunnel. This research also compared the prediction accuracy of the combined semi-empirical and numerical modeling to the ER-based approach.

Methodology
Two approaches were chosen to predict the squeezing potential along sedimentary rocks of the Beheshtabad and Golab tunnels: (1) a combination of numerical modeling and a semi-empirical approach and (2) geoelectrical prospecting. The observational data obtained during the excavation of the Golab tunnel were also employed for the comparison analysis. A relationship was defined between calibrated electrical data and in-the-field measured strains recorded in excavating the Golab tunnel, and the squeezing intensity along the Behehstabad study site was predicted by the developed approach.

Description of the study sites
The Beheshtabad tunnel is located in the Chaharmahal-Bakhtiari province in Iran, planned to convey water from the Beheshtabad dam to the Central Plateau of Iran. The tunnel has a northeast trend and ends in the Zayanderoud River upstream of Cham Aseman dam. The length (total) and diameter of the Behehstabad tunnel are 64.97 km and 6.8 m, respectively.
The Beheshtabad study site, chainage 7.6-16 km, is located in the Sanandaj-Sirjan zone, in the Zagros and Naien-Baft Orogen of Gondwanan provenance (Fergusson et al. 2016;Ghasemi and Talbot 2006). This zone (as a metamorphic belt) is affected by several geological and tectonic events: continental collision (Late Eocene-Early Miocene), Pan-African orogenesis, Late Palaeozoic rifting forming Neo-Tethys, Mesozoic convergence, and ophiolite obduction (Fergusson et al. 2016;Ghasemi and Talbot 2006). Figure 2a shows the location of the Beheshtabad tunnel as the dashed blue line and the study site as the red rectangle. The topographic and geophysical profiles of the study site are shown in Fig. 2b.
The Sanandaj-Sirjan zone (with a width of 150-200 km) extends 1500 km from Sanandaj (in north) to Sirjan city in the south. Paleozoic volcanic rocks, including Silurian, Devonian, and Permian rocks with intrusions of Mesozoic rocks and the metamorphism due to Cimmerian movements, are the main features of this zone. The zone is evolved dynamically by the Neotethys Ocean movements at the northeastern margin of Gondwana (Ghazi and Moazzen 2015), mainly composed of sedimentary rocks, including limestone, marlstone, siltstone, mudstone, shale, and conglomerate (Zayandab Engineering Consulting Company 2006, 2010. To estimate the behavior of the Beheshtabad study site during excavation, a chainage of the Golab tunnel (0-8.8 km), excavated in sedimentary rocks of the Sanandaj-Sirjan zone, was evaluated. Several instability failures, tunnel collapses, and squeezing damages were recorded in this tunnel; these observations were employed to predict the squeezing along the study site of the Beheshtabad tunnel (planned to be excavated).
The total length and diameter of the Golab tunnel are 17.02 km and 4.6 m, respectively. The support system was installed 1 m away from the face. The temporary and permanent support systems of the Golab tunnel were a 5 cm shotcrete layer reinforced with wire mesh and a 25 cm lining, respectively. In sections with severe squeezing problems, steel sets were also added to the temporary support system (Zayandab Engineering Consulting Company 2006). The mechanized tunneling using an open-shield TBM was planned for the Beheshtabad tunnel. The suggested support system was a 10 cm shotcrete layer reinforced with wire mesh as the temporary and a 25 cm lining as the permanent support system (Zayandab Engineering Consulting Company 2010). Therefore, the support systems of both study sites were designed almost similarly. Figure 3 illustrates geological profiles and drilled boreholes along the Golab and Beheshtabad tunnels. The axis of both tunnels is located below the water table. Eight boreholes were drilled along the Golab tunnel ( Fig. 3a) and four boreholes along the Beheshtabad tunnel (Fig. 3b).
For chainage 0-8.8 km of the Golab tunnel, various squeezing problems were recorded for formations with a high amount of clay, including siltstone and shale. A comprehensive set of data was recorded during the excavation of this tunnel along 2 m sections, including the joint properties, fault activations, groundwater condition, weathering, electrical resistivity, and squeezing intensity. The water table is about 150 m above the tunnel axis.
At the chainage of 0-1120 m, a conglomerate formation (E c ) (good rock based on RMR index) and hydraulic flux of 1.5 × 10 −7 to 1 × 10 −6 m/s was detected. The squeezing was limited to light cases (from the tunnel inlet to about 500 m). This squeezing was mainly tectonic-based, related to the branches of Fault Abrizan-Kamason. Water penetration (due to relatively high hydraulic conductivity) was also recorded (Zayandab Engineering Consulting Company 2006).
Light squeezing phenomena occurred in chainage 1210-1550 m, mainly related to fracture zones at limy shale. The squeezing in K lsh 4 (in sections 3050-6570 m, 6800-8100 m, and 8570-8800 m) was also recorded, especially at 6312-6570 and 7560-8100 m moderate squeezing cases were mainly related to fault activations (Zayandab Engineering Consulting Company 2006).
Along the Beheshtabad tunnel, chainages 7600-8120 m and 9630-10,790 are Razak marlstones ( M m r ) with relatively low flux (1 × 10 -8 -1 × 10 -10 m/s). The RMR class is fair to good. Possible risks of squeezing can be related to weathering of marlstones and the fault activations. At 7700 m, a fault was detected. Fault K2 can impact its surrounding marly rock mass at 9600-10,000 m (Zayandab Engineering Consulting Company 2010).
From 8120 to 9630 m, Asmary limestone ( OM L a ) (fair to good rock mass) is affected by a branch fault at 8400 m and Fault K2 at 9200-9600 m. Other parts of the rock mass do not show symptoms of weathering based on boreholes, and the flux of rock mass ranges between 1 × 10 -9 and 1.7 × 10 -7 m/s (Zayandab Engineering Consulting Company 2010).
10,790-12,110 m chainage mainly consists of Bakhtiary conglomerate ( PL-Q c b ), including sandstone ( Q s b ) at 10,790-11,100 m. The rock mass condition is fair to good, with no sign of weathering. The hydraulic flux of rock mass is relatively higher, 6 × 10 -8 to 1 × 10 -6 m/s, and water penetration can be possible. The most squeezing-prone part of this section is around 12,100 m, where a branch fault was detected (Zayandab Engineering Consulting Company 2010).
The remaining chainage (12,110-16,000 m) of the Beheshtabad study site is mainly composed of mudstone ( PL-Q ms b ) with a flux of 1.5 × 10 -7 -2.7 × 10 -7 m/s. At 15,600-16,000 m, conglomerate rock mass can increase the permeability. Fault Zagros impacts the chainage of 15,500 to 16,000 m (Zayandab Engineering Consulting Company 2010).
The existence of faults and rock masses with a high amount of clay shows that squeezing problems can be significant along the Beheshtabad tunnel. The groundwater level ranges between 40 and 60 m below the land surface (about 200 m above the tunnel axis).

Geophysical studies
Geophysical evaluations are essential for geoengineering projects, especially for deep tunnels. The geoelectrical surveys provide several advantages: (1) detecting geological, hydrogeologic, and geotechnical obscures. (2) Preparing electric resistivity maps. (3) Exploring tectonic factors and  Table 1. The profiles along the boreholes (Fig. 3a) present the RQD values along their lengths recognizing relations between physical and geological conditions. (4) Detecting faults/fractured zones.
The same effective parameters decrease ER while increasing the squeezing intensity (considering the literature review). Before defining a relationship between these two parameters, we should respond to two main challenges: (1) can we trust the accuracy of ER? (2) How a quantitative relationship can be defined between ER and squeezing?
We considered a novel idea: predicting tunnel squeezing based on ER variations. Nevertheless, the accuracy of this approach can be affected, because the structure depth (especially for deep tunnels) can affect the resolution of electrical prospecting (Bernabini et al. 1988). A long electrode distance can also undermine the efficiency of geoelectrical data for deep structures. Meanwhile, subsurface electromagnetic interferences (including conductive mineral masses, underground caves, and water-bearing layers) can influence the data quality due to being detected as sharp resistivity anomalies (Militzer et al. 1979;Rolia and Sutjiningsih 2018). A single geophysical method can recognize all these anomalies as potential squeezing areas. Therefore, some actions were required to evaluate the data accuracy and to cover the disadvantages of electrical prospecting before predicting squeezing.
Geoelectrical mapping of the Golab site includes a Schlumberger array (Fig. 4a) along three main profiles and continuous resistivity profiling (CRP) along seven directions. For the main profiles, the maximum flow line was 1000 m, the number of electrodes on each profile was 173, and the electrode distance was 50 m. The multi-electrode method with three electrodes was used for CRP at a 50 m distance (Zayandab Engineering Consulting Company 2006).
At the Beheshtabad study site, the Schlumberger array was performed along three main profiles (C, D, and E) in Fig. 2b; the CRP profiles, along eight directions. The distance of the main profiles was 250 m, the maximum flow line was 1000 m, and the electrode distance ranged from 50 to 100 m. The three-electrode array was considered for 75 locations for each CRP (Zayandab Engineering Consulting Company 2010).
The inversion method was used to obtain geoelectrical mappings. Res2DInv (a finite element software) inversed the data (Geotomosoft 2014) based on the Gauss-Newton approach (Sasaki 1992), and the apparent resistivity data built a two-dimensional model of the study sites. The ground layers were divided into rectangular blocks with constant resistivity values. Block thickness was smaller near the ground surface as the influence of the block thickness on results was lower in terms of the root mean squared error (RMSE) (Loke 2002).
In the Schlumberger array, the thickness of the first layer of the model is equal to half of the electrode space, i.e., 25 m and 50 m for the Golab and Beheshtabad tunnels, respectively. The thickness of each layer increased by 10% compared to its top layer (Fig. 4c). RMSE values for all profiles were almost constant and lower than 7% after eight iterations (Fig. 4b). Figure 5 illustrates the geoelectrical mappings of the study sites.
In section 7.6-16 km of the Beheshtabad tunnel (from the south in 7.6 km toward the north in 16 km), Asmari limestone (brownish color) was detected in about 400 m, overlaid by Razak marlstone with low resistivity values (the green areas in Fig. 5a). The resistant zone (the brownish color in chainage 10 to 12 km) is sandstone and Bakhtiary conglomerate, up to deeper than 400 m. A fault in the interface of Razak marlstone and the conglomerate was visible (12 km). Toward the north, the Bakhtiary conglomerate is covered by low-resistant formations of marly limestone and mudstone between 11,500 and 15,500 m. ER varies due to fault activations and material contacts toward the outlet. The geoelectrical data were also compared with geological observations on the surface and geotechnical data of the boreholes (Zayandab Engineering Consulting Company 2010).
For the Golab tunnel, in high elevations of the first kilometer of the tunnel route (from the east toward west), ER was lower than the value of conglomerate near the tunnel elevation due to the activation of Fault Abrizan-Kamason and the presence of marly layers. The resistivity of siltstone and shale was low (shown by the blue color), and the resistivity of shale was typically lower than limy shale. In Chainage 7.5-8.5 km, ER decreased due to fault activations (Zayandab Engineering Consulting Company 2006).
Before utilizing geoelectrical mappings for predicting squeezing, the efficiency and accuracy of electrical data were evaluated. Influences of the electrode distance on geoelectrical outputs were measured by a sensitivity analysis. The effects of depth on the accuracy of electrical prospecting were also calculated. The multiple geophysical approaches were considered to cover the inability of the geoelectrical method to detect unrelated anomalies. Gravimetric and magneto-metric methods were performed for Golab and Beheshtabad tunnels to detect misleading anomalies, e.g., conductive minerals, caves, and water-bearing layers (Zayandab Engineering Consulting Company 2006, 2010. The sensitivity analysis was performed based on data collected from different maximum flow lines, i.e., 80, 140, 200, 300, 400, 600, and 800 m. We also used data with an electrode distance of 10 m to determine the effects of electrode distance, depth, and saturation on electrical resistivity. The geoelectrical data were recorded for two sections: 0-24 km of the Golab tunnel and 0-64.9 km of the Behehstabad tunnel (Zayandab Engineering Consulting Company 2006, 2010. The geological units in Figs. 3 and 5 were also detected at other parts of the tunnels (i.e., 8.8-17 km for the Golab tunnel and 0-7.6 km for the Beheshtabad tunnel). The range of ER for geological units was evaluated for the whole lengths of both tunnels to compare the electrical data based on water saturation, depth, and fracture properties. Rock mass classification indices, overburden depth, electrical resistivity, weathering, and water saturation conditions were recorded for all sections. An example of such an analysis is shown for a 2 m section ( Fig. 6) ( K sh,s 4 in 1550-2730 m of the Golab tunnel).
According to Fig. 6, the factors were sorted from unfavorable to favorable conditions for squeezing potential. The range of ER* (Table 1) was considered for the whole section, and ER value was estimated based on each factor. The average value was calculated for 2 m sections with slight weathering (Fig. 6), e.g., GSI = 57 (30 Ohm m), rock strength of 1.9 MPa (30 Ohm m), overburden 275 m (145 Ohm m), permeability 2 × 10 −8 m/s (30 Ohm m), which was impacted by joints (30 Ohm m), the average value of ER was calculated 68 Ohm m. We compared this estimated value with the ER profile for the section and calibrated ER * values. We Considered methodology to estimate the accuracy of ER values. An example of the results provided by this method is presented. All rock masses were saturated in the studied parts. A lower permeability can also reduce the squeezing potential prepared similar charts for all 2 m sections and geological formations along the tunnels.
This method assumes linear relationships between ER and the factors (Fig. 6). Some of these factors were qualitative parameters. During the design phase, many qualitative scales (such as rock mass classification scales) are usually used based on some assumptions (accepting some errors). We utilized this methodology to control if the ER value was reasonable for each section or if it was required to modify ER*.
The measurements based on a 10 m electrode distance were compared for the chainage 0-3050 m of the Golab tunnel and 8120-15,360 m of the Beheshtabad tunnel (Zayandab Engineering Consulting Company 2006, 2010. These data were also used to calibrate the ER* of the whole lengths of tunnels. We compared sections with almost similar factors (according to Fig. 6) at 10 m and 50 m electrode distances and modified ER* values (in case it was required). The Beheshtabad tunnel was not excavated, and we used the information on geological units at the surface and borehole data to evaluate each section. The final results of this study were based on the Golab tunnel; therefore, the accuracy of this method was based on the data in the production phase. i.e., we used the calibrated data of the Beheshtabad tunnel for predicting squeezing intensity, not for developing the squeezing predictive criterion.
The average difference between the initial and calibrated (calculated by a 10 m electrode distance and the methodology presented in Fig. 6) ER values was 19.6% for Beheshtabad and 11.35% for the Golab tunnel. For example, for shale ( K sh 4 ), the electrical values changed from 30-210 to 20-200 Ohm m, i.e., an average error of 11% was detected.
The calibrated values were then calculated by averaging electrical values for all geological units along the studied sites, and the initial electrical resistivity values (ER*) were modified to the calibrated values (ER) (Zayandab Engineering Consulting Company 2006, 2010).
The maximum difference was detected for the Bakhtiary conglomerate in Beheshtabad and the Conglomerate (E c ) in the Golab tunnel. These two units have the maximum overburden, higher than 400 m, which confirmed the effect of depth on the inaccuracy of electrical results. Table 1 summarizes the geoelectrical data before and after the sensitivity analysis.
The proton precession magnetometer (Scintrex Co.) with a sensibility of 0.1 nT recorded the magneto-metric data (nanoscopic alterations of the magnetic field). The nuclear magnetic resonance of the earth was checked, reaching a 500 m depth. This device was used to locate the possible existence of igneous masses, faults, conductive minerals, and other anomalies (Howland-rose 1981). CG-3 automated gravity meter evaluated the gravimetric variations (Scintrex limited 1995) along three profiles (with a reading resolution of 0.001 mGal) with a 50 m distance to detect sharp density contrasts, e.g., the presence of the water-bearing structures or cavities (Zayandab Engineering Consulting Company 2006, 2010, up to the depth of 500 m.
The prementioned complementary and calibrating methods increased the accuracy of ER data. Then, we used the data to predict squeezing potential at the design stage. Therefore, we answered the second challenging question of our predictive model: ER alterations can be considered a relatively accurate scale to predict squeezing after calibration. Hence, the remained challenge was defining a relationship between ER and squeezing.

Geotechnical studies
A series of laboratory experiments were conducted on prepared samples according to the ISRM standards (ISRM 1977a(ISRM , b, 1978 to determine geotechnical characteristics of the study sites, including porosity, density, uniaxial and triaxial compression strengths. The properties of rock joints (e.g., weathering and roughness) were also recorded (Zayandab Engineering Consulting Company 2006, 2010. The rock mass characteristics and classifications indicesincluding quality index (Q) and geological strength index (GSI)-were calculated based on the properties of intact rocks and rock joints. The geological properties of the study sites are shown in Tables 2 and 3.
The rock mass deformation modulus was calculated using standard empirical equations (Khanlari et al. 2012;Shen et al. 2012). The effect of disturbance on rock mass properties for the sections of the Golab tunnel excavated by drilling and blasting was also considered (Sections 3 and 8) (Rockscience 2007). The rock mass deformation modulus was calculated by averaging the values estimated by different empirical equations to decrease the error, because each empirical equation depends on some rock mass properties (Khanlari et al. 2012).
Considering the intact rock properties as an upper limit and using the standard deviation values, the estimated values of rock mass characteristics were obtained close to actual values. The properties of intact rock and rock joints and the in-field observations were used to calculate the rock mass parameters, e.g., cohesion strength and internal friction angle (using RocLab software) (Rockscience 2007).
The study sites (the 0-8.8 km chainage of the Golab tunnel and the 7.6-16 km chainage of the Beheshtabad tunnel) were divided into sections based on lithology, structural features, overburden, and geotechnical characteristics. Tables 4 and 5 summarize the range and average value for all parameters.

Squeezing prediction
In this section, the observed intensity of squeezing along the Golab tunnel during excavation was related to ER values. A relationship was defined between ER and the measured normalized strains (based on in situ data collected by extensometers). Finally, a new strain-based classification system was developed. Throughout this text, N, L, M, S, and VS depict non-squeezing, light, moderate, severe, and very severe squeezing.

Combining numerical modeling and semi-empirical approach
Incorporating simple numerical modeling into a semi-empirical approach was previously considered an option to predict tunnel squeezing (Khanlari et al. 2012). The possibility of instability and support problems can be higher when tangential strain at the tunnel periphery passes the critical threshold (Aydan et al. 1993). The Aydan method (Eq. 1) can be linked to numerical modeling to calculate rock deformations around a tunnel. The squeezing index can be calculated as the ratio of calculated to critical strain: For SI < 1, tunnel squeezing is impossible; for 1 < SI < 2, L; for 2 < SI < 3, M; for 3 < SI < 5, S; and SI > 5, VS (Aydan et al. 1996). The critical strain can be calculated by experiment-based relations, such as Eq. 2 (Singh et al. 1997) and Eq. 3 (Barton 2002): Numerical modeling was used to evaluate the normalized radial closure, u a /a; the finite difference program FLAC 2D (ITASCA 2002) was used to calculate strain values. Input parameters were residual stress, dilation angle, and residual strength properties of the rock masses. Because the overburden of tunnels was high, the dilation angle can be considered zero (Mair et al. 2002).
With an assumption of the ideal elastic-plastic rock behavior, the residual strength parameters were set similar to the peak values. Therefore, the rock mass characteristics, including the deformation modulus, poison ratio, cohesion, and internal friction angle (i.e., the outputs of RocLab), were incorporated in FLAC 2D . For the stress state, vertical in situ stress (γz) was applied. The horizontal stresses were (1) SI = u a ∕a cr .
(2) cr = 31.1 1.6 ci E i 0.6 Q 0.2 , (3) cr = 5.84 0.88 ci Q 0.12 E 0.63 i calculated by multiplying K values (the ratio of horizontal to vertical stress, i.e., Eqs. 4 and 5) by γz. Maximum and minimum stress coefficients were determined by the following equations (Khanlari et al. 2012): In the numerical models, the horizontal fixed boundary condition was considered at the left and right sides of the domain, and the XY-fixed boundary was imposed on the bottom model boundary. Several models were simulated to evaluate the optimum distance from the excavation boundary to avoid boundary effects. This distance was calculated as ten times the tunnel diameter.
Each model side was considered 100 m based on several runs of numerical modeling. In each direction, 100 grids were set; a finer mesh size was used close to the tunnel periphery. After reaching the initial equilibrium (i.e., in situ stress condition), the excavation was simulated. Then, the support system was installed. The support system was installed at a distance from the excavation face (by permitting a certain amount of deformation before installing the support) (Rocscience 2004).
The evaluated closure values were implemented in the semi-empirical equation (Aydan et al. 1996), and the squeezing was predicted based on the calculated SI values. This method (the Aydan-FDM) combined simplified numerical modeling and semi-empirical methods (Fig. 7). We compared the results of this method to real squeezing cases along the Golab tunnel to determine its accuracy compared to the production-phase data.

Electrical resistivity and observations
A comprehensive data set was available for time-dependent deformations, support failures, financial costs, and rock mass conditions for the excavated Golab tunnel. The rock and stability conditions of the Golab tunnel were compared with the frequently used qualitative classification of tunnel squeezing (Aydan et al. 1993;Palmstrom 1996;Hoek and Marinos 2000), and the squeezing intensity was defined along the tunnel length. The main drawback of predicting squeezing with this approach was the qualitative presentation of the squeezing. Previous semi-empirical equations were also developed based on site-specific data and some, not all, geotechnical properties.
The normalized tunnel strain (ε) can be considered the only quantitative parameter to classify squeezing. However, usage of this parameter can be problematic: (1) calculating 150 z Fig. 7 Flow chart of the Aydan-FDM for predicting squeezing ε based on empirical equations or numerical models, e.g., (u a /a) in the Aydan-FDM introduced in "Combining numerical modeling and semi-empirical approach", can cause an error.
(2) For accurate measurement of ε, we required the recorded deformations during or after the excavation. However, the objective of this research work (similar to many other projects) was to predict squeezing in the design phase.
Correlating ε and another parameter with two intrinsic properties was the aim of this work. Such a parameter should be related to tunnel squeezing and can be measured before excavating the tunnel. ER was chosen as the parameter satisfying these two requirements, which can indirectly predict time-dependent behavior, based on what was mentioned in the Introduction section.
Calibrated ER values were related to the qualitative squeezing classification along the Golab tunnel. Then, an ER range was considered for each squeezing intensity class. Since the in situ closure data of the Golab tunnel were recorded during the excavation by a convergence tape (Moosavi and Khazaei 2003), accurate values of ε were measured at 60 sections. In each part, five pins were installed to measure the displacements with an accuracy of 0.01 mm (Zayandab Engineering Consulting Company 2006). Therefore, we related ER to the in situ values of ε. The ε-ER chart was prepared based on the recorded in situ data. Finally, the effects of various parameters on squeezing intensity were evaluated using the newly developed approach. The squeezing prediction based on the ER is summarized in Fig. 8.
The steps in Fig. 8, except the last one (the strain-based classification), are related to the data of the Golab Tunnel used for calibration. Then, the normalized strain-squeezing intensity was used to predict the squeezing of the Beheshtabad Tunnel. The quantitative approach was used to predict squeezing for various conditions of the rock type, depth, saturation, and joint properties.

Results
The in situ stress state of the study sites was evaluated by Eqs. 4 and 5 (Tables 6, 7). The strength parameters of some rock masses were lower (Tables 2, 3, 4) compared to the Golab tunnel; therefore, the squeezing problem can be more plausible along the Beheshtabad tunnel.
We considered the average values of Eqs. 2 and 3 as the critical strain in each section (Table 8). Each value was evaluated as a range, because we defined a range of intact rock properties based on laboratory experiments. The values of u a a (calculated by FLAC 2D ) were divided by the average value of ε cr (Table 8). Figure 9 shows displacement contours for Section 8 of the Golab tunnel in a 20 m × 20 m plot window. The values in Table 8 were incorporated into Aydan et al. (1996) semi-empirical equation to predict the squeezing intensity (the Aydan-FDM). Table 8 lists the in situ measured strains. A ε range was also presented, because several measurements were recorded along each section of the Golab tunnel. The difference between in situ ε and u a ∕a can illustrate the accuracy of the numerical modeling, i.e., the Aydan-FDM.
The qualitative classification based on previous studies (Aydan et al. 1993;Palmstrom 1996;Hoek and Marinos 2000) is shown in Table 9. This classification was used to determine the squeezing intensity along the Golab tunnel, and the range of ER values of the tunnel was recorded for each intensity class. Table 10 shows calibrated electrical values for each section and the variation of squeezing conditions along the Golab tunnel. The squeezing intensity based on the ER and  Aydan-modified approaches is also presented. At the GSI range of 60-70, the ER of saturated siltstone (the green bar) and conglomerate (the red bar) was decreased from 150 to 60 and 220 to 190 Ohm m, respectively (Fig. 10a). The decrease of ER due to the water saturation, under constant GSI (Fig. 10a), was higher for the rocks with higher amounts of clay, e.g., 60% and 41% for siltstone and shale, while 13.6% and 20% for the conglomerate and limy shale. At approximately constant depth, the effects of rock mass condition (or GSI) on the ER of the saturated tocks (Fig. 10b) were also noticeable. The resistivity decreases from 23.5 (for conglomerate) to 62.5% (for limy shale).   Figure 11 shows the occurrence of stability issues along some sections of the Golab tunnel. The squeezing problems based on the ER and Aydan-modified are presented in Table 11. The observed squeezing cases and ER variations are shown (Fig. 12) from 0 to 3.25 km along the Golab tunnel.
Based on the field observations, the length of N, L, M, S, and VS squeezing conditions was 1.105, 5.305, 1.745, 0.74, and 0.375 km, respectively. The squeezing calculated based on the Aydan-FDM was evaluated for each squeezing class. The total length of sections predicted with N, L, M, S, and VS conditions was 6.03, 0.34, 0, 0, and 2.43 km, respectively. Therefore, the Aydan-FDM predicted the squeezing with an average error larger than 100%.
The prediction accuracy of the Aydan-FDM may vary for other formations as it was developed based on the accuracy of geotechnical parameters at the study site. For the ER-based approach, some sections were divided into several squeezing classes, e.g., for Section 2 of the Golab tunnel, chainages 1210-1319, 1319-1495, and 1495-1550 m were detected as no, light, and moderate squeezing classes, respectively. According to the foreseen lengths for intensity, the accuracy of geoelectrical prospecting along the Golab tunnel was at least 94%. Therefore, the qualitative classification based on geoelectrical mapping was more accurate than the Aydan-FDM, and squeezing can be evaluated accurately in a section, not a location. Hence, we developed the relationship between the in situ measured strains and ER and prepared the new classification system based on ε. This methodology provided a quantitative relationship for predicting squeezing (Eq. 6). Figure 13a presents a chart based on the relationship of the in situ measured strains (Table 8) and ER, and a new strain-based classification system was introduced for predicting ε (and consequently squeezing) based on ER (Eq. 6). In the next stage, for each squeezing class, the ER-ε relationship was used to calculate the strain range based on ER (Table 12 and Fig. 13b). The stepwise shape of Fig. 13b was unrelated to the correlation between data, and it depicts the range of predicted normalized strain for each squeezing class in the collected data of this study.
This new approach (Eq. 6) provided quantitative results for predicting squeezing based on the strains calculated by ER values. The accuracy of the new classification system was 89.14 %, considering the R-square value of the ER-ε chart and the prediction error of the resistivity approach, 94%.
The squeezing potential along different sections of the Beheshtabad tunnel was predicted (Table 9)-presented in Fig. 14a. The squeezing potential can vary even along a geotechnical section. For instance, the squeezing changed from M to VS class along Section 2. A similar trend was observed during the excavation of the Golab tunnel. Variations of squeezing intensity can be related to (6) (%) = 3538 ER −1.523 30 < ER (Ohm m) < 270.  Fig. 14a shows the range of ER along a section of the Beheshtabad tunnel and presents the predicted squeezing intensity based on the geoelectrical qualitative approach. For example, the squeezing was defined as S from 7600 to 8400 m based on two ranges of ER (60-90 and 75-90 Ohm m).
For the comparison analysis, the predictive values of the Aydan-FDM are also presented in Fig. 14b. The predicted length of M, S, and VS squeezing classes based on ER approach along the Beheshtabad tunnel was 3425, 2355, and 2620 m, respectively. No zone was predicted as S or L squeezing conditions. The squeezing intensity along the Beheshtabad study site was also evaluated by Eq. 6 (Fig. 14c). Variations in Fig. 14a, c was related to the fact that a range of ER or strain was recorded for each squeezing condition along tunnel sections. For example, for the first section (7800-8400 m), ER and ε ranged from 60 to 90 Ohm m and 3-7%, respectively. However, the squeezing condition was predicted as severe for the whole range.
The stepwise format was related to variations of squeezing intensity along sections of the tunnel and unrelated to the accuracy of the methodology. For example (according to Fig. 14), in chainage 9200-10,000 m, ER and ε varied between 45 and 85 Ohm m and 2.7-10.5%, respectively. The squeezing intensity rose dramatically from M to VS conditions, with some variations of S and VS conditions between 9500 to 10,000 m.
The relation of ER with rock mass indices and rock mass properties along the Golab tunnel is represented in Fig. 15. We evaluated the effects of these parameters on ER using the newly developed quantitative strained-based system (Fig. 16). The higher the Q and GSI values were, the bigger the ER value was. Therefore, a direct relationship was confirmed between ER and rock mass condition (as predicted).
In Fig. 16, the effects of water saturation and GSI is shown. The depth was kept almost constant (depth difference < 30 m) while other parameters were compared (without depth effects). Then, the ER values recorded in these conditions were used to predict the normalized strain.  The strains increased by about 151% and 187% for mudstone and marlstone, respectively (for a constant GSI value and a depth difference of less than 30 m). The increase was 109% and 114% for the conglomerate and limestone. Therefore, the squeezing intensity changed: from L to M for conglomerate and limestone, from M to S for mudstone, and from M to VS for marlstone (see Table 12 and Fig. 16a).
At an approximate constant depth, the effects of GSI on ER-based strains of the saturated rocks were also noticeable. The measured ε of mudstone, marlstone, conglomerate, and limestone were increased by about 160%, 187%, 107%, and 167%, respectively, i.e., the squeezing intensity altered from severe to VS, M to VS, M to S, and M to VS. The maximum strain of saturated marlstone decreased from 10.7% to 6.9% for a depth increase of 300 to 360 m, while the maximum ε of saturated mudstone remained constant at 12.8% for depths 330 and 380 m.

Discussion
The ER-based approach can be a valid hypothesis to predict squeezing at the design phase. Meanwhile, technical measures at the production phase (e.g., installing heavy supports) should be evaluated based on complementary data achieved during the excavation, e.g., coupled three-dimensional numerical models (using more accurate geotechnical data, e.g., in situ measurements). We indirectly calculated the time-dependent rock mass deformations similar to the conventional approaches of predicting squeezing. The difference was related to using ER instead of classification systems or empirical equations. This assumption can result in more accurate outcomes for the design phase.
The advantages of predicting squeezing based on ER can be (1) geophysical methods are fast, non-destructive, and environmentally friendly. (2) Geophysical setups, e.g., the number of electrodes, are adaptable to cover the whole study area. Therefore, the final results are based on large-scale data compared to the limited geotechnical data gathered from boreholes/galleries. (3) The parameters affecting squeezing, e.g., the hydraulic and geologic characteristics, can influence ER likewise. Including more influential parameters (e.g., joint condition, degree of weathering) can increase the accuracy. (4) Geotechnical sections can be divided into more sub-sections with different squeezing classes, and the accuracy can be higher. The number of subsections for the Aydan-FDM and ER-based methods (Fig. 14a, b) was 14 and 6, respectively. (5) Geophysical studies (particularly in large projects) are common and without imposing extra financial costs to the design phase, like performing in situ tests.
All available geotechnical approaches present a universal classification system for all types of squeezing (Aydan et al. 1993;Palmström 1996;Hoek and Marinos 2000), including those related to the topographic or tectonic mechanisms. This classifying approach may result in erroneous results. However, the developed strain-based equation can provide more accurate data for squeezing in soft rocks (using the calibrated ER values).
Many squeezing cases are tectonic-based and can occur at a low depth (Shrestha and Broch 2008). The predicted squeezing intensity (based on measured strains) was detected unrelated to depth. Some previous works categorized the squeezing only on the overburden (Singh et al. 1992;Goel et al. 1995). The probable reason for this unclear relationship can be the complicated stress state of the study site and factors (such as the stress ratio, tectonic conditions, and fault activations). Tectonic conditions are partially understood and may alter based on findings during the production phase (Guan et al. 2012). Therefore, it can be impossible to transit these in-field characteristics into a quantified scale, and most previous studies did not consider a tectonic parameter to predict squeezing (Sun et al. 2018). However, detecting fracture zones and faults is one of the main applications of geophysical studies, and tectonic-based squeezing can be predicted at the design phase based on ER variations. Dramatic decrease of ER due to activation of faults showed that the ER-based system predicted tectonic-based squeezing; this ability was visible for the Golab tunnel, e.g., for 0-500 m and 7500-8100 m sections.
A crucial question must be answered: why a significant difference was recorded between the Aydan-FDM and geoelectrical predictive approaches? Several explanations can be presented: (1) the efficiency of the implemented geotechnical data employed in the numerical tool was effective. Determining the rock mass properties based on specific core samples gathered from a few boreholes can  cause erroneous results. This uncertainty was crucial in deep tunnels, where borehole drilling was challenging (and costly) and could not cover large parts of the study area. For example, all geotechnical data of the 7.6-16 km chainage of the Beheshtabad tunnel were evaluated only by four deep boreholes. Therefore, the accuracy of geotechnical input data can impact the numerical model. The considerable difference between calculated u a a and measured ε proved the inaccuracy of the numerical results. However, detailed geotechnical studies (e.g., drilling more boreholes or conducting in situ tests) may provide more accurate geotechnical data and less erroneous squeezing predictions by numerical tools.
(2) The impacts of the geotechnical data on the semiempirical equations were vital. Improper input parameters can result in poor predictions of the Aydan-FDM (i.e., the ε cr equations were also based on the geotechnical parameters of intact rocks, i.e., limited laboratory experiments).  (3) The intentionally simple numerical scheme was used to calculate deformations by FLAC 2D because of the limited geotechnical data in the design phase. Some factors (e.g., thermal data or chemical parameters) are rarely measured in the design stage. A fully coupled numerical model can improve the results, but running a coupled model specific to the study site is time-consuming. Several three-dimensional models should be built for long tunnels due to the parameter variations. Therefore, in the design phase of many projects, using numerical schemes is limited to simple models. (4) The Aydan equation was developed based on geotechnical data from tunneling projects in Japan; different site-specific conditions may result in different outputs of this approach.
Despite being more accurate, the main drawback of the ER approach was presenting squeezing based on a qualitative explanation (Table 9) instead of quantitative data. Table 9 shows some squeezing classes with shared values of ER, e.g., 80 Ohm m was introduced as a moderate or severe squeezing class. The electrical approach developed in this research cannot also predict the squeezing intensity of hard rocks in deep tunnels. A reason for this incapability was the chosen electrical range (30-270 Ohm m) in this study; this range was suitable for soft sedimentary rocks. Another limitation was using an almost similar support system for the studied tunnels. Using a heavier support system can affect the squeezing intensity significantly.
Considering the effects of different parameters on the squeezing, alterations of the intensity by mechanical stress (σ cm ) were almost invisible (despite increasing). Overburden did not impact squeezing (alterations did not follow a unique trend under depth variations at a constant GSI). The reason can be other factors (such as the stress ratio (K), tectonic parameters, and fracture properties).
The squeezing conditions along the Beheshtabad tunnel showed that most parts of the tunnel are prone to severe squeezing, and employing mechanized excavation can be problematic in these zones. Different geological formations in the tunnels, including shale, conglomerate, limestone, siltstone, marlstone, mudstone, and sandstone, were studied for squeezing problems. For both tunnels, high-intensity squeezing was mainly recorded for formations with a high amount of clay, e.g., siltstone and shale at the Golab tunnel. A high decrease of strains (based on ER) due to water saturation was visible for the rock masses with high content of clay, i.e., mudstone and marlstone, in the Beheshtabad tunnel.
In the case of complementary studies of predicting squeezing: (1) a more comprehensive database of measured strains (especially from other tunnels excavated in soft sedimentary rocks) can increase the R-square value of the ER-ε chart and the accuracy of the predictive approach. (2) Other exploration approaches (such as integrated seismic detection) can be combined with geoelectrical data to evaluate tunnel squeezing.
Due to active tectonic structures and high overburden, the stress condition of the Sanandaj-Sirjan zone in Iran is crucial. Regarding the enormous executed/planned tunneling projects in the Central Plateau of Iran, e.g., the Rokh valley water-supply tunnel and Kohrang tunnels, geoelectrical data can predict the squeezing for these projects at the design phase. In the case of deep tunnels excavated in soft rocks, the developed approach predicted the squeezing more efficiently comparing classical methods, although complementary investigations are required.

Conclusions
Predicting squeezing at the design stage is a challenging task. Accurate predictive methods can prevent technical difficulties, financial losses, and human injuries during or after the production stage. This study presents the ER-based approach to predict squeezing potential along a tunnel. The following conclusions can be highlighted: • The Aydan-FDM cannot accurately predict rock deformation in soft sedimentary rocks of the Sanandaj-Sirjan zone, maybe due to improper geotechnical parameters or specific conditions of the study sites. The extreme underestimation of squeezing is visible in all sections of the Golab and Beheshtabad tunnels. • The efficiency of electrical resistivity to predict the squeezing potential is noticeable for soft sedimentary rocks. • Having a high accuracy in predicting squeezing, the ability to cover a large study area, and being environmentally friendly are all advantages of the ER-based method over conventional approaches. • ER is related to the properties of intact rock, fractures, and rock mass. We predict squeezing considering the effects of water saturation, rock type, and joint condition (based on the geoelectrical approach). • We relate ER to in situ measured strains (a post-excavation parameter). Then, we develop a quantitative strainbased squeezing classification system for tunnels excavated in soft sedimentary rocks. • The accurate prediction of tectonic-based squeezing is another advantage of the newly developed method. • Effects of the rock type, water saturation, and GSI on squeezing intensity are noticeable.