Stability Analysis and Support Design of Imam Reza Tunnel in Ardabil Sarcham Road by Numerical Methods

Geological structures and performance of the geodynamic processes can affect engineering projects on their own. Hence, the stability analysis and designing methods for foreseeing the retaining and support system for tunnels are diverse and came from different points of view. So this study seeks to present stability analysis of Imam Reza tunnel in Ardabil Sarcham Road with a special focus on the impact of future earthquakes on its stability using numerical methods. In this study, first designing and operating the initial structure with the height of 5.5 m and a semi-circular cross section. Secondly, drilling with the height of 3m and the width of 7.34 m and with a rectangle cross section. For stabilization, Rock Mass Rating (RMR) geomechanical classification systems and methods used. At the stabilization level, the materials were examined in laboratory, regarding the properties of sides and roof of the tunnel and pressure on them. The results of physical and mechanical experiments shown that the compressive strength ranged from 400 kg/cm to 500 kg/cm on average. The elastic modulus is between 12 and 13 GPa for the rocks. The Cohesion (C) ranged from 4-5MPa to 5 MPa and the Angle of Internal Friction (φ) is between 60o and 50o.


Introduction
The tunnel industry considered that tunnels were naturally resistant to earthquake action for many years, as they did not experience the same high levels of shaking as surface structures.This perception was supported by the relative good historic performance of tunnels and underground structures, especially of tunnels in rock, during large earthquakes [1].
One of the most significant current discussions in tunneling is stability.In recent decade, numerical modeling has been introduced to check the stability.When the tunnel has been made, stability is an important factor due to advancing cycle and cost especially in soil or weak rock [2].
One of the most important factors for designing a tunnel or other underground spaces is to provide the stability.Due to the fact that boring a tunnel requires relatively high costs, any mistake in studying, analyzing, and operating the support systems as well as utilizing them improperly, leads to wasted capital and furthermore, it prevents us from reaching our goals.Regarding these explanations, it's necessary to suggest a single or mixed support system through a complete analysis of different geological and geomechanical parameters so that it will get the highest efficiency with the lowest cost.Full check of the properties of rock mass and its situation is one of the crucial issues of reaching this purpose.
Methods developed to the stability analysis of slopes during earthquakes can be divided into three categories: (1) pseudo-static methods, (2) dynamic sliding block methods, and (3) stress-strain methods [3].Terzaghi (1950) first presented the pseudo-static method, which is a simple method for evaluating of seismic stability of a slope [4].
A preliminary study of the engineering stability of the Shilou tunnel in the Hipparion red clay were studied by the establishing of finite element numerical model [5].
Up to now, many research workers who major in the tunnel project bring a large number of mathematical analysis methods into the analysis of the stability such as block limit equilibrium method [6][7], imbalance thrust method [8], Principle of red flat projection analysis [9].This study has two stages.First stage includes designing and operating the initial structure in the upper part with the height of 5.5m and a semi-circular cross section.The second stage includes drilling in the lower part with the height of 3m and the width of 7.34m and with a rectangle cross section.Also at the stabilization level, the used components and materials were examined in laboratory in each stage, regarding the properties of sides and roof of the tunnel as well as the pressure on them.
So this study seeks to present stability analysis of Imam Reza tunnel in Ardabil Sarcham Road with a special focus on the impact of future earthquakes on its stability using numerical methods.

