Analytical Method for Preliminary Seismic Design of Tunnels

Abstract

Buried structures are categorized based on their shape, size and location. These main categories are near surface structures (e.g., pipes and other facilities), large section structures (e.g., tunnels, subways, etc.), and vertical underground structures (e.g., shafts and ducts). Seismic assessments of these structures are important in areas close to severe seismic sources. Seismic design of tunnels requires calculation of the deformation in surrounding geological formations. The seismic hazard on a site is usually expressed as a function of amplitude parameters of free-field motion. Therefore, simplified relations between depth and parameters of ground motion are necessary for preliminary designs. The objective of this chapter is to study and review the main analytical seismic methods which are used to develop a simple relationship between maximum shear strain, maximum shear stress and other seismic parameters.

1. Introduction

A seismic ground motion poses a threat to urban infrastructure as well as human life. Individuals have a limited understanding of underground structures’ seismic resistance. Because of smaller deformations under the condition of encompassing rock or soil constraints, it is widely agreed that an underground structure is much more stable than a ground structure. Several communities have emerged in the United States of America to explain seismic behavior of underground opening under severe conditions since the 1990s. Numerous destructive seismic events, such as the Kobe, Chi-Chi, Kocaeli and Wenchuan earthquakes, have occurred since the 1990s, causing genuine harm to tram stations and tunnels, indicating that underground structures are still vulnerable to damage under intense seismic motions. A characteristic example of broad damage due to ground shaking and permanent displacements is the Hanshin earthquake caused liquefaction that contributed to the collapse of numerous underground structures in 1995, counting a tram station in Kobe, Japan, damages to highway tunnels during 1999 Chi- Chi and the collapse of the twin Bolu under construction tunnels, during the 1999 Kocaeli earthquake [1].

Owen and Scholl [2] characterized the deformation sorts of underground structures due to seismic excitation as axial compression/extension; longitudinal bending, ovaling, and racking deformations (Figure 1). Shear deformation of tunnels initiated by the vertically propagating shear waves has been broadly investigated by a number of researchers [3, 4], and it has been demonstrated to be the basic mode of deformation for tunnels under seismic loading. Ovaling and racking deformations are related to normal or nearly normal propagation of shear waves with respect to tunnel axes which cause distortion of tunnel cross section. Simplified seismic design approaches for tunnels are often favored by experts. They should be able to assess the general response of a tunnel system that has been subjected to seismic loading. As a result, simpler methods for measuring maximum shear strain (γ_max) in the tunnel depth are used [1].

Figure 1.

Types of deformations on tunnels under seismic actions (a) compression extension, (b) longitudinal bending deformation, (c) compression of tunnel section [2].

Many researchers proposed analytical solutions to estimate the seismic internal forces of tunnel linings under certain assumptions and conditions, such as elastic response of the soil and tunnel lining, and seismic loading simulation in semi-static construction, among others. Analytical solutions are useful, moderately fast, and easy to use for fundamental seismic design of tunnels, despite the fact that they are formed using relatively strict assumptions and simplifications. As a result, they’re commonly used in the early stages of design. With the improvement in technology and computer science, and consequently in numerical analysis of material deformation and stability, several methods are used for analysis of underground structures such as finite element, finite difference and discrete element method. Analyzing of axial and bending deformations can be best performed using 3-D models. In finite difference or finite element models, the tunnel is discretized spatially and the surrounding soil is either discretized or models by springs. Several computer codes perform these type of analysis such as FLAC, ABAQUS and so on [1].

2. Simplified estimation of ground deformations

The seismic design of tunnels is based on two approaches: (1) soil-structure interaction and (2) free field approach. In the first approach, the soil shear strains are affected by the deformation of the nearby underground structures and will conform to the structure strains. A reduction in the total mass of the soil and structure at the soil cavity may have a significant effect on the shear strain. In this case, shear strain of soil in the vicinity of structure will be greater than the free-field approach. In the free-field approach, the interaction between soil and structure is neglected and it is expected that structures accommodate the forced deformations from encompassing ground. These deformations are a function of maximum shear strain [1, 5]. The direct measurement of strains is not possible so it is correlated to other strong-motion parameters such as Peak Ground Velocity (PGV) [6, 7]. Newmark considered one-directional propagation of the harmonic wave in a homogeneous, isotropic, and elastic unbounded medium. According to Newmark, relationship between the maximum particle velocity (V_max) and (γ_max) is.

