Simplified Full Dynamic Analysis of a Railway Tunnel in Ethiopia ()
1. Introduction
Analysis and design of tunnel linings are often performed for static loading cases only [1]. This is mainly because fewer damages were recorded historically for underground structures than superstructures subjected to both static and dynamic loading cases. Deep seated tunnels are thus assumed to be seismic resistant. However, there will be high possibility of damage if the tunnel is located close to an earthquake fault or under difficult geological conditions.
[2] [3] [4] [5] reported cases of damages of tunnels due to seismic loading recorded in the past. Lessons learned from the damages of such underground structures indicated the significance of performing dynamic analyses for important structures, especially in earthquake prone regions. Such analyses also provide the possibility of assessing building response analyses, as tunnels in earthquake prone regions influence overlying structures [6].
Shortage of pertinent seismic records, especially in developing countries like Ethiopia, is also among the factors that such structures are designed for static load only. This research is aimed at generating such data for dynamic analysis of a tunnel in one of the seismic zones of Ethiopia.
Kinde et al. [7] cited Gouin [8] for reporting the existence of about 15,000 tremors, strong enough to be felt by humans that had occurred in Ethiopia and the Horn of Africa in the 20th century alone. A total of 16 earthquakes of magnitude 6.5 and higher were recorded in some of Ethiopia’s seismic active areas during that century. In 1961, an earthquake of 6.7 magnitude was experienced in Karakore, located in the escarpment seismic zone of Ethiopia. That was one of the largest earthquake records in the country. The railway tunnel at Karakore considered in this research is part of the 389 km line stretching from Awash to Haragebeya. It is among the new railway projects in the country, constructed by the Turkish Contractor, Yapi Merkezi (Figure 1).
Among the available methods of seismic analysis of tunnels, full dynamical analysis employing numerical methods are found to be more effective [4]. Even though three-dimensional analysis may simulate the actual conditions in a better way, the plain strain analysis is found to deliver results which are in good agreement with the 3D modeling [9]. To minimize the enormous computational time required for the 3D computations, this research followed the 2D plain strain analysis approach. Appropriate techniques have been followed for modeling some components of the lining that are difficult to be modeled two dimensionally.
Since Earthquakes of magnitude greater than 5.0 generate sufficiently large ground motions that are potentially damaging to structures [10], it is essential to perform full dynamic analysis for important structures like the tunnel at Karakore. To assess safety at the worst conditions, the tunnel section at the weakest ground condition and maximum overburden of 320 m (between stations 178 + 807 and 178 + 864) has been selected for the full dynamic analysis. The basic input for the dynamic analysis of the tunnel is the design earthquake, determined by performing seismic hazard analysis, which is considered in the next sections.
2. Seismic Hazard Analysis
The railway alignment under consideration is located in an area identified by the Ethiopian Geological Survey as “Were-Ilu” area, located between 10˚N - 11˚N
![]()
Figure 1. Location of the specific site from the railway network in Ethiopia (Source: https://www.skyscrapercity.com/).
and 39˚E - 40˚30'E, which is part of the escarpment seismic zone of Ethiopia [11]. Since acceleration time history records required for full dynamic analysis are not available in Ethiopia, both probabilistic seismic hazard analysis and deaggregation have been conducted to filter acceleration time history records from an existing earthquake database.
2.1. Probabilistic Seismic Hazard Analysis
Identification of the design earthquake for tunnel safety analysis is performed by conducting probabilistic seismic hazard analysis (PSHA). In PSHA, past seismicity is assumed to be suitable for prediction of future seismicity [12] [13]. Among the two levels of earthquake ground motions commonly considered in tunnel design and analysis [14]; namely, Maximum Design Earthquake (MDE) and Operating Design Earthquake (ODE), emphasis has been given to ODE due to the shorter return period as compared to MDE.
2.1.1. Earthquake Source, Magnitude and Distance Identification
Earthquake sources identified by the Institute of Geophysics, Space Science and Astronomy (IGSSA) in the region defined by 9.40˚ - 11.48˚ (latitude) and 38.93˚ - 41.10˚ (longitude) have been used for this study. These point sources, which are 74 in number, indicate that initiation point for earthquake around the selected site.
Earthquake magnitudes are distributed according to the well-known Bounded Guttenberg-Richter recurrence law [15] given by Equation (1).
(1)
where λm is the rate of earthquakes with magnitudes greater than m;
a and b are constants that depend on the earthquake magnitudes.
Depending on the contemporary earthquake data of IGSSA, the values of the constants a, b and m are taken to be 1.318, 0.93 and 6.86 respectively, as recommended by Ayele [16].
The distribution of distance between earthquake source and the tunnel was modeled using the 74 contemporary data having specific locations and source. Using the points along the tunnel alignment, where survey data were taken in combination with seismic source coordinates, the hypocentral distance was calculated by applying Pythagorean Theorem from the epicentral and vertical distances. Each survey point along the tunnel alignment was taken into account for each source with coordinates for distance determination and the point which gives the shortest distance is taken as the hypocentral distance for that specified source.
2.1.2. Ground Motion Intensity
Prediction of ground motion is the major theme of probabilistic seismic hazard analysis. Ground motion prediction models that are developed through statistical regression from observations are used to predict probability distribution of ground motion intensity. Among the three Ground Motion Prediction Equations (GMPE) for global earthquake in active shallow regions [17], the one whose input parameters agree with the available data [15] has been chosen for prediction purpose (Equation (2)).
(2)
Y = ground motion parameter,
FM(M) = magnitude scaling,
FD(RJB,M) = distance function,
Fs(Vs30,RJB,M) = site amplification term,
ε = fractional number of standard dev. of single predicted values of lnY,
σT = intra- and inter-event uncertainty coefficient.
All these terms have been determined by using the formulations provided by Boore and Atkinson [18] depending on the prevailing site conditions.
The ground motion analysis has been carried out transforming the actual layered ground into an equivalent single layer profile with suitable parameters. This single layer was used since the selected ground motion prediction equation can’t be applied for the shear wave velocity of the top 30 m, Vs30, of the layered soil.
The geological structure of the area is characterized by normal, strike and reverse slip faults [11]. Even if normal slip faulting mechanism is found to be the dominant one, peak ground acceleration (PGA) values of strike slip fault type have been used for the analyses to account for the worst case scenario.
Lognormal distribution is used to model many engineering data [19]. Specifically, peak ground acceleration distribution is characterized by lognormal distribution, provided that the annual probability of exceedance, P (PGA), is not less than 10−6 [20]. Accordingly, distribution of peak ground acceleration has been modeled by lognormal distribution for PGA values ranging from 0.05g to 0.8g with an increment of 0.05g. Annual exceedance rate of PGA, λ (PGA > x), for respective PGA values have thus been calculated to produce the hazard curve depicted in Figure 2.
2.2. Deaggregation
The determination of PGA values having a certain percent probability of exceedance based on aggregation of contributions from potential earthquake magnitudes and source to site distances was shown for seismic hazard analysis. Combinations of magnitudes and distances have been used so that this effect cannot be dedicated to a single or few pairs of magnitude and distance values.
According to the Federal Highway Administration of the U.S (FHWA), 10% and 2% probability of exceedance in 50 years are used to describe design seismic actions for damage limitation and no collapse requirement of MDE and ODE respectively [21]. These are equivalent to return periods of 475 and 2475 years, which have been reciprocated to result in PGA with a probability of exceedance of 0.0021 and 0.0004 for ODE and MDE respectively. Interpolation of the values
![]()
Figure 2. Hazard curve generated for the analyses corresponding to stations 178 + 807 - 178 + 864.
from the hazard curve gives a PGA of 0.36g for operating design earthquake (ODE).
Similar calculations of PGA have been performed by considering four different ground conditions and varying the top 30 m shear wave velocities, with the aim of addressing the worst effects. The considered ground conditions are:
· A homogeneous ground having the properties of the middle soil layer.
· Multi layered rock with intact properties.
· Multi layered rock with weathered properties.
· Single layer below the tunnel.
The summary of the major inputs used in each of these cases together with the corresponding outputs is demonstrated in Table 1.
Among the different ground conditions under consideration, the assumption of a single weathered homogeneous layer of appropriate material properties (shear modulus, G and density, ρ) resulted in the maximum PGA value of 0.36 g. No PGA calculation was performed for the case of multi-layered weathered rock since the value of the shear wave velocity deviates from the threshold ranges of Vs30 (between 180 m/s minimum and 1300 m/s maximum) to be used in GMPE of Boore and Atkinson [18].
According to [22], ground motion analysis often requires a single design earthquake where the earthquake threat is characterized by a single magnitude and distance. Thenhaus et al. [23] stated that earthquakes that contribute the most to the total hazard result from the discretized interval with the largest relative contributions. Using the deaggregation process, the most likely earthquake scenario that causes the peak ground acceleration with the specified probability of exceedance can be identified. Accordingly, 2D hazard deaggregation has been conducted for the PGA value of 0.36g, using Magnitude-Radius (M-R) distribution. The deaggregation calculations resulted in two modes of ODE; one with 5.6 magnitude & 20 km radius; and the other with 5.8 magnitude & approximately 33 km radius. This pair of values is considered as the most likely scenario which
![]()
Table 1. Ground conditions with the respective parameters andcorresponding calculated PGAs.
a. Not applicable as Vs30 is not within the threshold for GMPE by Boore and Atkinson.
corresponds to a 10% probability of exceedance in 50 years for the specific site being investigated.
2.3. Ground Motion Selection
The selection of appropriate ground motion is the basis for defining seismic load for dynamic analysis of a structure. Since no records of accelerograms exist for Ethiopian earthquakes, other records should be adapted for the dynamic analysis. The required ground motion for the dynamic loading is selected by using magnitude and distance values obtained from the deaggregation calculation [24]. Natural time histories from the database of Pacific Earthquake Engineering Research Center (PEER) were used by taking the first pair of the most likely magnitude and radius from the deaggregation calculation, shear wave velocity range, and fault mechanism. Accordingly, existing ground motion records have been selected. Ground motion with the normal oblique faulting mechanism has been entered into the PEER database search machine, due to the dominancy of these structures in the project area. The input parameters thus used as an input for the PEER database are shown in Figure 3.
The ground motion called L’Aquila Italy recorded at station Avezzano in 2009 with a magnitude of 5.6 was selected based on the above inputs. It was characterized by a Vs30 value of 199 m/s and the closest distance to the surface projection of the oblique normal fault, Rjb, of 27.38 km. This selected ground motion with a durationof 16.3 seconds has further been scaled by considering the effects of depth attenuation, PGA, and loads pertinent to the tunnel structure under consideration.
The structure is located 320 m below the surface indicating the need for depth attenuation correction of ground motion parameters. Accordingly, a depth attenuation factor of 0.7 has been taken as recommended by the Technical manual
![]()
Figure 3. Search parameters for PEER database.
for the design and construction of road tunnels provided by FHWA [21].
The PGA of the selected ground motion, which is equal to 0.03 g, shall be scaled to the PGA defined by the PSHA in the previous section, 0.36 g; leading to a scaling factor of:
0.36 g/0.03 g = 11.9.
Based on the recommended seismic loading for operating design earthquake, ODE, the load correction is determined using the load factor design-based formulation of Hashasha et al. [25] for circular tunnel lining. Accordingly, an earthquake load effect factor of 1.3 has been considered to account for the seismic load effects.
The selected ground motion has thus been scaled using the total factor K, determined by combining all the above three scaling factors as:
K = 0.7 × 9.81 × 11.9 × 1.3 = 106.23.
The corresponding time acceleration history, after being scaled with this factor, is presented in Figure 4.
3. Numerical Modeling
With the aim of considering the worst effects from the given conditions, a combination of the maximum overburden and the poorest ground conditions has been chosen for the numerical analyses. This very poor section is located between stations 178 + 807 and 178 + 864. Initial tunnel support was mainly provided by shotcrete, while rock reinforcement, lattice girders, and pipe umbrella were used at different sections in addition to the shotcrete.
The seismic response of the tunnel has been analyzed by using the Finite Element software PLAXIS 2D. The difficulty of modeling the diagonally extruding pipe umbrella was facilitated by considering its stiffening effect with assignment of an equivalent stiffness of the soil-umbrella composite. Among the three types of 2D continuum methods of numerical analysis recommended by [21], dynamic time history analysis which considers mechanically coupled soil and tunnel responses has been used by applying the ground motion time histories at the base of the numerical model as suggested by [26]. The standard earthquake boundary conditions with absorbent boundaries were used at vertical boundaries while the bottom boundary was restrained vertically and given a prescribed displacement of 1 m in the horizontal direction as depicted in the model geometry of Figure 5.
![]()
Figure 4. Scaled time acceleration history.
The soil-structure interaction has been modeled by employing the no-slip interface condition recommended by Pescara et al. [27]. This assumption will maximize the thrust forces acting on the lining [25] [28], which in turn makes the results more conservative. Linear elastic properties were used for modeling the lining while the Mohr-Coulomb constitutive model was employed for the soil (including the ballast and invert). Geotechnical parameters have been adopted from the data provided by the Contractor (Yapi Merkezi). The pipe umbrella, rock bolts, and shotcrete were considered as parts of the primary linings. The material parameters of the equivalent soil cluster at the crown,where the pipe umbrella is provided, were determinedby converting to equivalent values based on the proportions of the soil and umbrella cross sections. Components of the lining including primary and secondary supports, together with their parameters are summarized in Table 2.
The actual construction procedure of the tunnel that was carried out byemploying the New Austrian Tunneling Method (NATM) was simulated properly by using the staged construction scheme of PLAXIS. Since the construction sequence has critical effect on ground deformations [1], the following eleven calculation phases were followed for the realistic simulation purpose.
Phase 1: Initial stress condition,
Phase 2: Pipe umbrella construction,
Phase 3: Top heading construction,
Phase 4: Construction of lining at the crown & provisional invert,
Phase 5: Bench excavation,
Phase 6: Full primary lining,
Phase 7: Invert construction,
Phase 8: Final lining,
Phase 9: Ballast construction,
Phase 10: Sleeper construction,
Phase 11: Dynamic load application.
4. Discussions on the Critical Numerical Analyses Results
The numerical computations of the separate static and dynamic analyses indicated that the resulting stresses have extreme differences. This is attributed to the strong motion duration and reduced support capacity of the assumed unreinforced section used for the analyses, in addition to the maximum PGA value of 0.36 g used as input for the dynamic analyses. The consideration of the unreinforced concrete section together with the ignored lining components, i.e., fiber reinforcement, wire mesh or lattice girders account for the extremely reduced section capacity.
Since the extreme stresses correspond to the combined effects of the static and dynamic loads, more focus has been given to the stresses derived from that critical load combination. The maximum values of the axial force and bending moment shown in Figure 6 are used for the detailed assessment of the section capacity. Both maximum values are found to appear at the crown-side wall interface, where stiffness difference manifests especially due to the pipe umbrella and the rock bolts.
The combined effect of the axial force and bending moment induced a tensile stress in excess of the section capacity of the plain section, about 10 times larger than the unreinforced concrete capacity that was determined according to the Ethiopian standard [29]. The steel reinforcements ignored in the computational model, which are mainly expected to serve the purpose of resisting tensile stresses, will be expected to account for a major portion of this deficiency. The capacities
![]()
Figure 6. Axial force and bending moment diagrams on the lining for the critical load case.
of fiber reinforcement, wire mesh, lattice girders were also not considered in the computations, even though they were constituents of the tunnel support. On the other hand, the plane strain analyses do not account for the stress release associated with the deformations ahead of tunnel face. According to the British Tunneling Society and The Institution of Civil Engineers [30], 30% - 50% of the deformation during construction occurs ahead of the face which will result in significantly reduced pressures acting on the lining than predicted from a two- dimensional analysis. The conservative no-slip assumption of the interaction at the tunnel soil interface has also contributions to the increased stresses.
The superimposed effect of these factors, together with the additional safety embedded within the strong motion duration used for the determination of the design earthquake; will ensure the safety of the tunnel [31]. No further rigorous modeling with three dimensional numerical analyses are recommended due to the associated enormous computational efforts. It is thus demonstrated how to prove the safety of a tunnel in an earthquake prone region with limited input data. The minimum possible computational efforts were used to perform the seismic hazard analysis as well as the Finite Element Analysis by employing efficient and satisfactory simplification schemes.
5. Conclusions
This research demonstrates a professional method of carrying out dynamic analysis of a tunnel in regions with scarcity of seismic data by making rational assumptions in the determination of the loads as well as section capacity. The poorest rock mass condition from the longitudinal stretch of a tunnel at Karakore, located in the escarpment seismic zone of Ethiopia, was considered to represent worst case conditions. It was made possible to prove the safety of the tunnel by making reasonable simplifications of the complex reality with the minimum possible efforts and computational time.
Determination of the earthquake loading on the tunnel structure has been carried out by making appropriate conservative assumptions at different calculation stages. The seismic hazard analysis under operating design earthquake gives rise to a maximum PGA value of 0.36g as governed by the assumption of a single soil layer around the tunnel. Deaggregation process resulted in two pairs of design earthquakes: the first with magnitude 5.6 and radius 20 km and the second with magnitude of 5.8 and radius 33 km. The design earthquake has been used to select ground motion from the PEER ground motion database, which in-turn was scaled to incorporate specific site and loading conditions.
With the aid of the carefully simulated plane strain computations using the Finite Element software PLAXIS 2D, it was possible to prove the stability of the tunnel at the extreme load combination, with minimum computational efforts. The capacity of plain concrete alone was considered while calculating the strength of the lining. Even though the conservatively assumed plain concrete section was not capable of resisting the calculated tensile stress obtained for the worst combination of load and resistance, the reserves left in parameters considered and the resistance of the ignored lining components, together with the practical three dimensional stress release effects will potentially bridge over the capacity gaps. This is also supported by the findings of Jaramillo [5] stating that tunnels constructed in rocks with peak ground acceleration of 0.5g will not be susceptible to damages due to dynamic actions.
The simplified numerical approach followed in this research gives the opportunity of analyzing tunnel safety with minimum computational efforts. This procedure can be used especially in areas with limitations in input data and computational aid (software). It will also be pretty much appropriate for preliminary design of tunnels on difficult ground.
Acknowledgements
The Authors would like to thank the Turkish Contractor Yapi Merkezi, especially the project manager Mr. Abdullah Kilic, for facilitating site visits and availing the necessary documents required for the purpose of this research.