Literature Review
Generally, the stability of tunnels and the designing process of appropriate support systems for them have been of great importance in mining and construction projects since a long time ago.
The topic of modelling tunnels and slopes has drawn experts' attention to itself since 1970.Recently, accurate methods for calculating the probability of collapse of tunnel walls or its roof have been made which are developing in terms of applicability.Studies in this field which have been done in Iran are: Fasihnikoutalab et al. ( 2012) conducted a study entitled numerical stability analysis of tunnel by PLAXIS.The advantages is that we can have an appropriate calculation of behavior of soil or rock before drilling.Also in order to consider part of tunnel is under a main settlement risk, especially the tunnel is drilled into multi-layer soils.PLAXIS considered as most well-known software in the field of Finite Element method that can be analysis the behavior of tunnel and the amount of ground movement.According to the calculated thickness of lining, 35 cm thickness observed as optimum thickness of coverage system, based on the analysis the most displacements is obtained in ceiling and floor that is equal to 9.65E-03 meters.The main stress estimated in the model in addition to the lining is -478 kN/m 2. Stress is returned to normal condition in 13 meters distance from the center of the tunnel [10].
Ming Xiao and Chen Zhao (2017) in a paper, a coupling method for simulating the hydromechanical interaction of the reinforced concrete lining is adopted to investigate the effect of internal water exosmosis on the seepage field at the collapse zone and to determine the external water pressure on the steel lining.The findings reveals that a 10 m thick consolidation layer and the adjusted antiseepage measures confirm the stability of the steel lining at the collapse area under internal and external water pressure.These research findings show a reference for the design of treatment measures for large-scale collapses in steel-lined pressure tunnels [11].Cheng et al (2017) were established models of structure and fluid fields by considering the effect of a twodimensional viscoelastic boundary condition, and the stability of a cross-sea tunnel structure under seepage and a bidirectional earthquake was studied.The results showed that under this type of earthquake, plastic zones first appeared on both sides of the tunnel arch feet and in the peripheral areas of the arch shoulders; by contrast, the vault was secure, and no plastic zone appeared in this area.When all factors were considered, the plastic zones extended farther when the overlying rock thickness and sea depth were greater.When the overlying rock was thick and the rock mass surrounding the tunnel was only slightly weathered, the rock mass was close to intact, and as a consequence, the permeability coefficient did not significantly affect the development of the plastic zone or the stability-safety factor [12].Yang et al (2017) in a study investigated the stabilities and associated collapse mechanisms of deep buried unlined elliptical tunnels in cohesive-frictional soils with the action of soil weight by the "upper-bound finite element method with rigid translatory moving elements" (UBFEM-RTME).The soil masses are assumed to obey the Mohr-Coulomb yield criterion and a nonassociated flow rule.Upper-bound stability coefficients are deduced for different values of friction angles, dilatancy coefficients, and dimensionless spans.The obtained collapse mechanisms do not extend to the ground surface and are primarily composed of a series of mutually movable rigid blocks [13].
Wang et al (2017) analysis soil excavated hole portion of the hole and the hole body during deformation movement of surface soil slope.The main conclusions are: (1) on the stepped portion of the hole excavation on the scope of the maximum slope excavation slope deformation rate in step fastest, over time, the deformation gradually converge.(2) along the longitudinal direction of the tunnel, the more monitoring points away from the hole excavation unit distance, the smaller the amount of deformation of the slope obtained its monitoring.(3) by Yang slope settling cloud contrast slope under different cavity length of the inner body can be seen: With the entry of the mountain tunnel and increase the body, affecting Yang slope deformation area also increased to reflect the slope according to the monitoring point.Monitoring data, with the increase of the hole inside the mountain itself, the settlement value of the monitoring points are increased, and finally stabilized [14].Maghous et al (2017) conceived as a potential alternative to the classical design methods employed for analyzing the stability of underground structures driven in jointed rocks, the homogenization approach stems from the heuristic idea that, from a macroscopic point of view, a rock mass cut by a network of joints may be perceived as a homogenized continuum.This criterion is then applied to the investigation of stability analysis of a tunnel excavated in a jointed rock mass.Upper bounds estimated of the stability factor are derived from the implementation of the kinematic approach directly on the homogenized underground structure [15].Li and Zhang (2009) the computation result indicates that the explicit dynamic analysis program could simulate the whole process of the instabilization and failure of the weak rock layer, and avoid the convergence existing in the general finite element computation, and the analysis model and program could offer effective measure for the analysis of local stability of the tunnel surrounding rock mass with weak rock layer [16].
Yang and Zhang (2017) conducts the plane-stain stability analysis of unlined dual shallow-buried tunnel roof with considering settlement.Due to the geomechanical difference between different parts of twin tunnels, analytical expressions of the two different shape curves which describe two different parts of neighborhood tunnels are obtained with the variation calculus, which are verified by numerical modelling.The dimensionless solutions for the stability number of twin tunnels are obtained by virtual work equations and also compared with those calculated by other collapse mechanisms.The stability number is derived to examine the effects of spacing between tunnels on the stability of tunnels.The parameter sensitivity analysis is carried out to explore the scope of failure soil block and the upper bound limit supporting load in twin tunnels [17].Liu et al (2017) to further reveal the difference between the stress condition and the stability of the tunnel lining of the drainage system under two block patterns, Zhongliangshan Tunnel was selected for field tests.On the basis of the tests, the stress on various parts of the tunnel lining and its variation trend were studied.Then, the numerical model of the tunnel was established with the finite element method and the numerical simulation results were compared with the field test results.On this basis, the stability of the tunnel when the drainage system was partially blocked was analyzed.The findings show a significant difference between the two block plans of the drainage system in how they affect the lining stress.The tunnel stability analysis reveals that when the blocked lengths are 16 and 24 m, the decrease factors are approximately 10.3 and 9.7, respectively, and the tunnel stability is monitored by the symmetric block plans of the drainage system.The tunnel stability is controlled by the asymmetric block pattern.The length of the blocked part determines which block pattern controls the tunnel stability.In the case of asymmetric block, vault limit displacement is not recommended as an instability criterion [18].Goshtasbi and Dehghan (2001) reports the results of a research study which was carried out in order to analyze the stability of a tunnel at Glandroud coal mine, which is situated in northern Iran.Therefore the findings of the analysis reveal that instability in the tunnel resulted of falls of the roof and sidewalls.So, different support systems were evaluated.This result reveals that grouted rock bolt and shotcrete is the main proper support system in Glandroud coal mine tunnel [19].