$${\mathrm{\gamma }}_{\max }={\mathrm{V}}_{\max }/\mathrm{C}$$E1

Where C is the apparent wave velocity [8].

C cannot be estimated straightforward and is depended on wave type, the angle of incidence, and material property [9]. To calculate this parameter, some formulas are proposed. For instance, O’Rourke and Elhmadi [10] proposed a relation for calculation of longitudinal deformation on buried pipes:

$$\mathrm{C}=\mathrm{Vs}/\text{sinØ}$$E2

Where Ø is angle of the incidence at the ground surface and Vs is the shear wave velocity of the top layer. C is variant at different geological situations [10, 11, 12]. Ovaling and racking deformations are correlated with γ_max on a vertical plane, so C is close to C_s, which is the incident horizontal shear-wave velocity in geological layers. The consequent structural deformations are basically related to γ_max in the imperforated ground as shown in Figure 2 [13, 14, 15].

Figure 2.

Ovaling and racking deformation on buried structures [5].

Wang [13] considering ovaling deformation related C to effective shear modulus, G, and the mass density of the medium, ρ by.

$$\mathrm{C}=\surd \mathrm{G}/\rho$$E3

In the case of replacement of Eq. (3) in Eq. (1) some problems may arise such as the indeterminacy in the definition of deep depth or application of this formula for layered strata. Considering all these issues, they are still adopted by most of the available technical guidelines [6, 7, 12].

St. John and Zahrah [9] developed Newmark’s formula and proposed relationships to estimate longitudinal, normal and shear strains in the free field which is depicted in Table 1.

Wave Type		Axial Strain	Shear Strain	Curvature
P-wave		$\mathrm{\epsilon }=\frac{{\mathrm{V}}_{\mathrm{P}}}{{\mathrm{C}}_{\mathrm{P}}}{cos}^2\mathrm{\varphi }$	$\mathrm{\gamma }=\frac{{\mathrm{V}}_{\mathrm{P}}}{{\mathrm{C}}_{\mathrm{P}}}sin\mathrm{\varphi }cos\mathrm{\varphi }$	$\frac{1}{\mathrm{\rho }}=\frac{{\mathrm{a}}_{\mathrm{P}}}{{\mathrm{C}}_{\mathrm{P}}^2}sin\mathrm{\varphi }{cos}^2\mathrm{\varphi }$
		${\mathrm{\epsilon }}_{\max }=\frac{{\mathrm{V}}_{\mathrm{P}}}{{\mathrm{C}}_{\mathrm{P}}}\text{for}\mathrm{\varphi }=0^{\circ }$	${\mathrm{\gamma }}_{\max }=\frac{{\mathrm{V}}_{\mathrm{P}}}{2{\mathrm{C}}_{\mathrm{P}}}\text{for}\mathrm{\varphi }={45}^{\circ }$	$\frac{1}{{\mathrm{\rho }}_{\max }}=0.385\frac{{\mathrm{a}}_{\mathrm{P}}}{{\mathrm{C}}_{\mathrm{P}}^2}\text{for}\mathrm{\varphi }={35.27}^{\circ }$
S-wave		$\mathrm{\epsilon }=\frac{{\mathrm{V}}_{\mathrm{S}}}{{\mathrm{C}}_{\mathrm{S}}}sin\mathrm{\varphi }cos\mathrm{\varphi }$	$\mathrm{\gamma }=\frac{{\mathrm{V}}_{\mathrm{S}}}{{\mathrm{C}}_{\mathrm{S}}}{cos}^2\mathrm{\varphi }$	$\frac{1}{\mathrm{\rho }}=\frac{{\mathrm{a}}_{\mathrm{S}}}{{\mathrm{C}}_{\mathrm{S}}^2}{cos}^3\mathrm{\varphi }$
		${\mathrm{\epsilon }}_{\max }=\frac{{\mathrm{V}}_{\mathrm{S}}}{2{\mathrm{C}}_{\mathrm{S}}}\text{for}\mathrm{\varphi }={45}^{\circ }$	${\mathrm{\gamma }}_{\max }=\frac{{\mathrm{V}}_{\mathrm{S}}}{{\mathrm{C}}_{\mathrm{S}}}\text{for}\mathrm{\varphi }=0^{\circ }$	$\frac{1}{{\mathrm{\rho }}_{\max }}=\frac{{\mathrm{a}}_{\mathrm{S}}}{{\mathrm{C}}_{\mathrm{S}}^2}\text{for}\mathrm{\varphi }=0^{\circ }$
