Three-Dimensional Numerical Investigation of the Interaction Between Twin Tunnels

The construction of twin tunnels is a mandatory guideline and a prevailing practice in either conventional or mechanized tunneling. Nevertheless, most of the design methods for calculating the tunnel loads focus on single tunnels, neglecting thus the potential interaction between neighboring tunnels. The effect of such interaction can be significant, especially for closely-spaced twin tunnels. In this context, this paper investigates via parametric 3D finite element analyses the interaction between deep, parallel-twin, circular and non-circular tunnels excavated with a conventional (i.e., non-TBM) method and supported with shotcrete lining. The numerical investigation focuses on the axial forces acting on the primary support of the tunnels by examining the effect of a wide range of geometrical (pillar width, overburden height, tunnel diameter and section (shape), lagging distance), geotechnical (strength and deformability of the surrounding rockmass, horizontal stress ratio), structural (thickness and deformability of the shotcrete lining) and construction parameters (full- or partial- face excavation and support of the tunnels). The results of the analyses indicate that the construction of the subsequent tunnel influences the loads of the precedent. The stress state of the single tunnel is used as the reference for the quantification of the interaction effect. The output is presented in normalized design charts of the quantified interaction effect on the axial forces versus key geomaterial and geometry parameters to facilitate preliminary estimations of primary support requirements for twin tunnels. Furthermore, nomographs are provided for preliminary assessments of the optimum pillar width (spacing) between twin tunnels, which practically eliminates the interaction effect.


Introduction
It is common engineering practice to construct deep, parallel-twin tunnels with conventional methods (Sequential Excavation Method, Observational Method (Kovári and Lunardi 2018), etc.) within weak rockmasses in highway and railway networks. Cases of twin tunnels constructed within weak rockmasses and associated with challenging or ''squeezing'' conditions for tunneling construction, stress instability problems and potentially excessive radial displacements (convergence) are, among others, the Kallidromo Railway Tunnel (Diasakos et al. 2010) and some of the Egnatia Highway Tunnels . The twin tunnels are excavated and supported with a lag; thus, the second branch is typically constructed after the first has advanced sufficiently to maintain a longitudinal separation distance between the faces of the tunnels. Therefore, the advance of the subsequent tunnel mobilizes stress and strain redistribution in the area (pillar) between the tunnels, resulting in the additional loading of the precedent tunnel . However, most design methodologies calculate the loads for single tunnels , neglecting thus the potential interaction between closely-spaced twin tunnels, more accurately, for twin tunnels with pillar width less than or up to two tunnel diameters (i.e., W B 2D), indicatively according to Ghaboussi and Ranken (1977), Soliman et al. (1993), Kim et al. (1998), Elwood and Martin (2016), Fang et al. (2016), Shaofeng et al. (2018). The literature review on the under-study problem (presented in the next section) indicates that significant interaction phenomena occur between adjacent tunnels and further identifies the crucial effect of the pillar width on controlling their intensity Vlachopoulos et al. 2018). Most of the studies use 2D analyses, which inherently involve simplifications and have practical limitations (Vlachopoulos and Diederichs 2014) and cannot realistically capture the purely 3D nature of tunneling, especially if ''squeezing'' conditions prevail for tunneling construction (Cantieni and Anagnostou 2009). On the other hand, studies that utilize 3D analyses to model twin tunneling are either limited to a relatively narrow range of variable parameters or targeted to specific projects (i.e., case studies), since such simulations are computationally demanding. Moreover, most of the relevant research focuses on shallow, either conventional or mechanized, twin tunneling interaction.
Acknowledging the latter points, this paper includes an extensive series of parametric 3D-FE analyses to identify and quantify the interaction effect between deep, parallel-twin tunnels, constructed with conventional methods within weak rockmasses, by varying the following parameters: pillar width (W), overburden height (H), tunnel diameter (D) and section (shape), lagging distance (L), strength (r cm ) and deformability (E m ) of the surrounding rockmass, horizontal stress ratio (K o ), thickness (t sh ) and deformability (E sh ) of the shotcrete lining and construction sequence (full-or partial-face) of the tunnels. The results of the analyses are used to produce general-purpose normalized graphs (design charts) and propose analytical equations for the axial forces developing in the primary support of twin tunnels versus key geomaterial and geometry parameters. These charts and equations can be useful for preliminary estimations of primary support requirements for twin tunnels. Moreover, nomographs for the optimum pillar width (spacing) between twin tunnels, which practically obliterates the interaction effect, are provided. These nomographs consider the influence of the variation of multiple tunneling parameters and constitute a practical guide to limit or control the interaction effect to negligible or non-critical levels, therefore, ensuring the stability of twin tunnels.