Research Methodology
In this section we compared endurance time of structures after subjecting them to ET records based on the time of first plastic hinges formation and the time of entering to nonlinear region as well as their inter-story and residual drifts.In addition, maximum inter-story and residual drift results of structures under real earthquake records were evaluated.Based on total results, the best structure with lower weight or cost was determined.
In order to carry out a major construction project before anything else, it is necessary to know the area where the structure is to be constructed.Geotechnical studies are carried out for this purpose and with regard to some important features of rock and soil.Certainly, knowledge of these features requires visiting the site and familiarity with the issues of rock and soil and doing field tests and laboratory tests on the sample.The collection of information collected for this purpose, after the general geological survey of the area, is: A. Drawing geological profiles B. Determine the type of rock and its status C. Measure the physical and mechanical properties of the stone D. Determine the quality of the core drilling or the frequency of rock fractures E. Describe the status of the joints F. Groundwater conditions In this section, according to the above, to study and describe the general geological conditions of the region and its geotechnical studies to understand the behavior of the land and the geometric characteristics of the design are discussed.

Technical Specifications of the BH1 Borehole
According to the investigations carried out at the site of the BH1 borehole and the examination of manual and microscopic specimens at the site of this borehole, according to the specimen, 27 cores with a diameter of 112 mm up to a depth of 38 m are provided.For each core, the geometric characteristics of the layers, core recycling and rock quality index.Studies have shown that the site of this borehole is of igneous rock and it is from acid tuff.In the upper parts of the borehole, up to a depth of 5 meters, the rocks are altered and crushed and it can be seen, but from a depth of 12 meters later, the rocks are massive and uncut.Average percentage of core yield ratios and the stone quality index is 78 and 25, respectively.

Technical Specification of the BH2 Borehole
In this borehole, its drilling is done in a rotary manner, the drilling depth is 22 m and the diameter of the borehole is 112 mm, of which 20 cores are provided.The information obtained from the cores indicates that the site is an acidic igneous rock type (Rhytadasite tuff).In the upper parts of the crustal and severe alteration to a depth of 6 meters, from a depth of 6 to 10 meters of joints in the cores, but from a depth of 10 meters later, the rocks are massive and integrated.The texture of the rocks is from the depth of 6 meters to 19 meters, and is from the depth of 19 to 22 meters of the fragmental tissue.Also, in this borehole, the groundwater level is not encountered; the average percentage of the ratio of the cores' yield is the rock qualitative index of 84 and 46, respectively.

Geometrical Properties of Sections of the Tunnel and the Gallery
In this part, geometrical properties of the cross-section of the tunnel will be studied, as well as technical tips like directing the water around the tunnel, cross slope of the bed, and the necessity of insulation.

Figure 1. Geometric cross-section of the tunnel
As shown in the Figure 1, this typical cross-section consists of an arch with a radius of 5.225 m at the roof and straight inclined wall on the sides.Inclination of the wall has its pros and cons.regarding the construction, the moment in the middle of wall will be reduced and internal forces of the initial and final linings will lean towards the predominant axial force.It will form the most appropriate position for designing.In this case, if the bed rock is weak, accurate underpinning can be used.If the rock is average or good, the rock mass will do the underpinning alone.Consequently, the only weakness of inclination of the side wall, will unfold only when the tunnel is drilled in weak rock and the accurate underpinning is not used for some reasons.Directional drilling of the sideway is not a problem in terms of executive actions and leads to economizing on drilling as well as materials of temporary lining.
In the section of galleries which are designed in the ends of the tunnel, gallery wall construction must lean on alluvial bed that doesn't have enough strength, and creating foundation for it, seems to be necessary.On the other hand, and from executive perspective, inclined formwork leads to many problems.Consequently, it's reasonable to sustain the internal inclined shape of the gallery wall in order to sustain the shape of the mold.But if the external aspect is converted to perpendicular line, the width of the foundation will be reduced and constructing the gallery will become easier.Considering these, Figure 2 shows the section that is used for galleries.In the bottom of the tunnel, a special cross slope has been considered and the cross slopes of both sides depend on general properties of the cross slope of pathway pavement.
The watercourses have been considered in both sides of the tunnel through prefabricated construction method according to presented plans.These sections will be constructed on final basis after completion of the final lining.

