Seismic vulnerability of circular tunnels in sand

The paper is focused on the assessment of seismic fragility curves for circular tunnels under moderate to severe earthquakes with the aim of improving the reliability of the risk assessment of underground infrastructural networks. Non-linear two-dimensional dynamic analyses were performed on different tunnel and soil configurations by using the finite-difference method implemented in Flac2D software. The applied input motions were selected considering their amplitude and frequency content variability. The response accelerations and predominant frequencies computed at ground level, above the tunnel, were compared with the corresponding free-field results to distinguish the effects attributable to the tunnel presence from those due to the site amplification. Tunnel safety was assessed through fragility curves, taking into account the dependence of tunnel lining bending resistance on the axial force variation during the earthquake. The more effective intensity measure was investigated correlating the tunnel performance to peak ground accelerations and peak ground velocities computed at the ground level and at the bedrock depth, respectively. The resulting fragility curves showed a satisfying matching with the empirical ones, generated on the basis of the observed seismic damage on tunnels.


INTRODUCTION
Nowadays, many countries in the world are facing a continuous demand for urbanisation, structuring the underground space to develop their physical interconnectivity (Broere, 2016). Lifeline networks, in particular, including transportation infrastructures and utility systems, have a strategic impact on the socio-economic conditions of a country, hence their seismic risk prediction, prevention and mitigation in the perspective of seismic urban planning and resilience is a crucial issue in seismically hazardous areas. Owing to their spatial extent and variability, linear infrastructure systems are vulnerable to multiple geotechnical hazards (i.e. seismic ground shaking, liquefaction, landslides and fault ruptures). Among other elements of such systems, tunnels of urban and extra-urban road, railway and highway transportation systems, or for other utility networks (i.e. waste water network), are particularly exposed to such hazards.
The object of this work is to develop numerical fragility curves for circular bored tunnels subjected to seismic ground shaking. This can induce damage in the tunnel structure, as observed after recent and historic seismic events (e.g. Dowding & Rozen, 1978;Sharma & Judd, 1991;Power et al., 1998;O'Rourke & Liu, 1999;Corigliano, 2007). The earthquake effect on the tunnel lining is also shown by experimental evidence from centrifuge models of tunnels undergoing shaking in the transverse section (e.g. Cilingir & Madabhushi, 2011a, 2011bLanzano et al., 2012), indicating that kinematic soil-structure interaction produces an increase of internal forces in the tunnel lining under dynamic shaking. Therefore, a probabilistic seismic vulnerability assessment study has been conducted for some ideal tunnel configurations, based on a numerical approach for the fragility analysis. In the perspective of the seismic risk analysis according to the principles of performance-based earthquake engineering (PBEE), the fragility function expresses the probability of exceeding different damage states (DS) for a given level of earthquake intensity measure (IM), accounting for the various uncertainties associated with the seismic hazard and the structural response estimation .
In this work, the cloud method (Jalayer et al., 2015) has been adopted to generate the proposed fragility curves expressed in terms of peak ground acceleration (PGA) and velocity. The method considered requires fewer seismic analyses to be performed and the use of unscaled records to avoid unrealistic and undesired modification of the scaled signals. Finally, the tunnel structure fragility was expressed in terms of the ratio between the loading and resistance bending moment, assuming its dependence on the axial force variation during earthquake.

NUMERICAL MODELS
Most of the available experimental measurements of seismic increments of the tunnel lining internal forces were collected in centrifuge tests on circular tunnels. Mainly for experimental reasons, dry sand was generally used in the tests to model the ground surrounding the tunnel (e.g. Cilingir & Madabhushi, 2011a, 2011bLanzano et al., 2012). Such tests refer to shallow circular tunnel models (i.e. cover-to-diameter ratio, C/D, lower than 3) embedded in either dense (D r = 75%) or loose sand (D r = 40%). Although limited, the range of cases covered by such experimental evidence is large enough to represent several cases of real shallow tunnels, where the pore pressure build-up during shaking may be neglected. Consequently, for the purposes of this study, the seismic soil-tunnel interaction was evaluated with reference to geometrical schemes and ground conditions similar to the experimental cases available in the literature. Therefore, a 10 m outer diameter and 0·50 m thick circular cross-section embedded in a 60 m thick dry sand layer, placed on a soft rock bedrock was considered. It is worth noting that, although the effect of the construction stages on the seismic behaviour of a tunnel lining in sand has been evidenced in published numerical studies (Fabozzi & Bilotta, 2016;Fabozzi, 2017;Bilotta, 2018), in this work the tunnel was assumed as wished in place and the excavation process was not simulated in the analyses. Two geometrical layouts were investigated, varying the tunnel axis depth from 15 m (C/D = 1) to 30 m (C/D = 2·5), as shown in Fig. 1(a). The analyses were performed on two idealised soil models representative of loose and dense sand, characterised by a relative density D r = 40% and D r = 75%, respectively. The assumed D r values correspond to the void ratios, e, reported in Table 1 and calculated on the basis of the typical values of e min and e max for clean sandthat is 0·5 and 1, respectively (Cubrinovski & Ishihara, 2002). The dependence of the shear wave velocity, V S , on the void ratio and depth was defined through the empirical power law correlation adopted by Hardin (1978), assuming the coefficients proposed by Silvestri & Lanzo (1999). The resulting V S profiles are shown in Fig. 1(b), highlighting the correction (Tropeano et al., 2018) close to the ground surface to avoid the unrealistic convergence to zero (dotted lines). In more detail, a linear trend was assumed, anchored to the tangent to the V S value at a depth of 10 m. The depth-dependent small-strain shear stiffness, G 0 , and bulk modulus, K 0 , were computed from the V S profile, setting the constant soil density, ρ, and Poisson's ratio, ν, as reported in Table 1. The shear resistance of the cohesionless loose and dense sand was assumed to be governed by the friction angle ϕ = 30°and ϕ = 35°, respectively. The non-linear and dissipative soil response was inferred from the variation of the normalised shear modulus, G/G 0 , and damping ratio, D, with shear strain, γ, resulting from laboratory tests performed on loose to dense sand by Seed & Idriss (1970). Fig. 1(c) shows, for instance, the mean and the envelope curves of the experimental data.
A linear behaviour was assumed for the underlying soft rock, with the typical physical and mechanical properties of a seismic bedrock.
The concrete strength class C25/30 was adopted to model the tunnel lining characterised by the Young's modulus, E, the Poisson's ratio, ν, the uniaxial compression, σ cy , and tensile strength, σ t , reported in Table 2.   Note that, even if the considered value of C/D = 2·5 does not correspond exactly to a deep tunnel, in the present work this definition is assumed only to distinguish the studied cases.
The two-dimensional models reproducing the four reference schemes were generated through the finite-difference code Flac (Itasca, 2011). The soil domain is 140 m wide and 60 m deep, based on the top of a stiffer 10 m thick layer simulating the bedrock. Solid elements were used to model the tunnel lining located in the middle of the numerical model, with the axis depth equal to 15 m or 30 m according to the geometric layout considered (see Fig. 1(a)). A no-slip condition modelled the surface contact between the lining and the soil, which seemed to be appropriate for a concrete lining cast in situ. This contact condition affects the axial loading of the tunnel lining during ground shaking, and consequently its bending capacity.
The mesh size, Δ, was calibrated according to the rule Δ , V S /8f max proposed by Kuhlemeyer & Lysmer (1973), to allow the typical seismic wave frequencies up to f max = 25 Hz to propagate reliably through the layered soil. Smaller mesh elements, instead, were used to model the tunnel lining. The infinite extension in depth of the bedrock was simulated by means of dashpots, providing viscous normal and shear stresses proportional to the volume and shear wave velocities of the bottom layer. To minimise the numerical model size, the so-called 'free-field' boundary conditions were applied to the lateral sides, simulating an ideal horizontal soil profile connected to the main-grid domain through viscous dashpots. As an example, Fig. 2 shows the numerical model of the shallow tunnel in dense sand reference scheme.
As numerous non-linear dynamic analyses need to be performed to generate fragility curves, the constitutive model directly included in the Flac library was chosen, to limit the computational effort. The bedrock was assumed elastic, while the elastic-perfectly plastic Mohr-Coulomb   constitutive model with shear stiffness increasing with depth was adopted to model the sand behaviour. According to this model, yielding occurs when the maximum shear strength is attained, the latter being a function of the mean stress acting in each element of the soil model. Conti et al. (2014) proved that, at least for the limited number of experimental cases in dry sand that were analysed, the Mohr-Coulomb model with non-linear and hysteretic behaviour provided cyclic changes of bending moment similar to those predicted by using a more refined constitutive model defined within the framework of bounding surface plasticity and critical state soil mechanics (Papadimitriou & Bouckovalas, 2002). In fact, although the shear-volumetric coupling affects the residual changes of internal forces at the end of shaking in dry sand (Bilotta et al., 2014a(Bilotta et al., , 2014b, this is less important for the prediction of the reversible changes of the internal forces generated in each cycle. Fig. 3 compares the maximum variation of bending moments provided by analyses in a cycle with both the above-mentioned constitutive models with centrifuge tests performed on dense and loose sand. Even though the experimental data are underestimated, the numerical results are in good agreement. Soil models were calibrated according to the physical and mechanical properties reported in Table 1. The unloadingreloading soil behaviour was defined by the hysteretic model provided by the Flac library and based on the well-known Masing rules (Masing, 1926). The decay of the normalised shear modulus with the increase of shear strain is introduced through the sigmoid function expressed in equation (1) Cyclic shear tests were numerically simulated on a cubic soil sample at different strain levels, from γ = 0·0001% to γ = 10%, to calibrate the parameters m, n and p in equation (1) on the G/G 0 -γ and D-γ curves obtained by Seed & Idriss (1970). The parameters obtained for the two soil models are shown in Fig. 1(c), revealing that a good approximation of the mean target stiffness (m = 1, n = À0·5 and p = À1·4) can lead to an excessive overestimation of the damping, due to the application of the Masing rules (Pelecanos et al., 2015). A more reasonable behaviour is predicted by the parameters m = 1, n = À0·6 and p = À1·2, selected for the dynamic analyses.
The effect of the Mohr-Coulomb criterion on the hysteretic soil behaviour was evaluated by simulating cyclic shear tests under different values of confining stress p′ = 113 kPa, p′ = 320 kPa, p′ = 533 kPa and p′ = 166 kPa, p′ = 333 kPa, p′ = 555 kPa, reproducing the lithostatic condition of loose and dense sand at a depth equal to 10 m, 30 m and 50 m, respectively. The sensitivity analyses revealed that the introduction of the Mohr-Coulomb criterion implies the magnification of the shear stiffness degradation and the damping increase in correspondence of shear strain values higher than 0·3% for the shallowest soil samples, while no differences are recognised at lower strain levels. Since shear strains higher than 0·1% are expected to occur occasionally, for example under the peaks of input motions (see for instance Fig. 6, later), the cyclic behaviour is expected to be slightly influenced by the constitutive assumptions.
The low strain damping was introduced through the Rayleigh (1945) approach, assigning to the sand and bedrock layers the frequency-dependent damping ratio with the minimum D 0 = 0·5% at the soil predominant frequency. The adopted D 0 value is consistent with the low strain damping provided by Seed & Idriss (1970).
For the tunnel lining, typical values for a Class 25/30 concrete mix were adopted (BSI, 2004) and the damping D 0 was introduced through the Rayleigh (1945) approach, as had already been done for the soil low-strain damping.
It is worth noting how Kampas et al. (2019) previously reviewed a number of typical tunnel lining modelling approaches and the associated effect on axial forces and acceleration fields in the model of soil and lining pre-and post-cracking. A comparison between the results revealed that the lining axial forces and acceleration would not be significantly affected by the modelling approach. In this study, the structural material of the lining was assumed visco-elastic with the density, Young's modulus, Poisson's coefficient and damping ratio D 0 reported in Table 2.
Accelerations and shear strain-time histories along the vertical sections V1, V2 and V3, and horizontal and vertical stresses in the mesh elements of the tunnel sections in Fig. 2 were recorded. The former allowed the dynamic behaviour of the soil-tunnel system to be identified and the latter were used to compute the variation of the internal forces induced by the seismic motion.

FRAGILITY CURVES Method
The computation of fragility curves requires (a) the definition of the engineering demand parameter, EDP, that measures the performance of the analysed system, and (b) the choice of the seismic hazard IM. The EDP is defined by a threshold expressing the accepted damage levelthat is DE SILVA, FABOZZI, NIKITAS, BILOTTA AND FUENTES the limit stateto compute the demand to capacity ratio, DCR LS . The definition of such parameters and of their correlation depends on the failure mechanism expected to occur during the earthquake. In this study, the cloud method was chosen to generate the fragility curves. Compared with the more commonly adopted stripe analyses, this method requires a lower number of seismic analyses to be carried out and it allows the adoption of unscaled records, thus avoiding the unrealistic and undesired modifications of scaled signals (Jalayer et al., 2015).
Moreover, in the studies of the tunnel seismic fragility available in the literature (e.g. Argyroudis & Pitilakis, 2012;Argyroudis et al., 2017;Fabozzi et al., 2017), the probability of failure is computed as the probability of exceeding an IM value associated with the occurrence of the DS and a two-parameters fragility model is appliedthat is the logarithm of the mean IM value and the standard deviation. The upgrade provided by the cloud method draws on the fact that the probability of failure is computed as the probability of exceeding the DCR LS threshold on the basis of a lognormal distribution of DCR LS , whose mean value is linearly dependent on IM in the logarithmic scale. Based on such hypotheses, the mean system performance at each IM and its standard deviation can be computed as follows in equations (2) and (3), respectively The vector of the parameters χ = [log(a), b, σ log DCRLS|IM ] is derived from the best fitting of the data D = {(IM i , DCR LSi ), i = 1:N} resulting from the performed N non-linear analyses. Equation (3) highlights that in the cloud method the uncertainties in the seismic performance are evaluated as the standard error in the estimation of the regression linking IM to DCR LS .
The probability of failure is then computed from the standard normal distribution shown in equation (4) The regression predictions are reliable if the failure mechanism is mobilised in most of the dynamic analyses performed. This condition is readily satisfied when the input motions are scaled to a different hazard level, as in the multiple stripe method. However, as already stated above, the use of unscaled natural records was preferred in this work, provided that most of the selected records were of sufficient magnitude to mobilise the investigated failure mechanism.

Selected input motions
Ten, moderate to severe, natural accelerograms were selected from the Pacific Earthquake Engineering Research Center (PEER) ground motion database to perform nonlinear dynamic analyses. Fig. 4 shows the ten natural accelerograms ( Fig. 4(a)) and the associated acceleration response spectra (Fig. 4(b)).
The main characteristics of the selected signals are reported in Table 3: earthquake and station name, date, component, moment magnitude, M W , source-station epicentral distance, R epi , equivalent shear wave velocity in the shallowest 30 m below the station, V S30 , PGA, peak ground velocity (PGV) and predominant period, T P .
Input motions should be selected to cover a large range of IM values, in order to reduce the standard error in the estimation of the regression between IM and DCR LS . However, Jalayer et al. (2015Jalayer et al. ( , 2017 revealed that the improvement of the DCR LS -IM regression prediction tends to reduce for more than ten simulations. Consequently, additional simulations were not performed, owing to the huge effort corresponding to the execution of the dynamic analyses. Consequently, the selection was based on the variability of PGA, PGV and spectral accelerations associated with the first resonance period of the soil (T 1 in Fig. 4(b)), as well as of the frequency content, as proven by the variable predominant period, T P , in Table 3. Furthermore, the V S30 in Table 3 proves that the selected signals were recorded on sufficiently stiff soil, identified as A or B type in the Eurocode, thus avoiding the possible signal modifications induced by local site effects. A Butterworth fourth-order filter was applied to all signals to limit their frequency content to the range 0·1-25 Hz.

Soil response and site effects under non-linear dynamic analyses
The initial stress conditions were reproduced into two steps.
(a) Computation of the lithostatic stress acting in free-field conditions.   (Fig. 5(a)) and deep ( Fig. 5(b)) tunnels in loose and dense sand. Independently of the reference scheme, the amplification functions of the full models are close together along the three vertical sections, with the first resonant frequency equal to 1 Hz and 1·7 Hz for loose and dense sand, respectively. The numerical models are validated by the computed frequency values, which are very close to the theoretic predominant frequencies of a layer in which the stiffness increases with depth (Dobry et al., A* indicates that soil station is identified as rigid, but no field investigations were performed to verify its stiffness.   , resulting in f = 1·08 Hz and f = 1·68 Hz for loose and dense soil, respectively. The peaks associated with the predominant frequencies are not affected by the presence of the tunnel, nor is the whole dynamic response along the vertical sections V2 and V3. On the other hand, a drop at high frequencies is recognised along the vertical V1 with respect to the amplification functions computed before the excavation (dashed lines in Fig. 5). The shadow zone of loose and dense sand is, respectively, limited by a frequency f min % 9 Hz and f min % 14 Hz, independently of the tunnel depth. Since the attenuation of the motion occurs at moderate to high frequency, near-fault ground motions are expected to be reduced by the presence of the tunnel. On the contrary, smaller effects are expected at the far field due to the loss of the high-frequency content of the propagated signal. Figure 6(a) summarises the results of the non-linear dynamic analyses in terms of amplification factorsthat is the ratio between the PGA computed at the ground level (PGA gl ) and on the top of the bedrock layer (PGA im ). Independently of the reference scheme, the amplification factor reduces when increasing the amplitude of the input motion, due to the soil hysteretic behaviour. As expected, the effect of non-linearity is so pronounced that an attenuation of the seismic motion is observed in the loose sand under the strongest five earthquakes. Even though the same decay of shear stiffness (see Fig. 1(c)) was applied to the soil profiles, no attenuation of the seismic motion is observed in the dense sand, due to its higher initial stiffness. As an example, the contours of the maximum shear strain, γ max , developed around V1 under the earthquake number 1 are shown in Figs 6(b)-6(e), demonstrating that the shear strains are much larger and more diffused in the loose sand models.
Linear regressions were fitted through the computed amplification factors along the three control vertical sections as shown in Fig. 6, for an easy comparison of the different behaviours. In most of the cases for a deep tunnel, the highest amplifications occur along the vertical V2 (light grey continuous line). When the tunnel is deep and the sand is dense, the presence of a void zone (i.e. V1dashed line) reduces the ground motion with respect to the free-field results (V3continuous black line). This trend most probably confirms the 'shadowing effect' observed under the low-amplitude excitation (see Fig. 5 in the earlier section entitled 'Numerical models'), due to the loss of highfrequency content when the seismic wave crosses the tunnel.
A different trend is recognised for the shallow tunnel in dense sand for all PGAs and in loose sands for PGA values greater than 0·3g approximately. In these cases, the maximum amplifications correspond to V1 due to soil non-linear behaviour. As evidenced by the contours shown in Figs 6(b) and 6(c), the higher strain level mobilised around the lining reaches the ground level only when the tunnel is shallow. As a consequence, the ground motion is differently modified depending on the co-existence between the 'shadowing effect' of void and the soil non-linearity.

Earthquake-induced forces in the lining
During the seismic analyses, the stresses acting in correspondence to the tunnel lining control sections   Fig. 2 were recorded to compute the evolution of axial force and the bending moment. Fig. 7 compares the static load (black lines) with the maximum and minimum profiles (grey lines) induced by each earthquake in the lining of the four reference schemes. In all cases, the maximum variation of the internal forces is recognised in sections 45°, 135°, 225°and 315°, with significant positive and negative peaks in the bending moment. This trend is indicative of the expected ovalisation of the cross-section under the exciting cyclic shear strain. A larger increase of internal forces affects deep tunnels, in turn affecting the safety. Significant reductions of the axial forces are shown in shallow tunnels, critically approaching zero. Since the bending resistance typically increases with axial force, its reduction influences the tunnel safety.
These numerical results have been compared with the experimental results of the centrifuge tests carried out by Lanzano et al. (2012) and, when possible, with their numerical interpretations available in the literature. Fig. 8 for instance, shows the comparison in terms of maximum cyclic change of bending moment, ΔM max , per unit length of tunnel lining, reached during the seismic excitation at the most loaded tunnel sections, as a function of the peak value of the input time history of acceleration applied at the base of the model, PGA im . Fig. 8(a) refers to the case of shallow tunnel (C/D = 1), for dense and loose sand, corresponding to T1 and T2 models Lanzano et al. (2012), while Fig. 8(b) concerns the case of deep tunnel (C/D = 2) for both dense and loose sand cases, corresponding to T3 and T4 centrifuge models in Lanzano et al. (2012), respectively. In the figure, ΔM max has been made non-dimensional by dividing it by τ eq r 2 À Á , where r is the radius of the tunnel and τ eq is the equivalent shear stress. The latter is representative of the shear stress transmitted to the tunnel lining during shaking and it was calculated from the horizontal equilibrium of a soil column undergoing inertial forces through equation (5): a cc,max is the maximum acceleration at the mid-depth of the tunnel axis; g is acceleration due to gravity; γ soil is the unit weight of soil; and H axis is the depth of the tunnel axis.
The numerical results shown in Fig. 8 Wang (1993) (which extend Hoeg (1968) are also shown, assuming equation (6) ΔM max τ eq r 2 ¼ where to consider the non-linear soil stiffness decay during shaking, the flexibility ratio, F, is calculated as in equation (7) F ¼ G  where G soil;mob is the mobilised shear modulus that changes case by case and, for the sake of simplicity, is assumed either equal to its maximum value, G 0 (continuous line), which provides a low estimate, or equal to a conventional value of 0·3G 0 , which provides a realistic estimation of the highest possible value in all the cases (dashed line). Despite the influence of soil plasticity on the accumulation of residual bending moments at the end of shaking (experimentally confirmed by Cilingir & Madabhushi (2011a) and Lanzano et al. (2012)), the comparison shown in Fig. 8 of numerical predictions (that include different plasticity models) with the classical Wang's solutions (that assume an equivalent linear behaviour) is legitimate since the calculated steady-state cyclic changes of bending moments occurring during shaking are considered.

Centrifuge tests Present study
The analytical solution largely overestimates the centrifuge tests results (black markers) and the associated numerical predictions (grey markers). This is expected, according to Kontoe et al. (2014), considering that the tunnel lining tested in the centrifuge was extremely flexible (F ranges between 644 (dense sand) and 322 (loose sand) for the initial shear modulus G 0 , and between 193 (dense sand) and 96 (loose sand) for the mobilised shear modulus 0·3G 0 ).
The numerical results of this study (empty markers), that take into account stiffer linings (F ranges between 37 (dense sand) and 11 (loose sand) for G 0 , and between 11 (dense) and 3 (loose sand) for 0·3G 0 ), are closer to the analytical solutions, and in many cases fall within their observed range.

Discussion of results
For strategic linear infrastructures, minor damage should be prevented and, hence, the fragility curves of tunnels were computed with reference to the first onset of cracks in the lining. Such a limit state is associated with reaching the concrete tensile strength, σ t (see Table 2). The resistance was consequently calculated as the moment (M res (t) in Fig. 9(a)) induced in the resistant section by the stress diagram plotted in Fig. 9(a), in which σ c is the maximum compression stress acting in the section, while σ ru and σ rl is the stress acting in the upper and lower steel reinforcements, respectively. Since the response of the structural materials is still elastic at the onset of the assumed limit state, σ ru and σ rl are typically normalised through the homogenisation coefficient n c = 6·7, calculated from the ratio between the steel and concrete Young's moduli. All the stress values can be easily calculated from geometric considerations, once the depth of the neutral axis (x(t) in Fig. 9(a)) is known from assuming equilibrium with the axial force N(t). As observed in the section 'Soil response and site effects under non-linear dynamic analyses', the axial force changes during the seismic excitation, hence the resistant moment M res (t) is also time dependent. This is a critical aspect of the seismic tunnel performance considered in this study, where the static equilibrium was solved at each time step using a routine implemented in Matlab. As an example, the evolution of M res (t) under earthquake number 1 for the shallow tunnel in loose soil is compared with the loading moment M load (t) in Fig. 9(b); for the sake of clarity, the negative resistant moment is omitted, as it is equal and opposite to the positive one due to the symmetry of the resistant section. As expected, the numerous and instantaneous decreases of N(t) reduce the resistant moments and, in most cases, the lowest values occur when the loading moment shows the most significant peaks. The seismic performance of the tunnel was 'measured' by the maximum ratio between the loading and the resistance momentthat is the maximum mobilised strength during earthquake Since the bending loads are responsible for the tunnel damage, demand and capacity are generally identified with the loading and resisting bending moment. In the existing studies on the fragility of tunnels (e.g. Argyroudis & Pitilakis, 2012;Argyroudis et al., 2017), only the peak loading values are assumed in the calculation of the DCR LS and the evolution of the capacity associated with the variation of the axial load during the whole excitation is neglected. Table 4 reports the mean values and the standard deviations, σ, computed in the four tunnel sections (45°, 135°, 225°and 315°), in which the highest increment of internal forces was recognised during the dynamic analyses (see Fig. 7). In all cases, 30% of strength is mobilised on average, with values exceeding 40% in the case of a shallow tunnel in dense sand.
Interpreting the bending moment as an EDP, the limit state was assumed to be reached when DCR LS = 1. Thus, the probability of failure was calculated through the cloud method (see the section entitled 'Method'), assuming the PGA and PGV of the input motion (PGA im and PGV im ) and recorded at the surface in free-field conditions (PGA gl and PGV gl ), as seismic IMs.
Even though the numerical fragility curves available in literature assume the PGA of the ground surface as IM, the  correlation of the demand-to-capacity ratio was checked, in addition (a) against the PGV at the surface, because the ovalisation of the transverse section of the tunnel lining is governed by the shear strains close to the tunnel, which is correlated to the soil ground velocity (b) against the PGA and PGV of the input motion, because fragility curves expressed as a function of the input motion parameter are more useful for real-time analyses, since signals provided by the seismic stations are recorded on stiff rock outcrop and consequently do not include site effects.
The envelopes of the fragility curves computed for the control sections of each reference scheme are shown in Fig. 10 as a function of the PGA (Fig. 10(a)) and PGV (Fig. 10(c)) of the input motion (i.e. of a stiff and flat rock outcrop) and of the PGA (Fig. 10(b)) and PGV ( Fig. 10(d)) at the ground level (i.e. accounting for the site effects). Table 5 reports the determination coefficient of the DCR LS -IM regression in the logarithmic scale together with the three parameters of the fragility model. As testified by its highest R 2 values, PGV gl captures the relevant features affecting the seismic performance of soiltunnel systems.
When the parameters of the input motion PGA im and PGV im are assumed as IMs, there is a scatter between the probability of failure of tunnels in dense and loose sand associated with the different site amplifications. In fact, high amplitudes of input motion are amplified by dense sand, but they are dampened by loose soil (see Fig. 6), thus the seismic loads acting in the lining are reduced and, consequently, result in a lower fragility of the tunnel embedded in loose sand. As a matter of fact, the scatter is evidently reduced, when the curves are expressed in terms of ground motion parameters PGA gl and PGV gl . If the occurring earthquake is strong enough to produce the same high amplitudes at the ground level of the different soil profiles, damage is expected to tunnels in dense as well as in loose soil.
The computation of fragility curves in terms of PGA im was repeated assuming zero axial force acting in the lining and consequently the lowest bending resistance of the section, which is superimposed on Fig. 8 through the horizontal line. The resulting curves are plotted using thinner lines in Fig. 10, confirming the beneficial effect of the axial forces in the reduction of the fragility.
The reliability of the numerical curves was checked through comparison with empirical fragility curves generated from the damage observed on tunnels after earthquakes. The empirical curve produced by ALA (2001) in terms of PGA gl and relevant to minor damagethat is minor cracking of tunnel lining, minor rock falls, spalling of shotcrete or other supporting materialsis superimposed on the plot in Fig. 10(b). The comparison is satisfactory.
Moreover, superimposed on Fig. 10(d), there are the curves computed by Corigliano (2007) as a function of the PGV gl for deep tunnels corresponding to no or slight damage (i.e. cracking of the concrete lining) and moderate damage (e.g. spalling of the concrete lining, liner steel exposed, cracking and crushing of the concrete lining, etc.). Consistently with the class of the deep tunnel for which it was obtained, the empirical curve associated with moderate damage is practically coincident with the numerical curve corresponding to deep tunnels in dense soil and very close to the ones corresponding to deep tunnels in loose soil. On the other hand, the empirical curve for slight damage is far from the numerical predictions, probably because the first is associated with no or negligible effects, whereas the others are associated with the damage being sustained. The numerical fragility curves showed a satisfactory agreement with the empirical curves, although the latter were obtained (a) by applying a different fragility model, for which the probability of failure is computed as the probability of exceeding the IM when the DS was observed (b) by collecting damage data felt by different types of tunnel lining embedded in heterogeneous soil types.

CONCLUSIONS
This paper investigates the dynamic soil-tunnel interaction mechanisms under different earthquakes with amplitudes expressed in terms of PGA, ranging from 0·05g to 0·515g. Soil stiffness and tunnel depth were also varied to investigate their effects.
The presence of the tunnel was found to attenuate the surface motion since the high-frequency content of seismic waves is filtered by the void. This effect is only relevant under weak motions and tends to decrease when strong motions mobilise the soil non-linearity.
Earthquake-induced variations of internal forces were recognised in the lining, with greater changes in the deep tunnel, where stress levels are higher. Since critical reductions of the axial forces occurred in the shallow tunnels, the safety of the tunnel was assessed during the whole duration of the seismic excitation, taking into account the dependence of the bending resistance on the axial force.
The results were expressed in terms of fragility curves, considering different seismic IMs. The best correlation with the seismic performance of tunnels was obtained by the PGV recorded at the surface in free-field conditions, owing to its direct correlation with shear strains. Even though its effectiveness is shown, such a parameter is affected by the local site response, which is unknown immediately after an earthquake. For this reason, the same curves were expressed as a function of the input motion parameters to facilitate their application in real-time analyses in the post-earthquake management of infrastructures (Fabozzi et al., 2018;Fabozzi et al., 2019). The comparison revealed that tunnels in loose sand are less vulnerable to severe input motions, thanks to the seismic attenuation of soil non-linearity.
Furthermore, the satisfactory comparison with the empirical fragility curves confirms the hypothesis that the variation of axial force plays an important role in the seismic performance of tunnels.
Overall, the suitable evaluation of the seismic vulnerability of strategic infrastructure like tunnels has undoubtedly an important impact on the socio-economic resilience of a specific system or community. It should therefore be an integral part of the risk identification and management of linear infrastructure. The current availability of fragility curves for tunnels is very limited. This has a negative impact on the inclusion of these structures in streamlined risk management tools using geographical information systems or loss-driven earthquake early warning and rapid response systems (Fabozzi et al., 2018).
It seems worth noticing that one of the most recent developments of performance-based earthquake engineering is the expression of the annual rate of failure, obtained from the combined calculation of the fragility of the system under investigation and the seismic hazard on site. In its most refined form, the approach is intended to include in the definition of seismic risk site-specific effects that are neglected by the traditional uncoupled formulation of hazard and fragility. This approach is becoming increasingly common for assessing the seismic risk on structures (Jalayer & Cornell, 2009;Iervolino et al., 2017), while a few pioneering applications have been developed by geotechnical researchers (Kramer, 2008(Kramer, , 2014. In this vein, fragility curves like those presented in the present work, generated from full dynamic analyses, improve the reliability of risk assessment, as they directly take into account site amplifications and soilstructure interaction.