Literature Review
The literature review indicates that the interaction effects between twin tunnels are investigated via processing monitoring data (field-measurements) recorded during their construction phase, using analytical and empirical methods, performing physical modeling tests, conducting numerical analyses and carrying out back analysis of failures occurring during their excavation and support sequence. The main points of the relevant literature are presented in the following paragraphs. Ghaboussi and Ranken (1977) summarize the main conclusions of relevant previous research about twin tunneling interaction and conduct 2D-FE analyses, identifying that for pillar width more than two tunnel diameters (W/D C 2), the interaction is eliminated. Leca (1989) studies the interaction between twin tunnels constructed with conventional or mechanized methods in soft ground (clayey soil) relative to the tunnel face stability, the design of the support system and the ground surface settlements. Adachi et al. (1993) perform model tests to investigate the behavior of shallow, circular, twin tunnels in sand and quantify the interaction via the overburden over the pillar width ratio. Respectively, Kim et al. (1998) perform model tests to examine the response of shallow and deep, circular, closely-spaced tunnels in clay, concluding that for W/D C 1.5, the interaction is obliterated. Chapman et al. (2007) perform experiments with physical models of deep, circular, multiple (side-by-side) tunnels in clay and propose a modified gaussian method to calculate the ground surface settlements. Choi and Lee (2010) conduct experiments with physical models for shallow, circular, twin tunnels in a material simulating a rockmass with low strength, studying the evolution of the displacements and the propagation of the cracks in the material, suggesting that W [ 3D is required to control the interaction. La et al. (2018) study the interaction between deep, non-circular, twin tunnels with different dimensions in a material simulating a rockmass with low strength via both physical modeling tests and numerical simulations by assessing the distributions of the principal stresses and the strength/ stress ratios of the material for varying geometrical configurations (spacings and diverging directions) of the tunnels. Jiang et al. (2021) study the overall failure process of deep, circular, twin tunnels (located in the same or different depth) via physical modeling simulations based on the 3D printing of sandstone analogs models and present a safety factor method for evaluating the stability of the multi-tunnel system, based on the propagation of the cracks in the material. The results of the experiments for the failure process of the neighboring tunnels are verified with numerical back simulations. Potts (1996, 2001) utilize 2D-FE analyses to study the response of shallow, circular, neighboring tunnels (either side-by-side or one above/below the other) in London clay, indicating that the relative position and the spacing of the tunnels significantly modify the profile of the ground surface settlements due to the interaction phenomena, which become insignificant for W [ 7D. Koungelis and Augarde (2004) carry out 2D-FE analyses for the same parameters and conditions, concluding that for W [ 4D, the interaction is drastically reduced. Chapman et al. (2002) use 2D-FE analyses to identify the influence of shallow, circular, multiple closely-spaced tunnels in clay on the evolution of the ground surface settlements. Suwansawat and Einstein (2007) and Yang and Wang (2011) combine the superposition technique with empirical and analytical methods to acquire the accumulated ground surface settlement trough for twin tunnels. Ng et al. (2004) investigate via 3D-FE analyses the influence of the lagging distance between the faces of shallow, non-circular, twin tunnels, conventionally constructed, and discuss the load transferring mechanism between the tunnels focusing on the ground surface settlements, the radial displacements and the internal forces developing in the linings of the tunnels. Similar research with 3D-FE analyses targeted to twin mechanized tunneling interaction in soft ground (clayey soil) is presented by Do et al. (2014Do et al. ( , 2015Do et al. ( , 2016. Furthermore, among others, Yamaguchi et al. (1998), Hage Chehade andShahrour (2008) and Chakeri et al. (2011) carry out 2D or 3D-FE analyses for analogous parameters and conditions focusing on the effect of the ground, the size, depth and relative position of the tunnels on the ground surface settlements, highlighting that the influence of the new tunnel on the existing decisively depends on the relative position and the spacing of the tunnels. Chang et al. (1996), based on a case study of twin tunnels in Taiwan, propose a simplified method to estimate the safety factor for the stability of the pillar. Chen et al. (2009a, b) investigate the effects of the rockmass pillar width on the behavior of three and four parallel-branch, deep, circular and non-circular tunnels for a case study in Taiwan, via 2D and 3D-FE analyses, defining W/D = 2-4, as the critical pillar width. Liu et al. (2008) examine via 3D-FE analyses the interaction between shallow, non-circular, conventionally constructed tunnels in rockmass for a case study in Sydney, pointing out that the advance of the new tunnel significantly affects the existing, with the extent of the influence depending on the relative distance between the tunnels. Vlachopoulos et al. (2013) and Connor Langford et al. (2016) present a case study including excessive deformations and primary support failures encountered during the construction of the two-branch Driskos Tunnel in the Egnatia Highway, associated with weak rockmasses. Detailed insights into how tunnels designed in such geological and geotechnical conditions can be realistically analyzed by employing sophisticated 2D and 3D-FE models and evaluating field measurements are provided and the optimization of the primary support via a quantitative risk approach is discussed. Respective intense interaction phenomena between twin tunnels, also related to weak rockmasses, resulting in uncontrolled deformations and severe primary support failures, are also recorded among others in the Kallidromo Tunnel, as presented by Diasakos et al. (2010). Based on case studies of the construction of shallow or deep twin tunnels, the interaction effects are evaluated by other research via conducting 2D or 3D-FE analyses. For instance, Shaofeng et al. (2018) discuss the stability of twin tunnels according to the principles of the convergence-confinement method, Ghorbani and Ajalloeian (2019) present a method to define the optimum distance between twin tunnels, Fang et al. (2016) and Das et al. (2017) study the ground surface settlements for twin tunnels with different geometrical arrangements, Elwood and Martin (2016) deal with the response of the surrounding geomaterial for twin tunnels constructed in heavily overconsolidated clayey soil. Siahmansouri et al. (2012) propose, with the aid of 2D-FE analyses, a method to calculate the required rock pillar width, with criterion its stability, between deep, circular, twin tunnels. Fortsakis et al. (2012) condense the main results of relevant previous research about twin tunneling interaction and conduct 3D-FE analyses for marginally shallow, circular, twin tunnels conventionally constructed in weak rockmasses, concluding that for W/D C 2, the interaction is negligible and that the effect mainly concentrates on the precedent tunnel than the subsequent. The average additional load, developing in the primary support, varies from 10 to 60 % for the existing tunnel, depending on the pillar width and the geotechnical conditions, while it is constant and oriented at practically negligible levels, up to 10 %, for the new tunnel, despite the pillar width and the geotechnical conditions. Kumar (2013, 2018) use the upperbound (UB) finite element linear analysis (FELA) to study the stability of circular, twin tunnels in cohesive and cohesive-frictional soil. The required internal support pressure (s i ) to maintain the stability of the tunnels is calculated as a function of the dimensionless variables S/D, H/D, cD/c, u (S = net spacing between the tunnels, H = tunnel cover, D = tunnel diameter, c = cohesion, u = friction angle). The optimum net spacing (S opt ) between the tunnels required to eliminate the interaction effect increases versus increasing H/D and decreasing u. The required internal tunnel support (s i ) increases versus increasing H/D and cD/c or decreasing S/D and u. The S opt lies approximately in a range of 1.5D-3.5D for H/D = 1 and 7D-12D for H/D = 7. Yamamoto et al. (2013) and Yamamoto et al. (2014) present analogous research, using the same analysis method, for the stability of circular and rectangular, twin tunnels in cohesive-frictional soil subjected to surcharge loading, while Zhang et al. (2016), Yang et al. (2017) and Zhang et al. (2018) for the stability of circular, elliptical and horseshoe, twin tunnels in cohesive-frictional soil subjected to gravity, respectively. In this context, Rahaman and Kumar (2020) utilize the lower-bound (LB) finite element linear analysis (FELA) to study the stability of horseshoe, twin tunnels in rock subjected to surcharge loading. The maximum permissible surcharge pressure on the ground surface, which does not compromise the stability of the tunnels, is calculated. The loading factor q per /r ci depends on the following dimensionless parameters H/B, S/B, r ci /cB, D, GSI and m i . The effect of the net spacing (S) between the tunnels is also evaluated. The calculations indicate that the optimum net spacing (S opt ), beyond which the tunnels respond as quasi-single tunnels, is almost independent of the properties of the rock and increases continuously versus increasing H/B. The S opt varies from 2.5B-3B for H/B = 1 to 13B-14B for H/B = 6. Zhang et al. (2019) and Xiao et al. (2019) present analogous research, utilizing the same analysis methodology, for the stability of circular and rectangular, twin tunnels in rock subjected to surcharge loading. In a relatively similar research, Guo et al. (2021) propose an elastic analytical solution for the stresses developing in the geomaterial surrounding circular, twin tunnels in isotropic stress field and validate it with numerical simulations.
Irrespective of the type of the conducted research, it is derived that the pillar width between the tunnels is a factor governing the level of the interaction and, consequently, an essential parameter in controlling it. The majority of the research mainly focuses on the radial displacements (convergence) developing in the sections of the tunnels, the stress and strain alterations occurring in the surrounding geomaterial of the tunnels (for deep and shallow tunnels) and the ground surface settlements (for shallow tunnels). On the other hand, the minority of the research studies the internal forces acting on the primary support or the final lining of the tunnels. Therefore, based on the investigation areas of the previously mentioned research, this paper aims to enrich the available literature and the tunnel engineering arsenal with a computing framework for the interaction between deep, circular and non-circular, conventionally constructed twin tunnels, focusing on the loads imposed on the primary support of the tunnels. In addition, it provides insights on the required spacing between twin tunnels that eliminates the load alteration due to the interaction, regarding as reference the loads of the single tunnel.

Numerical Models and Analyses
The present investigation includes parametric 3D numerical analyses with the Simulia Abaqus 6.14 Finite Element Code (Abaqus Documentation 2014). The tunnels are assumed to be either circular with a diameter equal to D = 8 and 12 m or non-circular with an equivalent diameter equal to D eq = 7.6 m. The excavation length of the tunnels is L exc = 10D. Figure 1 presents 3D numerical models with pillar width to tunnel diameter ratio (henceforward referred to as pillar width ratio for brevity) W/D = 2 for circular and non-circular tunnels. Except for the pillar width, which is a variable parameter, the basic dimensions (excavation lengths, limits of external boundaries) of the models are constant regardless of the tunnel diameter (D). The basic dimensions of the models are increased versus those assumed by Fortsakis (2012) and Fortsakis et al. (2012) and either comparable or increased versus those adopted by Vlachopoulos et al. (2013) and Vlachopoulos et al. (2018), in similar 3D numerical simulations that calculate the loads of single and twin tunnels, to minimize the boundary effects. Benchmark checking, not presented herein, has shown that an additional increase of the limits of the external boundaries of the models by 2D-4D has a computational/time cost, while it has a non-critical influence on the tunnel loads since their differentiation is up to 1-2 %. Also, sensitivity analysis of either coarser or finer meshes has oriented the deviations in the tunnel loads to noncritical levels (± 1-2 %). The overburden height (H), which is also a variable parameter, is simulated by applying an additional uniform vertical load on the surface of the models. The models use solid elements for the rockmass and shell elements for the primary support (shotcrete).
The geomaterials used in the analyses with uniaxial compressive strength (UCS) of the intact rock varying from 10 to 30 MPa, according to the classification of rocks as presented by Brown (1981), correspond to soft rock, weak fractured or heavily fractured rock and marginally medium-strong rock. For the surrounding rockmass, the geological strength index (GSI) (Marinos and Hoek 2018) varies from 10 to 40 and the ''global'' UCS (r cm ) (Hoek and Brown 2019) varies from 0.5 to 5.5 MPa. Examples of such geomaterials can be the following: sedimentary soft rocks, mainly silty/argillaceous rocks (e.g., mudstones, siltstones, claystones, clayey shales, marls-marlstones), highly/completely weathered/altered rocks (e.g., weathered granite, peridotite, basalt), some qualities of peridotite-serpentinite rockmasses in ophiolitic complexes [fractured peridotites or schistose serpentinites, poor or very poor quality of sheared serpentinite ], heavily jointed and/or tectonically sheared rocks, often forming a chaotic structure (e.g., silty/clayey flysch, ophiolites, heavily broken-disintegrated limestones-dolomites), some types of flysch formations [from type V to type XI, (Marinos 2014)], some qualities of molassic formations [heavily broken or brecciated in fault zones (Marinos 2017)]. The failure mechanism-behavior type in tunneling of such geomaterials is dictated by stress-controlled failures [i.e., shear failures (Sh), squeezing (Sq)] according to the assessment of the behavior of the rockmass in tunneling, as presented by Marinos (2012) and Marinos (2014). In this context, the behavior of the rockmass in tunneling can be regarded as quasi-isotropic; therefore, herein, the rockmass is simulated as a linear elastic-perfectly plastic material following the Mohr-Coulomb (MC) failure criterion. The strength parameters of the surrounding rockmass for the MC failure criterion are calculated according to the fitting process between the Mohr-Coulomb and the Generalized Hoek-Brown failure criteria for tunnels, as proposed by Hoek et al. (2002). Respectively, the elasticity modulus of the surrounding rockmass is calculated, as proposed by Hoek et al. (2002) and Hoek and Diederichs (2006). The parametrically examined ratios of the rockmass strength over the initial (in-situ) stress (r cm /p o ) vary from 0.1 to 1, corresponding to tunneling within weak rockmasses. In such conditions, the tunneling issues encountered are ranked from minor support problems to major ''squeezing'' problems of escalating level, as r cm /p o reduces, according to the relationship between the anticipated radial strain for unsupported tunnels within ''squeezing'' rockmasses and the degree of excavation and support difficulty, as presented by Hoek and Marinos (2000). The shotcrete lining is modeled as a linear elastic material and is assumed to be relatively stiff, with thickness t sh = t = 0.4 m and elasticity modulus E sh = E = 20,000 MPa, despite the variation of r cm /p o . The influence of a more stiff shotcrete shell (t sh = 2t and E sh = E, t = t sh and E sh = 1.5E) on the interaction effect is also evaluated. Benchmark checking, not presented herein, has shown that a less stiff shotcrete shell (i.e., t sh = 0.5t and E sh = E) has a non-substantial (practically negligible) influence on the percentage differentiation of the tunnel loads due to the interaction.
A typical sequential excavation and support method is simulated for both tunnels, assuming either a full-or a partial-face excavation with an advance rate equal to 1 m per step and installation of the primary support at a distance of 1 m behind the tunnel face. The length of the unsupported tunnel span is assumed to be constant, despite the variation of r cm /p o . The full-face construction sequence is applied for either circular or non- Fig. 1 3D numerical models with W/D = 2: a circular and b non-circular tunnels circular tunnels, while the partial-face is implemented exclusively for non-circular. The adopted construction sequence assumes that the excavation and support of the right (denoted as ''second'') tunnel starts after the left (denoted as ''first'') tunnel has been excavated, supported and reached a steady-state. This sequence differentiates only if the effect of the lagging distance (L) is examined. To quantify the influence of the right over the left branch of the tunnels, the results corresponding to the stabilized state of the excavated and supported ''first'' branch (left tunnel), before constructing the ''second'' branch (right tunnel), are denoted as ''single'' and used as the reference. Figure 2 presents the clockwise orientation of angle H at the sections of the left (''single''/''first'') and the right (''second'') tunnel used in representing the spatial distribution of the axial forces. Due to the correlation of the (normal stress) load (p) with the hoop axial force (N) acting on the shotcrete lining via the equation N = p Á D/2, the presented results focus exclusively on the axial force. In this context, the term load represents the axial force. Table 1 summarizes the values of the examined parameters for the rockmass and the primary support and presents the notations used in this paper to facilitate the interpretation of the results. In total, approximately 1000 parametric 3D numerical analyses have been conducted. Figure 3 presents the axial force acting on the primary support of the ''single'' and the ''first'' tunnel for indicative numerical analyses with r cm /p o = 0.2, thus, corresponding to ''squeezing'' geotechnical conditions Barla 2001;. Figure 3a shows the spatial distribution of the axial force at the typical section of the tunnel; at a monitoring section sufficiently behind the tunnel face where the profile of the axial force has been stabilized. The results are presented for three pillar width ratios (W/D), namely 0.5, 1 and 2, indicating an increase in the axial force acting on the ''first'' tunnel. This increase develops asymmetrically and maximizes at the springline (H = 90°), which neighbors with the area of the pillar. It minimizes at the springline (H = 270°), which is on the opposite side of the area of the pillar. Thus, the interaction effect is more intense in the half part of the section of the tunnel (H = 0°-180°), which is adjacent to the area of the pillar. In this context, Fig. 3b shows the longitudinal distribution of the axial force along the springline (H = 90°) of the ''first'' tunnel and compares it with that of the ''single'' tunnel. Both graphs identify a strong dependency on the degree of the interaction effect with the pillar width ratio (W/D). In particular, the relationship between the interaction effect   (additional net loading) versus the pillar width ratio (W/D) is inversely proportional. Figure 4 presents the evolution of the axial force acting on the primary support of the left and the right tunnel at monitoring sections in the middle of their excavation lengths (S/D = 5). Thus, it examines the effect of the pillar width ratio (W/D) on the axial force developing at the springlines (H = 90°or 270°) of the twin tunnels during their construction sequence. The construction of the left tunnel starts at step 1 and finishes at step 80, while the correspondent initial and final steps for the right are 81 and 160, respectively. Therefore, the time-history output presented for steps 1 to 80 corresponds to the ''single'' tunnel, while that presented for steps 81 to 160 refers to the ''first'' and the ''second''.

Effect of the Pillar Width Ratio (W/D)
The results indicate that the load profile of the ''second'' tunnel is practically identical to that of the ''single'', irrespective of the pillar width ratio (W/D). Towards that end, the effect of twin tunneling on the load state of the ''second'' tunnel is not further systematically evaluated in this paper, which mainly focuses on the ''first''. On the contrary, a considerable effect is observed on the load distribution of the ''first'' tunnel. The latter is attributed to the fact that the ''first'' tunnel has already been excavated and supported during the construction of the ''second'' and therefore, the stress and strain redistribution triggered by the advance of the ''second'' tunnel, in the area of the pillar, is directly ''transformed'' into developing additional load on the primary support of the ''first''. More accurately, the stiffness of the shotcrete lining of the ''first'' tunnel prevents the development of additional convergence in the section of the tunnel due to the stress and strain alterations developing in the area of the pillar during the construction of the ''second'' tunnel and condenses the interaction effect into a load transfer. In particular, the load profile of the ''first'' tunnel indicates that the interaction effect begins at step 100, where the advancing face of the ''second'' tunnel is approximately 1.25 diameters (1.25D) behind the monitoring sections. As the ''second'' tunnel advances, the load of the ''first'' tunnel increases, while the interaction effect stabilizes at step 140, where the advancing face of the ''second'' tunnel is about 1.25 diameters (1.25D) in front of the monitoring sections.  at the typical section for the ''single'' and the ''first'' tunnel for indicative numerical analyses with r cm /p o = 0.2, thus, corresponding to ''squeezing'' geotechnical conditions. For K o = 0.5, the stress regime prevailing in the ''first'' tunnel (compared to that of the ''single'') due to the interaction effect is that of purely additional loading at the entire section of the tunnel. This response is irrespective of the pillar width ratio (W/D). For K o = 1 and 1.5, the stress regime developing in the ''first'' tunnel splits into two regions. The first-purely loading-region (I) is that of additional loading (marked with a white fill), despite the pillar width ratio (W/D). The second-mixed loading-region (II) is that of either additional loading or unloading (marked with a grey fill) depending on the pillar width ratio (W/D). The first region (I) develops at the springlines of the tunnel and their neighboring areas, while the second region (II) at the crown/invert of the tunnel and their adjoining areas, respectively. As K o increases, the extent of the second region (II) and the unloading trend within the region are intensified, while also the unloading trend prevails, despite the pillar width ratio (W/D). Therefore, Fig. 5 indicates that the most adverse loading conditions for the ''first'' tunnel develop for K o = 0.5, while, on the contrary, the most favorable for K o = 1.5, in terms of additional net loading. It also demonstrates that K o influences both the extent of the change and the shape of the distribution of the axial force for the ''first'' tunnel. Before the construction of the ''single'' tunnel, for K o = 0.5, the initial (in-situ) horizontal stress (r h ) is lower than the vertical (r v ); thus, the horizontal is the minimum principal stress (r 3 ) and the vertical is the maximum (r 1 ). Therefore, at the springlines of the tunnel, the horizontal stress is the radial (r h = r R ) and the vertical stress is the tangential (r v = r h ), while at the crown/invert of the tunnel, the stress condition is the opposite. After the construction of the ''single'' tunnel, at the springlines of the tunnel, the horizontal/radial stress reduces and the vertical/tangential increases. The opposite redistribution pattern develops at the crown/invert of the tunnel. Therefore, the deviatoric stress, the accumulation of plastic strain (thus, the extent of the plastic zone) and the load acting on the tunnel at the springlines are higher than those encountered at the crown/invert. The latter response is also intensified by the deformational pattern of the section of the tunnel due to its excavation and support, being inwards at the crown/invert and outwards at the springlines. Thus, the higher convergence (indicating higher deconfinement) occurs at the crown/invert and the lower (indicating lower deconfinement) at the springlines. For K o = 1, the initial (in-situ) stress field is uniform leading to a practically uniform redistributed stress and strain field and load for the ''single'' tunnel. For K o = 1.5, the response is the opposite of that for K o = 0.5.
As a result, the area of the pillar is more significantly affected by the construction of the ''single'' tunnel, as K o decreases, exhibiting higher load and plastic strain. It is also subjected to secondary stress and strain redistribution due to the construction of the ''second'' tunnel, leading to a comparatively higher net load transfer on the ''first'', as K o decreases. For instance, for K o = 0.5, during the construction of the ''second'' tunnel, the evolution in the plastification of the already severely affected area of the pillar magnifies the interaction effect on the primary support of the ''first''. On the contrary, for K o = 1.5, the construction of the ''second'' tunnel causes secondary stress and stain redistribution at the comparatively not so intensely affected area of the pillar. Thus, the interaction effect is mitigated and can partially act favorably by decreasing the reference load of the ''single'' tunnel at some regions of the section of the ''first''.   Figure 6 presents the quantified net interaction effect on the axial force developing in the primary support of the left tunnel due to the construction of the right. In particular, it shows the distributions of the ratios of the axial force at representative locations in the typical section [springlines (H = 90°and 270°) and crown (H = 360°)] of the ''first'' over the ''single'' tunnel, relative to r cm /p o , for the numerical analyses category and parameters shown in the graphs heading and legend. The ratios exhibit differentiations at the representative locations in the typical section since their relative arrangement versus the area of the pillar differs. At the springline neighboring with the area of the pillar (H = 90°), the interaction effect is maximized. It causes a purely additional loading increase at this location, despite K o or W/D. However, the level of the additional net loading increase is determined by the combination of r cm /p o , K o and W/D; it is amplified, as r cm /p o , K o and/or W/D decrease. At the springline, located at the opposite side and not adjacent to the area of the pillar (H = 270°), the interaction effect exhibits an analogous trend and is similarly affected by the combination of r cm /p o , K o and W/D. However, the additional net loading increase is comparatively lower than that of the other springline. At the crown (H = 360°), the interaction effect exhibits a mixed response depending primarily on K o and secondarily on W/D and r cm /p o . For K o = 0.5, a purely additional loading increase is encountered, despite W/D. However, the level of the additional net loading increase is also determined by the combination of r cm /p o and W/D. For K o = 1 and 1.5, either loading increase or unloading (loading decrease) is encountered, depending on r cm /p o and W/D. For K o = 1, the interaction effect mainly causes unloading, except for some combinations of r cm /p o and W/D, which cause additional loading. For K o = 1.5, the interaction effect mainly causes unloading, despite the combinations of r cm /p o and W/D. Furthermore, for both K o = 1 and K o = 1.5, the level of unloading is mitigated, as W/D increases.

Design Charts (Effects of r cm /p o ,K o ,W/D)
To sum up, the ratios of the ''first'' over the ''single'' tunnel exhibit an increasing trend versus the decrease of the geotechnical conditions ratio (r cm /p o ), the horizontal stress ratio (K o ) and the pillar width ratio (W/D). In particular, for r cm /p o C0.5, the maximum interaction effect exhibits a constant distribution versus the improvement of the geotechnical   1.5 0.5 1.7-0.9 (1) 0.9-1 (2) 1.3-0.9 (1) 1-0.9 (2) 1-0.8 1-0.8 1 1.4-0.9 (1) 0.9-0.8 (2) 1.2-0.9 (1) 1-0.9 (2) 1-0.8 1-0.8 2 1-0.9 (2) 0.9 (2) 1-0.9 (2) 1-0.9 (2) 1-0. conditions; the constant level of the ratios depends on the variation of K o and/or W/D. Respectively, for r cm / p o B 0.5, it exhibits an increasing distribution versus the degradation of the geotechnical conditions; the variable level of the ratios is enhanced by reducing K o and/or W/D. Furthermore, the maximum interaction effect is strongly influenced by the pillar width ratio (W/D) and more accurately, it is magnified, as W/D decreases. The effect of W/D exhibits a significant relation to K o . In more detail, assuming as the criterion for the reduction of the interaction effect at practically negligible levels the ratios to be N ''first'' /N ''single'' B 1.2, for K o = 1 and 1.5, for W/D = 2, the ratios are marginally 1.2 or lower, indicating the elimination of the interaction. However, for K o = 0.5 and r cm /p o B 0.2, a further increase of the pillar width ratio (W/D) to 3 or 4 is required for the obliteration of the interaction, while for r cm /p o [0.2, a pillar width ratio (W/D) of 2 is adequate.
Complementary to the ratios presented at the representative locations, Table 2 summarizes the range of the maximum and average ratios of the axial force at the typical section of the ''first'' and the ''second'' over the ''single'' tunnel. Thus, an additional quantification and evaluation of the net interaction effect are provided, leading to the observation that the interaction does not crucially influence the ''second'' tunnel (0.8 B N ''second'' /N ''single'' B 1.2). Figure 7 presents the distribution of the axial force N ''first'' (H) at the periphery of the typical section of the ''first'' tunnel (H = 0°-360°), relative to angle H, normalized versus the maximum axial force N ''first' ',max (H = 90) , for all the examined values of r cm /p o , for the numerical analyses category and parameters shown in the graphs heading and legend. The maximum, average and minimum envelopes of the distributions are plotted in the graphs, providing the range of the potential loading conditions. The limits of the distributions depend on r cm /p o ; the maximum envelopes correspond to the lower values of r cm /p o , while the minimum to the higher values of r cm /p o , respectively. The scatter of the distributions is approached, via nonlinear regression, by analytical equations with the following lorentzian general expression:

Analytical Equations (Effects of r cm /p o , K o , W/D)
The N ''first'',max (H=90) of the ''first'' tunnel can be calculated by the output presented in Fig. 6 with required input the N ''single'' (H=90) of the ''single'' tunnel and the geotechnical conditions ratio (r cm /p o ). The constants (a 1 , b 1 , x 1 , y 1 ) of the analytical equations are summarized in Table 3. The analytical equations differentiate for H = 0°-180°and H = 180°-360°and are provided for K o = 0.5 and 1. For K o = 1.5, the scatter of the distribution is more complicated and cannot be strictly approached using analytical equations (although it can be used for calculations).
The ratios of the quantified net interaction effect on the axial force of the ''first'' tunnel versus r cm /p o , W/D and K o presented in Fig. 6, together with the distributions and the equations presented in Fig. 7; Table 3, respectively, constitute a tool for preliminary calculations of the additional net loading (mainly) and unloading (rarely) on the ''first'' tunnel. Despite the method used to define the axial force acting on the ''single'' tunnel (i.e., empirical or closed-form analytical solutions, two-or three-dimensional numerical analyses, etc.), the ratios in conjunction with the equations can provide the required coefficients to modify the distribution of the axial force of the ''single'' tunnel and estimate the correspondent of the ''first''. Figure 8a presents the effect of the diameter of the tunnels on the distributions of the maximum and average ratios of the axial force at the typical section of the ''first'' over the ''single'' tunnel, relative to r cm / p o , for the numerical analyses category and parameters shown in the graphs heading and legend. The diameter of the tunnel is an essential tunneling parameter with a complicated role to be quantified/''decrypted'' since it influences the development of loads in multiple and contradictory ways. Indicatively, it influences the overburden height ratio (H/D), which resembles the ability of the surrounding rockmass to develop the arching response around the section of the tunnel and affects the deconfinement before the application of the primary support in conjunction with the distance of the latter from the tunnel face (x/D). Moreover, it alters the rigidity/flexibility of the primary support, resulting in either increasing or decreasing the loads. It also has a crucial effect in increasing the loads (especially, as the geotechnical conditions deteriorate), since it determines the volume of the mobilized part of the surrounding rockmass, where the plastic zone develops and loads the primary support. The latter issues trace the multiparametric effect of the diameter on the loads of the ''single'' tunnel. However, the distributions plotted in the graphs of Fig. 8a indicate that the quantified net interaction effect is not substantially influenced by the variation of the diameter, despite r cm /p o and W/D. Comparing the ratios for different diameters is not examined for a more expanded range of overburden height ratios (H/D) than those presented. Hence, the proportionality incorporated in the relative variation of the overburden height ratio (H/D) due to changing the diameter (D) for constant overburden height (H) is constant, reasonably presuming that it does not critically affect the level of the interaction.

Effect of the Deformability of the Surrounding
Rockmass (E m ) and the Stiffness of the Primary Support (t sh , E sh ) Figure 8b presents the effect of the deformability of the surrounding rockmass (E m ) and the stiffness of the primary support (t sh , E sh ) of the tunnels on the distributions of the maximum and average ratios of the axial force at the typical section of the ''first'' over the ''single'' tunnel, relative to r cm /p o , for the numerical analyses category and parameters shown in the graphs heading and legend. The distributions plotted in the graphs indicate that increasing the stiffness of the primary support by increasing either the thickness or the elasticity modulus causes a slight decrease in the quantified net interaction effect. The latter response is comparatively enhanced, as W/D and/or r cm /p o decrease. In more detail, increasing the stiffness of primary support decreases the deconfinement of the surrounding rockmass (at the stage where the interaction of the system surrounding rockmassprimary support reaches a steady-state), leading to increased axial forces for the ''single'' tunnel. The latter effect on the deconfinement of the surrounding rockmass slightly mitigates the intensity of the interaction between the twin tunnels; thus, it slightly controls the additional stress and strain redistribution developing in the area of the pillar, resulting in also decreased axial forces for the ''first'' tunnel. Therefore, the net interaction effect exhibits a decreasing trend. Respectively, the distributions plotted in the graphs indicate that increasing the stiffness of the surrounding rockmass by increasing the elasticity modulus causes an increase in the quantified net interaction effect. The latter response is comparatively enhanced, as W/D decreases, while it is constant, despite r cm /p o . More accurately, increasing the elasticity modulus of the surrounding rockmass, via its calculation according to Hoek et al. (2002) than Hoek and Diederichs (2006), reduces the maximum potential convergence of the surrounding rockmass (corresponding to full deconfinement), in terms of absolute magnitude. In other words, it alters the slope of the reaction curve of the surrounding rockmass. Therefore, it imposes lower loads on the primary support of the ''single'' tunnel, since the latter undertakes or counterbalances lower maximum potential convergence of the surrounding rockmass. The latter effect mitigates the intensity of the interaction between the twin tunnels. As a result, it blunts the additional stress and strain redistribution developing in the area of the pillar, thus inducing also lower loads on the primary support of the ''first'' tunnel. Therefore, the net interaction effect exhibits a constant or mainly increasing trend. However, the loads of both the ''single'' and the ''first'' tunnel decrease, in terms of absolute magnitude, as the stiffness of the surrounding rockmass increases.

Effect of the Lagging Distance (L)
The effect of the lagging distance ratio (L/D) is also examined, via three scenarios with different construction sequence of the tunnels: (i) for L/D = 10, the construction sequence of the right tunnel starts after the excavation and support of the left tunnel has been completed (reference scenario, already discussed), (ii) for L/D = 5, the construction sequence of the right tunnel starts when half of the excavation and support of the left tunnel has been completed, (iii) for L/D = 0, the construction sequences of both (left and right) tunnels are assumed to be simultaneous (''fictitious'' scenario not applicable in tunneling practice, exclusively aiming to examine its potential effect).
In this context, Fig. 9 presents the evolution of the axial force acting on the primary support of the left and the right tunnel at monitoring sections in the middle of their excavation lengths (S/D = 5), for L/D = 0, 5, 10 and W/D = 0.5, 1. For L/D = 10, the construction of the left tunnel starts at step 1 and finishes at step 80, while the correspondent initial and final steps for the right are 81 and 160, respectively. For L/D = 5, the construction of the left tunnel starts at step 1 and finishes at step 80, while the correspondent initial and final steps for the right are 41 and 120, respectively. For L/D = 0, the correspondent initial and final steps for both (left and right) tunnels are 1 and 80, respectively. Figure 10 presents the effect of the lagging distance ratio (L/D) on the distributions of the maximum and average ratios of the axial force at the typical section of the ''first'' over the ''single'' tunnel, relative to r cm /p o , for the numerical analyses category and parameters shown in the graphs heading and legend. The plots in the graphs indicate that the quantified net interaction effect practically coincides or negligibly differs for L/D = 10 and 5, despite r cm / p o and W/D. On the contrary, the simultaneous construction of both (left and right) tunnels (L/ D = 0) proves favorable for the left and adverse for the right, respectively, since the quantified net interaction effect coincides for both tunnels, despite r cm /p o and W/D. However, the level of the quantified net interaction effect, for L/D = 0, is decreased compared to that for L/D = 10 and 5, depending on r cm /p o and W/D; the level of reduction exhibits an increasing trend, as r cm /p o and/or W/D decrease.

Circular Versus Non-circular Tunnels-
Effects of the Section (Shape) and the Construction Sequence (Fullor Partial-Face) of the Tunnels 4.4.1 Effect of the Lagging Distance (L) Figure 11 presents the effect of the section (shape) and the construction sequence (full-or partial-face) of the tunnels on the evolution of the axial force acting on the primary support of the left and the right tunnel at monitoring sections in the middle of their excavation lengths (S/D = 5), for W/D = 0.5. For circular and non-circular tunnels and full-face construction sequence, (C*) and (N-C*), respectively, the construction of the left tunnel starts at step 1 and finishes at step 80, while the correspondent initial and final steps for the right are 81 and 160, respectively. For non-circular tunnels and partial-face construction sequence (N-C** and ***), the construction of the top-heading phase (N-C**) of the left tunnel starts at step 1 and finishes at step 80, while the bench phase (N-C***) starts at 81 and finishes at 160. The correspondent initial and final steps for the topheading phase (N-C**) and the bench phase (N-C***) of the right tunnel are 161 and 240, 241 and 320, respectively. The distributions plotted in the graphs indicate that the maximum net load transfer is either comparable or slightly reduced for non-circular than circular tunnels if a full-face construction sequence is applied. Furthermore, the maximum net load transfer is slightly mitigated for partial-than full-face construction sequence for non-circular tunnels. The latter points are valid for the indicative numerical analyses presented Fig. 11 corresponding to ''squeezing'' geotechnical conditions (r cm /p o = 0.2) and a narrow pillar width ratio W/D = 0.5. 4.4.2 Effects of r cm /p o and W/D Figure 12 presents the effect of the section (shape) and the construction sequence (full-or partial-face) of the tunnels on the distributions of the ratios at representative locations in the typical section [springlines (H = 90°and 270°) and invert (H = 180°)] of the ''first'' over the ''single'' tunnel, relative to r cm /p o , for the numerical analyses category and parameters shown in the graphs heading and legend. The distributions plotted in the graphs indicate that at the springline with H = 90°, the interaction effect is either comparable or slightly reduced for non-circular than circular tunnels if a full-face construction sequence is applied, with the differentiation level intensifying, as W/D and r cm /p o decrease. Respectively, for noncircular tunnels, at the springline with H = 90°, the partial-compared to the full-face construction sequence causes mainly either comparable or increasing interaction effect, with the deviation degree exhibiting a relation to the combinations of W/D and r cm /p o and the construction phase (top heading and bench). However, the distributions of the maximum quantified net interaction effect, developing for H = 90°, despite the assumed section and the adopted construction sequence of the tunnels, orient distinct, well-established, scatter clouds versus the pillar width ratio (W/D). At the springline with H = 270°and the invert (H = 180°), the interaction effect is systematically more intense for non-circular tunnels than circular, despite the construction sequence, with the latter response relatively mitigating, as W/D increases. Figure 13 presents the optimum pillar width ratio (W/ D) optimum that drastically reduces the interaction effect, relative to r cm /p o and K o . The assumed criteria for the elimination of the interaction effect and the identification of the optimum pillar width ratio (W/ D) optimum are the ratios of the axial force to be N ''first'',H=90 /N ''single'',H=90 \ 1.1 and N ''first'',max / N ''single'',max \ 1.1. The maximum, average and minimum level of the scatter of the distributions is approached, via non-linear regression, by analytical equations with the following exponential general expression:

Circular and Non-Circular Tunnels-Optimum Pillar Width Ratio (W/D) optimum (Nomographs and Analytical Equations)
The constants (a 2 , b 2 , y 2 ) of the analytical equations are summarized in Table 4. For K o =0.5 and 1, the two adopted criteria produce identical or slightly differentiated results for (W/D) optimum since the maximum load for both the ''single'' and the ''first'' tunnel develop at the side with H = 90°. On the contrary, for K o =1.5, the two criteria lead to significantly differentiated results for (W/D) optimum since the maximum load for the ''single'' tunnel develops at the crown/ invert (H = 360°/180°), while that for the ''first'' at the side with H = 90°.

Conclusions
This paper utilizes parametric 3D Finite Element analyses to calculate the interaction effects between deep, circular and non-circular, parallel-twin tunnels, excavated and supported with a sequential method, within weak rockmasses. The investigation focuses on the axial forces developing in the primary support of the tunnels. The results of the analyses indicate that the construction of the new (''second'') tunnel triggers interaction effects, which lead to an increase of the loads acting on the shotcrete lining of the existing (''first'') tunnel. The load distribution of the ''first'' tunnel is asymmetric and exhibits its maximum at the springline with H = 90°neighboring the area of the pillar between the tunnels. On the contrary, the load distribution of the ''second'' tunnel is practically identical to that of a ''single'' tunnel, with the latter serving as a reference to quantify the interaction effects.
The quantified interaction effect, reflected via the ratios of the axial force of the ''first'' over the ''single'' tunnel, exhibits an increasing trend versus the reduction of the geotechnical conditions ratio (r cm /p o ), the horizontal stress ratio (K o ) and the pillar width ratio (W/D). In particular, for the r cm /p o C 0.5 regime, the distribution of the ratios is constant with their level depending on the variation of K o and/or W/D. Respectively, for the r cm /p o B 0.5 regime, the distribution of the ratios is variable, increasing versus the deterioration of the geotechnical conditions with an additional escalating trend for their magnitude, as K o and/or W/D reduces. The effect of W/D exhibits a substantial dependency on K o . More accurately, for K o = 0.5 and r cm /p o B 0.2, the orientation at practically negligible levels of the interaction effect (establishing as criterion N ''first'' /N ''single'' B 1.2) requires a pillar width ratio (W/D) of 3 or 4. On the contrary, for K o = 0.5 and r cm /p o [ 0.2, K o = 1 and 1.5, irrespective of r cm /p o , a pillar width ratio (W/D) of 2 is sufficient. Therefore, the most adverse loading conditions for the ''first'' tunnel prevail for K o = 0.5, while, on the other hand, the most favorable occur for K o = 1.5, in terms of additional net loading. Several relevant parameters, as shown in Table 1, are varied to produce dimensionless graphs (design charts) of the ratios of the axial force, relative to r cm / p o , W/D and K o , also incorporating the effect of other tunneling parameters such as the diameter (D) and the overburden height (H). The aforementioned normalized plots constitute a tool for preliminary calculations of the net load change at representative locations of the ''first'' tunnel for a known correspondent load of the ''single''. Furthermore, analytical equations are provided for calculating the distribution of the axial force along the periphery of the typical section of the ''first'' tunnel, normalized versus its maximum axial force. Therefore, the ratios combined with the equations can provide the required coefficients to acquire the distribution of the axial force for the ''first'' tunnel based on the correspondent for the ''single''. In addition, nomographs and analytical equations are provided for estimating the optimum pillar width ratio (W/D) optimum , which reduces the interaction effect in practically negligible levels, relative to r cm /p o and K o .
Furthermore, two longitudinal separation distances between the tunnel faces are simulated, equal to L/D = 10 and 5, raising identical effects on the magnitude of the interaction. On the contrary, the simultaneous construction of both tunnels, representing the case of zero lagging distance (L/D = 0), although it is an idealized scenario, results in being favorable for the ''first'' and adverse for the ''second'' tunnel, especially, as r cm /p o and/or W/D decrease, by comparing their loads with those of the ''single''.
Moreover, for full-face construction sequence of twin tunnels, the maximum interaction effect (thus, at the springline with H = 90°) for a non-circular tunnel shape can be adequately provisioned using a circular section. On the contrary, for partial-face construction sequence (and non-circular tunnel shape), the correspondent maximum interaction effect is either comparable or slightly increased, with the differentiation being associated with the combinations of W/D and r cm /p o and the construction stage (top heading and bench).
Finally, the applicability area of the proposed methodology is oriented for (disturbed, disintegrated, sheared) weak rockmasses exhibiting a quasi-isotropic behavior in tunneling; thus, their response is dominated by stress-controlled potential failures. On the other hand, for jointed rockmasses, whose response in tunneling is dictated by anisotropy and gravity-driven potential failures controlled by discontinuities, numerical analysis simulating joints is suggested to be carried out. Also, it is pointed out that the proposed methodology does not consider the potential timedependent deformations (creep effect) ensuing in some types of rockmasses, reasonably presumed to deteriorate the interaction between twin tunnels. K o = 1.5 N ''first'',max /N ''single'',max \ 1.1 Avg 2.10 9.80 0.50