Drainage
Longitudinal profile of the tunnel shows that it's necessary to check tunnel drainage in underground excavation parts and end galleries separately in terms of the amount of intruding water and also because of the executive considerations.Tunnel is drilled in a moist area and sometimes with thin streams.Water intrusion into tunnel or back of the temporary lining during the utilization is mostly predictable and embedding an appropriate drainage system for it, is required.Suggested method for rocks and soils has been shown in Figures 2 and 3.According to these figures, drains have been foreseen in excavation section with longitudinal distances of 2m or 3m across the width by which the intruding water in back of the lining will be directed to the internal watercourses.In the ending galleries part that the surrounding areas are alluvial, the amount of water intrusion around gallery is considerable too, and this drainage problem can be easily solved through embedding longitudinal drains on both sides of gallery walls.Also, it's required to fill the back of gallery walls with coarse-grained materials that facilitate streaming the water into the drain.It's necessary to execute, control and finalize the dimensions and the type of drain materials in real situation, due to the fact that the discharge in soil and rock is uncertain at this stage of studies.Considering that the excavation surfaces are moist, an insulator can be used between the initial and final linings in some parts of tunnel longitude in case of observing massive water intrusion, and details about each area are presented on executive plans.Figure 5 shows the insulation section of II area as an example.Pictures of trench walls and accumulated water are shown in Figure 5.This figure shows the insulation within the tunnel too.As shown in the figures, pyroclastic masses are relatively heterogeneous in this part and clean sand lenses are seen without any adhesion at some points, or on the contrary, they are seen with high amount of fine-grained adhesive particles.
Considering this heterogeneity and also the role which clean sand lenses can play in trench sidewalls surface instability, computational cohesion for pyroclastic mass was selected to be 50 KPa in order to determine stable slope of the sidewalls.Indeed, regarding that the impact of local lenses on forces exerted on initial lining is of much less importance because of arch effect and redistribution caused by crusty structure with long longitude, it seems that it's more realistic to select the cohesion to be 100 KPa for designing the initial lining.Obviously, this assumption needs to be controlled through installation of an accurate tool and conducting behavior assessment during performance and also it needs to be corrected if necessary.

Results
In order to assess the seismic properties of the area, time distribution, magnitude distribution, epicenter distribution and focal-depth distribution of earthquakes occurred in past 115 years were analyzed using collected data from tunnel site surroundings within a radius of 100 km.

Time Distribution of the Occurred Earthquakes
In order to study the time distribution of the occurred earthquakes at tunnel site surroundings within a radius of 100km, two dimensional graphs were depicted using Excel Software and data related to year, month, day, hour and number of occurred earthquakes which indicate the time distribution in different years (Figure 6), different months (Figure 7), different days (Figure 8) and different hours and magnitudes (Figure 9 and 10, respectively).

Epicenter and Hypocenter Distribution of the Occurred Earthquakes
The data related to the occurred earthquakes were used in order to assess the epicenter distribution of earthquakes in the tunnel site within a radius of 100 km and also maps of epicenter distribution were created which indicate that no earthquake occurred within a radius of at least 16 km, the closest earthquake to the site in the last 86 years was in a distance of 16 km, and most of the earthquakes occurred in north east of the site with distances of more than 50 km.To evaluate hypocenter of the occurred earthquakes in the surroundings of the tunnel site, two dimensional graphs were depicted using Excel Software and data related to the date and focal depth of occurred earthquakes which indicate their hypocenter distribution within a radius of 100 km.  Figure 11 indicates that in the past 86 years, most shallow earthquakes that occurred around the site of the tunnel, were within a radius of 100 km and at a depth of 10km underground.Also, the deepest earthquake has occurred at a depth of 79 km.Most of these earthquakes occurred at a focal-depth of 14 km.The strongest earthquake occurred at a focal-depth of 45 km and its magnitude was 5 on the Richter magnitude scale.