R-wave	Compressional Component	$\mathrm{\epsilon }=\frac{{\mathrm{V}}_{\mathrm{R}}}{{\mathrm{C}}_{\mathrm{R}}}{cos}^2\mathrm{\varphi }$	$\mathrm{\gamma }=\frac{{\mathrm{V}}_{\mathrm{R}}}{{\mathrm{C}}_{\mathrm{R}}}sin\mathrm{\varphi }cos\mathrm{\varphi }$	$\frac{1}{\mathrm{\rho }}=\frac{{\mathrm{a}}_{\mathrm{R}}}{{\mathrm{C}}_{\mathrm{R}}^2}sin\mathrm{\varphi }{cos}^2\mathrm{\varphi }$
		${\mathrm{\epsilon }}_{\max }=\frac{{\mathrm{V}}_{\mathrm{R}}}{{\mathrm{C}}_{\mathrm{R}}}\text{for}\mathrm{\varphi }=0^{\circ }$	${\mathrm{\gamma }}_{\max }=\frac{{\mathrm{V}}_{\mathrm{R}}}{2{\mathrm{C}}_{\mathrm{R}}}\text{for}\mathrm{\varphi }={45}^{\circ }$	$\frac{1}{{\mathrm{\rho }}_{\max }}=0.385\frac{{\mathrm{a}}_{\mathrm{R}}}{{\mathrm{C}}_{\mathrm{R}}^2}\text{for}\mathrm{\varphi }={35.27}^{\circ }$
	Shear Component		$\mathrm{\gamma }=\frac{{\mathrm{V}}_{\mathrm{R}}}{2{\mathrm{C}}_{\mathrm{R}}}cos\mathrm{\varphi }$	$\frac{1}{\mathrm{\rho }}=\frac{{\mathrm{a}}_{\mathrm{R}}}{{\mathrm{C}}_{\mathrm{R}}^2}{cos}^2\mathrm{\varphi }$
			${\mathrm{\gamma }}_{\max }=\frac{{\mathrm{V}}_{\mathrm{R}}}{{\mathrm{C}}_{\mathrm{R}}}\text{for}\mathrm{\varphi }=0^{\circ }$	$\frac{1}{{\mathrm{\rho }}_{\max }}=\frac{{\mathrm{a}}_{\mathrm{R}}}{{\mathrm{C}}_{\mathrm{R}}^2}\mathrm{\varphi }=0^{\circ }$
where:
VP = soil particle velocity caused by P-waves
aP = soil particle acceleration caused by P-waves
CP = apparent propagation velocity of P-waves
VS = soil particle velocity caused by S-waves
aS = soil particle acceleration caused by S-waves
CS = apparent propagation velocity of S-waves
VR = soil particle velocity caused by R-waves
aR = soil particle acceleration caused by R-waves
CR = propagation velocity of R-waves
1/p = curvature

Table 1.

Strain and curvature due to body and surface waves [9].

If the shear waves propagate vertically in a uniformly elastic half space, γ_max for a specific ground motion is a function of d/Vs, the ratio of depth below free boundary to shear-wave velocity in medium [16]. In layered medium, the equivalent travel-time concept proposed by Imai et al. [17] for estimation of maximum shear-stress (τ_max) may be used. Consequently, γ_max, can be calculated by dividing τ_max by the secant shear modulus of material G_sec, representing the average stiffness in a range of shear strain.

$${\mathrm{\gamma }}_{\max }={\mathrm{\tau }}_{\max }/{\mathrm{G}}_{\max }$$E4

For calculation of ovaling deformation, v_max is frequently assumed to be equal to the Peak Ground Velocity (PGV) in free field [10, 18]. A reduction coefficient (r_d) is proposed to reduce the ratio of ground motion at tunnel depth to motion at ground surface as it is shown in Table 2. This correlation is based on earthquake databases gathered from accelerograms [6, 7].

Tunnel Depth (m)	Ratio of Ground Motion at Tunnel Depth to Motion at Ground Surface (rd)
≤ 6	1.0
6 to 15	0.9
15 to 30	0.8
> 30	0.7

Table 2.

Ratios of ground motion at tunnel depth to motion at ground surface [6, 7].