Computing the Acceleration of Future Earthquakes in Imam Reza Tunnel Site and Assessing the Shaking Potential of the Faults in this Area
The acceleration of future earthquakes in the tunnel site was computed as (  = 1080 ×  0.58× × ( + 25) −1.32 ) through computing accelerations of future earthquakes in their epicenter and using the suggested empirical equation.Table 2 shows the obtained results.
Based on the studies conducted by different researchers, it has been cleared that there is a direct relationship between the shaking potential of the faults and the length of their straight part.On the other hand, the geometrical set up of the faults can influence their shaking potential.Maximum magnitude of earthquake and the expected horizontal and vertical accelerations that are caused by different faults can be computed in every site by determining the shaking potential of the faults.The shaking potential of the faults can be determined using the following equation which is an example of the mentioned equations.Considering the data in Table 2, it's obvious that horizontal acceleration is higher than vertical acceleration and this factor indicates high seismicity in this area.

Continuous Region Analysis (Numerical Method)
Numerical methods are very convenient tools for analyzing the ground behavior caused by boring the tunnel as well as interaction modelling of retaining system and the ground in complicated boundary conditions.The real crosssection of the tunnel has been modelled using Plaxis software in order to control the impact of its shape on the stress field, the deformation occurred around the tunnel and the forces made in the retaining system.This software that has been created for modelling continuous region based on the finite element method, is able to consider the complete soilstructure interaction in the form of a continuous region with the possibility of modelling all types of retaining systems.
In these analyses, ground behavior has been considered as complete elastic-plastic deformation with Mohr-Coulomb yield criterion.Tunnel cross-section excavation is modelled through modelling the elements in the excavation area, and the 3D effects caused by the movement of the excavation face is modelled by reduction of equivalent nodal forces on excavation border to a certain extent that is determined from the corresponding convergence curve and the displacement curve of the wall based on the distance of excavation face ( − ).The remained equivalent nodal forces on excavation border will be reduced to zero through installing the retaining system.Thus, the effect of distancing of excavation face and the exerting force on retaining system will be modelled.5 stages have been considered for the tunnel boring in the upper and lower parts.The cross-section has been analyzed using Plaxis 7.2 and the numerical method (finite element method) with joint modelling of the mass and the retaining system.The borders of the model are extended from both sides and the bottom 3 to 5 times the tunnel diameter and have been considered from top to the earth surface.Support conditions in vertical borders of both sides have been modelled using horizontal rollers that limit the movement of marginal points in horizontal direction, and also have been considered to be closed in lower horizontal border of the horizontal and vertical displacements.Finite element network has been created automatically using the triangular elements with 15 nodes.Figure 12 shows the computer model of the tunnel and the finite element network.As can be seen in Figure 13, this curve has been resulted by removing the elements within the excavation border and exerting unbalanced nodal forces on it.The rate of stress relaxation is equal to 73% regarding the ground interaction curve in the first upper part.This value has been considered to be 70 to make sure of the interaction analysis.Therefore, 70% of the unbalanced nodal forces caused by removing the elements within the excavation border have been balanced before activating the elements of the retaining system in order to consider the 3D effects of the advancement of excavation face on the upper part.The remained unbalanced nodal forces have been balanced after activating the elements of the retaining system and changing the stiffness matrix to consider stiffness of the elements of the lining and thus, the effects of excavation face displacement and exerting forces on the elements of the lining have been modelled.The field of deformation in the surrounding rock mass and the exerting forces on the retaining system are determined once the balance is established.Hence, the stages of excavation modelling and installation of lining in the first part of the upper cross-section excavation have been divided into two phases:  In the second stage of analysis (Phase 2), according to Figures 14 and 15 stress-deformation field of ground mass and the exerting forces on retaining system at the end of creating and completing processes of initial lining in the upper part are obtained by activating the elements of the lining in the model and completely releasing the unbalanced nodal forces within the excavation border (Figures 16, 17 and 18).
According to Figure 16, it shows the total change occurring in the Earth's mass due to the approach of the drilling front to the cross section and the frontal distance from the section to the average distance, in addition the displacement resulting from the installation of the cover.According to Figure 18, it shows that the soil mass around the opening location has elastoplastic behavior.Continuing the completion of the excavation, maximum length without lining in the second stage of upper excavation has been considered equal to 1.5 m regarding the dimensions of the section and stability duration without lining.The initial lining distance from the face is considered equal to 0.5m before starting the excavation.The rate of stress relaxation for simulating the transition and relaxing the stress caused by approaching excavation face to the crosssection and passing it up to a certain distance (equal to 1m on average in here), is determined to be 78% with the purpose of modelling the 3D effects of advancement of excavation face.This value has been considered to be 75% to make sure of the interaction analysis.Therefore, 75% of unbalanced nodal forces caused by removing the elements within the excavation border have been balanced before activating the elements of the retaining system in order to consider the 3D effects of the advancement of excavation face in excavation of the second upper part.After activating the elements of the retaining system and changing the stiffness matrix for considering the stiffness of the elements of the lining, the remained unbalanced nodal forces have been balanced and thus, the effects of excavation face displacement and exerting forces on the elements of the lining have been modelled.The field of deformation in the surrounding ground mass and the exerting forces on the retaining system are determined once the balance is established.Hence, the stages of excavation modelling and installation of lining in the second part of the upper crosssection excavation have been divided into two phases: Phase 3: Removing the elements within the excavation border and releasing 75% of the unbalanced nodal forces.Phase 4: Activating the elements of the lining and releasing the remained unbalanced nodal forces until the establishment of full balance.
Results obtained from analysis of Phase 3, shows the situation of the stresses and deformations in the soil mass after the second stage of excavating the upper part and before installing the retaining system and removing lining of the middle part.In the fourth stage of the analysis (Phase 4), the temporary middle lining is removed and unbalanced nodal forces are completely redistributed in the excavation border by activating the elements of the lining (Figures 21 and 22).Final lining of the upper part is constructed after installing the initial lining of the upper part and before starting the excavation of the lower part by foreseeing the necessary arrangements for concrete placement and its connection to the linings of the lower part.Then, the excavation of middle section of the lower part is done.Modelling the stages of the excavation and installing the lining on the side of the cross-section lower part require determination of the rate of stress relaxation for modelling the 3D effects of the advancement of excavation face.Maximum length without lining in the excavation of the side of the lower part has been considered equal to 4m regarding the dimensions of the section and stability duration without lining.The initial lining distance from the face is assumed to be 2 m before starting the excavation.Therefore, the rate of stress relaxation for simulating the stress transition and relaxing the stress caused by approaching excavation face to the cross-section and passing it up to a certain distance (equal to 3m on average in here), is determined to be 77% with the purpose of modelling the 3D effects of advancement of excavation face similar to the previously applied method.This value has been considered to be 73% to make sure of the interaction analysis.Thus, 73% of unbalanced nodal forces caused by removing the elements within the excavation border have been balanced before activating the elements of the retaining system in order to consider the 3D effects of the advancement of excavation face in the stage of excavating the side of the lower part.After activating the elements of the retaining system and changing the stiffness matrix, the remained unbalanced nodal forces have been balanced and thus, the effects of excavation face movement and exerting forces on the elements of the lining have been modelled.The field of deformation in the surrounding soil mass and the exerting forces on the retaining system are determined when the balance is established.Hence, the stages of excavation modelling and installing the lining in the second part of the lower excavation have been divided into two phases as: Phase 6: Removing the elements within the excavation border and releasing 73% of the unbalanced nodal forces.The horizontal displacements of the ground mass that have been caused by excavating the side of the lower part and distancing of excavation face from the cross-section are 45mm maximally.In order to provide an appropriate supporting for bearing the final upper lining structure, it's necessary for the final lining structure of the wall to be constructed before excavating the side of opposite lower part.After creating the final lining of the wall, excavation is done on the opposite side like before.In this case, maximum length the excavation without lining is 4m like before and maximum of the initial lining distance from the face is considered equal to 2m.The rate of stress relaxation for simulating transition of the stress and the stress relaxation caused by approaching the excavation face to the cross-section and passing it up to a certain distance (equal to 3m on average in here) is determined equal to 75% similar to the previously applied method.Thus, the stages of modelling the excavation and installing lining in the lower excavation are finished by conducting the following phases:  The horizontal displacements of the ground mass that have been caused by excavating the second side of the lower part and distancing of the face from the cross-section, maximally are 48mm at the end of Phase 10.Selected lining in this area for conducting the Continuous Region Analysis includes 30 cm of shotcrete with a strength of 250 Kg/cm along with steel frame IPE200 in distances of 0.5 m.The next results of the computations is presented assuming the distribution of axial moment in proportion to compressive stiffness between the steel frame and shotcrete and bearing the profile up to the maximum of the bending capacity.Selected lining includes 15 cm of shotcrete with a compressive strength of 250 Kg/cm along with steel frame IPE200 in distances of 0.5 m.The next results of the computations is presented assuming the distribution of axial moment in proportion to compressive stiffness between the steel frame and shotcrete and bearing the profile up to the maximum of the bending capacity.

Shotcrete Elastic Modulus: n=10
The ratio of axial moment distribution in composite section: r1 =28%, r2 =72% The interaction curve of axial and bending moment for a shotcrete with a width of 15cm is shown in Figure 37 with two layers of mesh 6@100×100.Therefore, the selected lining for the temporary middle member of the upper part in area I, includes 15 cm shotcrete with a compressive strength of 250 Kg/cm 2 along with two layers of mesh 6@100×100 and steel frame IPE200 each 0.5 m.
Maximum of the shear force within the middle temporary lining with a width of 15 cm, is equal to 5.5 ton.Regarding the shear capacity of the steel profile, there is no problem in terms of the shear strength.

Designing the Lining of the Wall in the Lower Part
Selected lining in this area for conducting the Continuous Region Analysis includes 20 cm of shotcrete with a strength of 250 Kg/cm along with steel frame IPE160 in distances of 0.5 m.The next results of the computations is presented assuming the distribution of axial moment in proportion to compressive stiffness between the steel frame and shotcrete and bearing the profile up to the maximum of the bending capacity.

Shotcrete Elastic Modulus: n=10
The ratio of axial moment distribution in composite section: r1 =18%, r2 =82% The interaction curve of axial and bending moment for a shotcrete with a width of 20cm is shown in Figure 38 with two layers of mesh 6@100×100.Thus, the selected lining for the wall of the lower part of area I, includes 20cm shotcrete with a compressive strength of 240 Kg/cm 2 along with two layers of mesh 6@100×100 and steel frame IPE160 each 0.5 m.