For tunnels with shallow burial depths, maximum shear stress can be estimated by the product of Peak Ground Acceleration (PGA) in ground surface and overburden pressure [7]. This product is corrected by an empirical depth-reduction factor (r_d) due to the deformability of medium [19]. In this method, maximum shear stress (on a horizontal plane) at depth d is.

$${\tau }_{\mathit{max}=}\mathit{PGA}.\mathrm{\rho }.d.r_d$$E5

such that ρ is the density of the shallow geological formation, and d is the depth of interest. Then, maximum can be estimated by Eq. (3).

Penzien [20] also suggested closed-form solutions for seismic analysis of deep rectangular and circular tunnels, with the seismic loading being better replicated as a uniform shear-strain dissemination, τ_ff, forced on the soil boundaries of the soil-tunnel system, away from the tunnel. Penzien’s solutions, on the other hand, ignore the impact of typical stresses generated during loading along the soil-tunnel interface. They decided that the deformation of the tunnel could be approximated by the deformations of a circular cavity (e.g. through significant consideration of parameter β in Figure 3). Huo et al. [21] proposed improved arrangements by considering the genuine deformation example of rectangular-molded cavities and representing both the ordinary and shear stresses at the the soil-tunnel interface.

Figure 3.

Deformation of W × H rectangular cavity subjected to a uniform shear strain distribution γ_ff: (a) with free-field shear stress distribution applied to cavity surface; (b) with free-field shear stress distribution removed from cavity surface [20] (G: soil shear modulus, γc: shear distortion of cavity without the application of shear stress distribution around the cavity, β = γ_c/γ_ff).

Analytical solutions usually presume that the soil has a linear elastic behavior and therefore do not take into account the strain-dependent soil shear modulus. Bobet et al. [22] compensated for the reduction in shear modulus by iteratively adjusting the soil shear modulus as a function of shear strain magnitude before shear strain convergence was achieved. The analytical solution was then used to estimate the soil deformation using the compatible shear strain shear modulus [21]. The effect of soil saturation was overlooked in the production of all of the above closed-form solutions. Bobet [4] suggested circular tunnel solutions in saturated soil, assuming a non-slip interface. Bobet [23] went on to extend the previous solutions to look at the response of rectangular tunnels under no-slip and fully-slip interface conditions, as well as drained and undrained soil conditions. Park et al. [24, 25] looked over the previous solutions and proposed a new approach for considering future sliding along the soil-tunnel interface. The majority of the above-mentioned suggested analytical relationships are for shear S-waves propagating upward in the tunnel’s transverse direction. Kouretzis et al. [26, 27, 28, 29] proposed a set of relations for compressional P-wave tunnels as well.

The assumptions on which the analytical solutions are based limit their applicability (Table 3). Researchers started comparing the results of analytical solutions with the predictions of sophisticated numerical models after the rapid growth of computational power in the last two decades to recognize the shortcomings of these analytical solutions. For example, Kontoe et al. [15] compared four different analytical models (i.e. [13, 20, 23, 24]) and validated them against finite element simulations (FE). Tsinidis et al. [33] compared the results of analytical solutions (i.e. [13, 20, 24]) with numerical predictions for extreme lining flexibilities, i.e. very flexible or very rigid tunnels compared to the surrounding soil. Kontoe et al. [14] and Tsinidis et al. [33] found that the analytical solution of Penzien [20] underestimates the thrust added to the tunnel structure for a slip-free interface, which is consistent with previous findings [34]. As a result, using this solution for a rough soil-lining interface is not recommended.

Solution	Tunnel lining	Soil type	Saturation conditions	Soil layering	Soil-tunnel interface			Cross-section
	Elastic	Elastic	Dry	Homogeneous	No slip	Frictional Slip	Full Slip	Circular
St.John C.M. and Zahrah T.F [9]	Yes	Yes	Yes	Yes	Yes	No	Yes	Yes
Wang, J.N., [13]	Yes	Yes	Yes	Yes	Yes	No	Yes	Yes
Penzien and Wu [31]	Yes	Yes	Yes	Yes	Yes	No	Yes	Yes
Penzien [20]	Yes	Yes	Yes	Yes	Yes	No	Yes	Yes
Bobet [4]	Yes	Yes	Yes	Yes	Yes	No	Yes	Yes
Hou, et al. [21]	Yes	Yes	Yes	Yes	Yes	No	Yes	No
Park et al. [25]	Yes	Yes	Yes	Yes	Yes	Yes	Yes	Yes
Bobet [32]	Yes	Yes	Yes	Yes	Yes	No	Yes	Yes
Kouretzis [27]	Yes	Yes	Yes	Yes	Yes	No	Yes	Yes
Kouretzis [28]	Yes	Yes	Yes	Yes	No	No	Yes	Yes
Kouretzis [29]	Yes	Yes	Yes	Yes	No	No	Yes	Yes