Designing the Final Lining
Regarding the geomechanical properties of the soil mass and exerting forces on the initial lining, constructing the final lining of the upper part before the lower excavation is inevitable due to stabilizing the tunnel crown in the excavation and thus, as mentioned in the analysis part of this report, the final lining bears a considerable part of dead loads caused by the constructing stages and the ones after them.Hence, the forces that have been considered in constructing the final structure of the tunnel are: 1. Dead loads caused by earth lateral pressure at different stages of the construction and at the end of excavation stages.
2. Forces caused by the weight of the structure.

Bending Design
Considering the variable entity of moment in the different sections and its changing thickness at the bottom, wall, and the crown, it's necessary to set the most critical mode of force and capacity of the sections as a criterion for designing the latitudinal bending armatures.
The results of analyzing and designing latitudinal (bending) armatures for critical sections have been shown in Table 2. Figures 40 to 42 show the interaction curve of the critical sections.

Discussion
Generally, designing the underground spaces is designing the retaining systems for them.This topic forms a diverse range that includes processes from temporary mine borings that are conducted in good rocks and don't need  retaining system to big permanent underground systems in civil engineering that are consisted of retaining systems composed of full injection along with steel meshes and shotcretes.This two extents can be considered as high and low extents of designing the retaining systems in underground spaces.The necessities of designing every underground space must be based on using the rock as main materials.The amount of fracturing and weakening the surrounding rocks must be minimized in the boring process to reduce the need for steel or concrete retaining system.Usually, the hard rocks are much stronger than concrete before being under pressure stresses and also when they are in an untouched mode.Many of them can even be considered equal to steel in terms of strength.Consequently, it's not cost effective to change materials like rocks that may be strong enough with materials that their excellence is uncertain.The extent and scope of utilizing the rocks as retainers, depends on geologic conditions in the underground space and also, these depend on that to which extent the designing is sensitive to the conditions and how much it wants to include them in the design.Thus, accurate interpretation of geologic issues in the boring area is one of the necessary conditions for reasonable designing.

Conclusion
The horizontal displacements of the ground mass that have been caused by excavating the second side of the lower part and distancing of the face from the cross-section, maximally is 48mm at the end of Phase 10.The occurred deformations indicate the elastic-plastic behavior of ground mass around opening area.73% of the unbalanced nodal forces caused by removing the elements within the excavation border have been balanced before activating the elements of the retaining system.The type of soil in the tunnel area is fine-grained soils (slits and clays) which can usually be found in the beginning parts of the tunnel.The type of rocks surrounding the tunnel is Igneous and they are related to the Tuff and Rhyodacite rocks.They are gappy with cataclastic formation and aphanitic texture colored from creamy to pink.The gaps and deformed areas are full of ferruginous minerals.Compressive strength is in the range of 400 kg/cm 2 to 500 kg/cm 2 .Elastic modulus of these rocks is between 12GPa and 13GPa.Poisson's ratio is 0.35 on average.The Cohesion (C) ranged from 4-5MPa to 5MPa and the Angle of Internal Friction (φ) is between 60ᵒ and 50ᵒ.In the tunnel, by increasing the depth, the porosity and water absorption of the rocks will be reduced gradually, and also the density of the rocks will increase.In the area of the bored holes, there was no indication of water that shows the lower level of ground water in comparison to the tunnel site.Final score of the surrounding rocks in RMR system was 38 which put them on the weak rocks category.The score was 1.52 using Q classification method that also is on the weak side of the classes.Retaining system in RQD method is the utilization of extended rock bolts with distances of 0.6 to 1.2 and with steel frames that weighed from average to heavy.Suggested method of RMR classification system is to install arranged rock bolts with lengths of 4m to 5m and distances of 1m to 1.5m in the roof and walls as well as installing wire meshes and shotcrete 100mm to 150mm and 100mm on the walls, respectively.The Q classification system suggests single rock bolts with a length of 1.9m and distances of 1m to 1.5m.However, the equal dimension in this system is 9.3m for the tunnel.Maximum opening of the tunnel without retaining, obtained as 3.48 m in Tuff (Rhyodacite) rock unit.Boring the tunnel is performed using Drill-and-Blast operation and it has two stages.First stage includes designing and operating the initial structure in the upper part with the height of 5.5m and a semi-circular cross section.The second stage includes drilling in the lower part with the height of 3 m and the width of 7.34m and with a rectangle cross section.Hereafter, the stabilization will begin.At the stabilization level, the used components and materials were examined in laboratory in each stage, regarding the properties of sides and roof of the tunnel as well as the pressure on them.Then, proceedings of stabilization were done considering the results of the aforementioned actions.

Figure 2 .
Figure 2. Cross-section of the entrance and exit galleries

Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .
Figure 6.Time distribution of tunnel site earthquakes within a radius of 100km in different years

Figure 11 .
Figure 11.Epicenter and hypocenter distribution of the occurred earthquakes in Imam Reza Tunnel area within a radius of 50 km

Figure 12 .Figure 13 .
Figure 12.Computer model of the tunnel and the finite element network

Phase 1 :
Removing the elements within the excavation border and releasing 70% of the unbalanced nodal forces.Phase 2: Activating the elements of the lining and releasing the remained unbalanced nodal forces until the establishment of full balance.

Figure 14 .Figure 15 .
Figure 14.Mean of the main stresses in excavation of the first upper part, before installing the lining (End of Phase 1)

Figure 16 .Figure 17 .Figure 18 .
Figure 16.Ground mass displacement after activating the lining and advancement of the excavation face in the first upper part (End of Phase 2)

Figure 19 : 3 )Figure 20 .
Figure 19: Mean of the main stresses in the completion of excavating the upper part and before installing the lining (End if Phase 3)

Figure 21 .Figure 22 .Figure 23 .
Figure 21.Ground mass displacement after activating the lining and advancement of the excavation face in the upper part (End of Phase 4)

Phase 5 :
Activating the elements of final lining in the upper part, removing the elements within the lower excavation border (middle part), and releasing all of the unbalanced nodal forces until the establishment of full balance (Figures24 and 26).

Figure 24 .Figure 25 .Figure 26 .
Figure 24.Ground mass displacement after activating the final lining of the upper part and advancement of the excavation face in the middle of the lower part (End of Phase 5)

Phase 7 :
Activating the elements of the lining and releasing the remained unbalanced nodal forces until the establishment of full balance.

Figure 27 .Figure 28 .
Figure 27.Mean of the main stresses in excavating the side of the lower part before installing the lining (End of Phase 6)

Figure 29 . 7 )Figure 30 .
Figure 29.horizontal displacements of the ground mass that have been caused by excavating the side of the lower part (End of Phase 7)

Phase 8 :
Activating the elements of the final lining of the wall on one side.Phase 9: Removing the elements within the side excavation border and releasing 75% of the unbalanced nodal forces.Phase 10: Activating the elements of the initial lining of the wall and releasing the remained unbalanced nodal forces until the establishment of full balance.Phase 11: Activating the elements of the final lining of the wall and completing the final lining.

Figure 31 . 9 )Figure 32 .
Figure 31.Mean of the main stresses in excavating the second side of the lower part before installing the initial lining (End of Phase 9)

Figure 33 .Figure 34 . 10 Figure 35 . 3 . 5 . 3 . 5 . 1 .
Figure 33.The horizontal displacements of the ground mass that have been caused by excavating the side of the lower part and distancing the face from the cross-section up to ∞ after installing the initial lining (End of Phase 10)

Figure 36 . 3 . 5 . 1 . 1 .
Figure 36.Interaction curve of axial and bending moment for a shotcrete with a width of 30cm with two layers of mesh 6@100×100

Figure 37 .
Figure 37. Interaction curve of axial and bending moment for a shotcrete with a width of 15cm with two layers of mesh 6@100×100

Figure 38 .
Figure 38.Interaction curve of axial and bending moment for a shotcrete with a width of 20cm with two layers of mesh 6@100×100

Figure 39 . 1 Figure 41 . 2 Figure 42 .
Figure 39.Critical sections of designing the final lining in area I