Table 3.

Summary of assumptions and applicability of analytical solutions for the analysis of tunnels under ground shaking [30].

Since the soil response is often assumed to be linearly elastic, the solutions are usually more reliable only when the soil undergoes minor deformations, such as for very rigid clays and rocks at low shaking levels, with the exception of Bobet et al. [22]. The solutions for the transverse earthquake response are derived in the plane strain condition and therefore cannot be used for complex ground plans. In most cases, the contact interface is limited to two extreme states, full or no slip, while the lining is assumed to be continuous; therefore, a suitable representation of the segmental lining by an equivalent continuous lining is mandatory.

3. Application of random vibration theory in estimation of γ_max

Random vibration theory (RVT) relates the statistical properties of the random behavior of a dynamical system to the system properties or those of the random excitation. Therefore, RVT can be used to statistically estimate the random response of a system by representing the ground motion by a power spectral density (PSD) function.

Simplified theoretical conclusions are possible by assuming that ground motion is a stationary (i.e., the statistical properties of the motion are constant in time) Gaussian process. Although earthquake excitations are not stationary, the strong phase of such motions can be assumed to be stationary [35]. In this approach, the excitation is first defined by a PSD. The response PSD is either expressed theoretically or calculated using transfer functions. Then the statistical properties of the response are estimated using its PSD.

A well- known example of the use of RVT for the development of theoretical solutions is the Complete Quadratic Combination (CQC) method, which is useful for estimating peak displacements or forces within a structure [36]. CQC is also used for analyzing the nonstationary random responses of complex structures that are in an inhomogeneous stochastic field [37]. The analysis of the seismic response of linear multicolumn structural systems can be formulated by RVT, which takes into account the multicolumn input [38]. The steady-state filtered white noise model proposed by Kanai and Tajimi [39, 40] provides a well-known PSD in the field of earthquake engineering. White noise is a stationary random process that has a mean of zero and a constant spectral density for all frequencies. In the Kanai-Tajimi spectral model, the rock acceleration is assumed to be white noise and the overlying ground deposits are simulated by a linear one-degree-of-freedom system. Modified Kanai-Tajimi models are also proposed in the literature [41]. Therefore, RVT can be used to generate simple theoretical solutions. On the other hand, these simple solutions are limited to linear systems.

The theorems of random oscillation can be used to derive theoretical relationships between the parameters of dynamic response and ground motion. The theoretical analysis of the random response can be simplified by two assumptions. The first is that the excitation is statistically stationary in a broad sense. The second assumption is that the probability distribution of the excitation is Gaussian, so that each linear operation on this random process produces a different Gaussian process [42]. Although the properties of transient seismic motions obviously contradict these assumptions, the simplification can lead to reasonable theoretical functions that reflect the characteristic properties of dynamical systems. The applications concerning the combination of maximum modal displacements in structural dynamics [36, 43] and transfer functions for kinematic soil-structure interaction [44, 45] are well-known examples.

4. Conclusion

Analytical methods are implemented for analyzing underground structures by a numerous researchers. Though these methods have some shortcomings because of simplifying the design conditions, they provide a good approximation for preliminary analysis of such structures. Analytical methods are divided into two main categories: (a) soil-structure interaction and (b) free-field methods. In this chapter, free-field method, which ignores interaction between structure and encompassing soil, is being studied and its development has been discussed. For the practitioner, the simplified techniques are useful tools for preliminary studies. They make it simple to identify the variables that influence the severity of the prejudices, providing insight into the structure’s actions. Furthermore, the simplified approach and its solutions are invaluable in better understanding the relationship between dynamic loads, viscoelastic foundations, and tunnel structures, defining the most important parameters for the problem, and providing preliminary estimates or even a design. They also have the advantage of being able to conduct sensitivity analyses with little effort. The simplified approach may not be able to capture the responses and damage in structural specifics, components, or positions of possible failure due to the simplified assumptions for the tunnel layout and soil-tunnel interaction.

Conflict of interest

No potential conflict of interest was reported by the author.