Amplification effect of near-field ground motion around deep tunnels based on finite fracturing seismic source model

Qiankuan Wang,Shili Qiu,Yao Chng,Shaojun Li,Ping Li,Yong Huang,Shirui Zhang

a State Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan,430071,China

b Faculty of Engineering,China University of Geosciences,Wuhan,430074,China

c University of Chinese Academy of Sciences,Beijing,100049,China

d South-Central University for Nationalities,Wuhan,430074,China

e China Railway First Survey and Design Institute Group Ltd.,Xi’an,710000,China

Keywords:Near-field ground motion Amplification effect Seismic waves Deep tunnel Rockburst

ABSTRACT Dynamic failure of rock masses around deep tunnels,such as fault-slip rockburst and seismic-induced collapse,can pose a significant threat to tunnel construction safety.One of the most significant factors that control the accuracy of its risk assessment is the estimation of the ground motion around a tunnel caused by seismicity events.In general,the characteristic parameters of ground motion are estimated in terms of empirical scaling laws.However,these scaling laws make it difficult to accurately estimate the near-field ground motion parameters because the roles of control factors,such as tunnel geometry,damage zone distribution,and seismic source parameters,are not considered.For this,the finite fracturing seismic source model (FFSSM) proposed in this study is used to simulate the near-field ground motion characteristics around deep tunnels.Then,the amplification effects of ground motion caused by the interaction between seismic waves and deep tunnels and corresponding control factors are studied.The control effects of four factors on the near-field ground motion amplification effect are analyzed,including the main seismic source wavelength,tunnel span,tunnel shape,and range of damage zones.An empirical formula for the maximum amplification factor (αm) of the near-field ground motion around deep tunnels is proposed,which consists of four control factors,i.e.the wavelength control factor (Fλ),tunnel span factor (FD),tunnel shape factor (Fs) and excavation damage factor (Fd).This empirical formula provides an easy approach for accurately estimating the ground motion parameters in seismicityprone regimes and the rock support design of deep tunnels under dynamic loads.

1.Introduction

With rapid developments of underground works,seismicity hazards induced by high in situ stresses have become a challenging issue,which can pose a substantial threat to the safety of engineering construction.On one hand,excavation unloading can induce stress redistribution and associated areas characterized by high in situ stresses.On the other hand,the crack initiation,development,and expansion in the surrounding rocks of the tunnel will be induced,resulting in a decrease in rock strength,which can form an excavation damaged zone (EDZ) and excavation fracture zone (EFZ).This would provide favorable conditions for the occurrence of seismic hazards (Kaiser et al.,1996; Vazaios et al.,2019; Sun et al.,2021).Under this condition,the overstressed rock masses around tunnels are prone to fracturing,resulting in seismic events.As seismic waves propagate to the boundary of the tunnel,the ground motion increases due to site amplification effects,and the additional vibration energy can then be fed into rock masses.This will directly shake loose rocks or evoke rock ejection in these high burst-prone areas,thereby triggering three types of rockbursts (Wang and Cai,2014,2015).The three types of rockbursts are bulking caused by surrounding rock fracturing,rockfall induced by ground motion,and ejection caused by seismic energy transfer(Kaiser et al.,1996).The rockburst mechanisms induced by the seismic event around the deep tunnel are significantly complicated.Two factors are dominant accounting for the strong ground motion parameters,e.g.the peak particle velocity (PPV) of seismic events and the amplification effect of the excavation structure.

The parameter PPV is widely used in strong ground motion prediction,commonly as an indicator of rockburst induced by a seismic event to identify the damaged zones of rock masses or rock supports and to delineate their corresponding damage influence radius (Feng et al.,2012; Morissette and Hadjigeorgiou,2019).Generally,the PPV is estimated by simple empirical scaling laws proposed by McGarr et al.(1981),McGarr (1984),Spottiswoode(1984),Hedley (1988,1992),Kaiser et al.(1996),Durrheim et al.(2005) and Potvin and Wesseloo (2013),as listed in Table 1.Most of the PPV empirical scaling laws(Nos.1-5 in Table 1)are proposed using far-field ground motion monitoring data without considering the near-field saturation and amplification effect.As a result,they may not be suitable for estimation of the near-field ground motion effects around deep tunnels(Kaiser et al.,1996;Wang and Cai,2014,2015).Through combining the far-field PPV scaling law with the near-field saturation effect,Durrheim et al.(2005) and Potvin and Wesseloo (2013) proposed two empirical formulas applicable to both far-and near-field PPV calculations(e.g.Nos.6-7 in Table 1).However,these two relations do not consider the dynamic fracturing process of seismic events,uneven radiation of seismic waves,and dynamic effects between seismic waves and underground excavation structures.

In other words,by using empirical scaling laws in Table 1,the estimated near-field PPVs and their distribution around deep tunnels and caverns usually may not be precise.Based on the study of dynamically triggered rockbursts induced by seismic events in deep tunnels,due to the existence of site effect,there are some significant differences between the measured ground motions in the field and the results calculated from the empirical scaling laws.Milev et al.(1999) and Durrheim (2012) analyzed the rockburst damage on the tunnel surface in South African mines,and they found that the field PPV is 4-10 times that calculated by the empirical scaling law.Ortlepp (1993) and Stacey and Rojas (2013)also found that the initial ejection velocity of rock blocks on tunnel boundary exceeds 10 m/s and is much higher than that estimated by the empirical scaling law indicating the existence of amplification effect.For this,Kaiser et al.(1996) pointed out that the increase of ejection velocity of rock blocks is derived from the release and transformation of stored strain energy around the tunnel and that the amplification factor(α,which refers to the PPV ratio as between the excavation model and the unexcavated background model at the same location) should be 1-4.However,Durrheim et al.(1996) attributed the amplification effect to resonance,and Linkov (2000) believed that the amplification effect is caused by energy release due to rock softening.

Table 1 A summary of the scaling laws to assess the PPV.

For Long-Victor Mine,Mikula (2013) stated that the average amplification factor can be set to 3.Zhang et al.(2015)used a onedimensional (1D) model with a fracture zone to describe the fractured rock mass.The amplification factor in terms of PPV shows that the PPV of the model’s free end is 2-3.6 times that of the model’s input end due to the magnification of the fracture zone.Based on the rockburst cases in South African mines,Wesseloo(2018) proposed PPVR=αPPVBfor assessing the real peak velocity(PPVR)of the tunnel surface where PPVBis the theoretical value of the peak velocity without considering the interaction between seismic waves and tunnels.In the absence of engineering parameters,considering the wave reflection on the tunnel surface,α is the theoretical amplification factor usually with a value of 2.

There are significant differences between the actual and the above-mentioned amplification factors for practical cases.This is due to multiple factors,including the radiation patterns,complex interactions of seismic waves with geological structures and excavation structures,influence of surface waves,frequency of incident seismic waves,and the shape and span of excavation structure(Kaiser et al.,1996;Tang and Xia,2010;Diao et al.,2011;Chen et al.,2013; Potvin and Wesseloo,2013; Wang and Cai,2015; Almadani et al.,2020; Wang et al.,2021).Based on the literature review of the site effect of ground motion and the authors’ investigations,itseems that the amplification effect is affected mainly by structural factors and seismic source factors.The influence of seismic sources is primarily manifested by the wavelength (St John and Zahrah,1987; Hashash et al.,2001; Wang and Cai,2014,2015),while structural factors mainly include three factors:tunnel span,tunnel shape and distribution of damage zones (Kaiser et al.,1996; Milev and Spottiswood,2005; Wang and Cai,2014,2015; Zhang et al.,2015).For this,the study focuses on the control effect of four factors (i.e.wavelength,tunnel span,tunnel shape,and range of damage zones) on the amplification effect and distribution characteristics of ground motion.

The strong ground motion in rocks around deep tunnels occurs usually near a fracture source,or even within the fracture source region.Hence,a near-field finite fracturing seismic source should be introduced to analyze the strong ground motion characteristics around the tunnel.However,it is often assumed in the empirical scaling laws(see Table 1)that the seismic sources are far-field point sources with spherical radiation patterns.Only the geometry and material attenuation of seismic waves are analyzed(i.e.the PPV is a function of the seismic source distance).In other words,the site effects such as the amplification and shielding of the excavation structure,which are created by reflection,refraction,and diffraction of seismic waves around the excavation structure,are ignored(Milev et al.,2001;Durrheim,2012).The assumptions are incorrect in the near-field ground motion analysis.The distance between the concentrated seismic events and the tunnel structure is usually 2-3 times the seismic source size.The seismic source size of largemagnitude seismic events is generally relatively large,even comparable to the above-mentioned distance (Hudyma et al.,2003;Mutke,2008; Zhang et al.,2012; Ma et al.,2015).It means that when analyzing near-field ground motions induced by near-tunnel seismic events,the controlling role of source sizes should be considered,rather than simply using a point source model to estimate ground motion parameters.For this,the source size effect should be considered and a finite-size fracture source should be used to characterize the focal mechanism.Furthermore,the dynamic fracturing process of a finite-size seismic source,the propagation law of seismic waves excited by a finite-size seismic source in the near-field area,and the amplification effect of ground motion in the excavation structure area should be considered.However,in previous studies,the finite fracturing seismic source was rarely used in near-field ground motion analysis in the underground engineering field because associated theory and methods were insufficient.

Given this,this study introduced a recently proposed finite fracturing seismic source model (FFSSM) to characterize the radiation seismic wave field evoked by dynamic fracturing processes of finite dimension sources around deep tunnels.Based on the radiation seismic wave field,the amplification effect of seismic waves,caused by its interaction with the underground excavation structure,was simulated by the spectral element method (SEM) software Specfem2D.The near tunnel’s strong ground motion distribution law was also analyzed and revealed.The emphasis focused on the influence of four control factors on the amplification effect of near-field ground motion.Finally,an empirical formula was established to estimate the amplification factor of near-field ground motion around deep tunnels.

2.A brief introduction to FFSSM and Specfem2D

2.1.The FFSSM and parameters calculation procedures

When analyzing the ground motion of far-field seismic events,the effect of source size on the seismic response between seismic waves and tunnels is almost negligible,so a point source model can be used to characterize the seismic source.However,when simulating wave propagation induced by near-field seismic events,the seismic source must be expressed as a finite dimension seismic source.The size of the seismic source,dynamic fracturing processes,and wave interactions must be considered.For this,the FFSSM was developed,which references the point seismic source model and the finite fault model(Beresnev and Atkinson,1997;Aki and Richards,2002; Li et al.,2014).Considering the dynamic fracturing process of the seismic source,this model can effectively characterize the seismic wave field characteristics of the near-field seismic sources.

According to seismic wave propagation,the observable seismic energy of the pressure-shear vertical wave system(P-SV system)is much greater than that of the pressure-shear horizontal wave system (P-SH system) (Hashash et al.,2001).Thus,a twodimensional (2D) plane approximation of the P-SV system can be used to analyze the seismic wave propagation in deep tunnels.Fig.1a is the basic principle of a 2D FFSSM:the dynamic fracturing process of the main seismic source is characterized by dividing a fault or fracture(main seismic source,the length is L0)into a finite number of seismic sub-sources(see SS#1-SS#5 in Fig.1a)with the same length (Le),different dynamic corner frequencies (f0) and time shifts (ti).tirepresents the moment of excitation of the ith seismic sub-source.The number of seismic sub-sources (n) should be determined according to the size of the main seismic source.

Fig.1.A conceptual diagram of the FFSSM:(a)A conceptual diagram for the FFSSM-2D with the characteristics of central epicenter excitation and bidirectional rupture propagation in P-SV system;and(b)A conceptual diagram for the FFSSM-3D in the X-Y-Z coordinate system.φs,γr and δd represents the strike,rake,and dip of the fault,respectively.R1 and R2 are two receivers at a certain distance away from the FFSSM.Ux,Uy,Uz represent the pressure(P),shear horizontal(SH),and shear vertical(SV)component of the displacement at the position of the receivers.

The FFSSM has two special modes: bilateral fracturing and delayed fracturing modes.The bilateral fracturing mode means that the fracture propagates from the initial fracturing seismic subsource (SS #1 in Fig.1a) to the follow-up activated seismic subsources (SS #2-SS #5 in Fig.1a) on both sides simultaneously.In other words,the fracturing sequence is SS#1→SS#2-#3→SS#4-SS#5.The red and orange arrows in Fig.1a show the main direction of fracturing propagation of the seismic sub-source.The delayed fracturing mode means that the seismic sub-sources fracture sequentially,with different time shifts according to the fracturing sequence.The initially activated sub-source (SS #1,with crack initiation time t1=0 s)fractures firstly and starts to radiate seismic waves.As the fracture spreads towards other seismic sub-sources,these follow-up seismic sub-sources are activated one after another(For example,SS#2-#5 are excited at t2-t5).After all the seismic sub-sources are activated,the dynamic fracturing process is completed.In particular,it is noted that the FFSSM utilizes a moment tensor point source model to quantify the seismic source characteristics and seismic wave radiation patterns of a seismic sub-source.

Although only a 2D FFSSM(FFSSM-2D)is used in this paper,the FFSSM can naturally be extended to a three-dimensional (3D)FFSSM(FFSSM-3D),as shown in Fig.1b.The fracture surface of the main source is divided into many rectangular source sub-fracture surfaces of the same size.The red seismic sub-source fracture surface,where the dark red hexagram is located,is the fracturing initiation site,and the fracturing start time is t1.The orange and yellow arrows indicate the fracturing direction of seismic source fracture;that is,the fracturing sequence of the rectangular seismic sub-sources is from red rectangular sub-source (SS #1) to orange rectangular sub-sources(SS#2-#9)and then to yellow rectangular sub-sources (SS #10-SS #21),as shown in Fig.1b.The activating and fracturing of these rectangular sub-sources are used to simulate the dynamic fracturing process of the fracture surface of the main seismic source.Compared with the FFSSM-2D,the FFSSM-3D has many advantages: (1) The 3D model can study more complex geometric conditions of the seismic source,i.e.a spatial relationship between the seismic source and the tunnel;(2)More complex geological conditions such as faults,joints,soft and hard rock contact surfaces can be characterized so that the control factors of site effect of ground motion can be analyzed; (3) The 3D model research object is no longer just an excavation section,but a tunnel or cavern group with a certain length.It is more convenient to assess the ground motion risk along the tunnel.However,the FFSSM-3D also has its limitations.Most importantly,the increase in the number of computational elements and the complexity of geological conditions and spatial relationships of the 3D model will inevitably lead to a substantial increase in modeling period and calculation time.Moreover,for multi-factor study,the FFSSM-3D model modification also is time-consuming and laborious.This paper adopts a FFSSM-2D model considering time cost and research requirements.

The parameters to be determined in the FFSSM-2D include the seismic sub-source length (Le),the number of seismic sub-sources(n),the sub-seismic moment (Me),the dynamic corner frequency(f0),and the time shift (ti).There are two different calculation procedures for the parameters of the FFSSM that are selected according to the different research contents and interpretation methods(Fig.2).This paper focuses on the forward simulation and parameter research,setting different seismic cases to study the dynamic wave interaction between seismic waves and underground excavation structures(including the tunnel and the damage zones formed during the excavation of deep tunnels).Thus,Path#1 in Fig.2 is followed.

2.2.SEM and Specfem2D

Many numerical methods can be used to simulate the propagation of seismic waves,including the SEM,finite element method(FEM),finite difference method (FDM),and pseudo-spectral method (PSM).The SEM proposed by Patera (1984) was first used in fluid mechanics calculations,and then gradually applied to the numerical simulation of seismic wave propagation (Seriani et al.,1995; Faccioli et al.,1997; Komatitsch,1997).The SEM combines the ideas of the FEM and PSM and has the advantages of both methods,that is,the flexibility of FEM and the accuracy of PSM.

This paper uses an advanced numerical code,Specfem2D,to explore the dynamic wave interaction between seismic waves and deep tunnels.The Specfem2D software package based on the SEM has two main features,i.e.the tensor-product Lagrange interpolants and Gauss-Lobatto quadrature (Komatitsch and Tromp,1999,2003).Specfem2D is a powerful tool for simulating forward and adjoint seismic wave propagation in 2D acoustic,(an)elastic,poroelastic,or coupled acoustic-(an)elastic-poroelastic media,with convolution perfectly matched layer(PML)absorbing layers as well as higher-order time schemes (Basabe and Sen,2007).It has been proven that Specfem2D has high accuracy and convergence properties for simulating seismic wave propagation and can accurately capture the ground motion around tunnels (Seriani and Oliveira,2008; Lin et al.,2014; Wang and Cai,2015; Pan et al.,2020;Rohnacher et al.,2021).

3.Simulation schemes and numerical models

3.1.Numerical model configuration and layout of receivers

The numerical model of the P-SV system constructed in this paper is shown in Fig.3a.The size of the computation domain is 104 m 104 m,the tunnel is located in the center of the model(X=50 m,Z=-50 m).The four edges of 2 m thickness are set as the PML absorbing layers to eliminate the unreal seismic wave reflection.This study uses the engineering background of the headrace tunnel#2 at Jinping II Hydropower Station.The rock mass in the numerical model is homogeneous Jinping marble (Zhang et al.,2012; Qiu et al.,2014).The Young’s modulus,density,and Poisson’s ratio of Jinping marble are 56 GPa,2780 kg/m3and 0.27,respectively.

Three sets of receivers are placed in the model to capture the PPV distribution around the tunnel,covering the entire simulation area around the tunnel(Fig.3b).The Set#A receivers are arranged in a 25 m 25 m square region centered on the tunnel.The interval between two adjacent receivers is 0.5 m; The Set #B receivers are arranged on the tunnel surface.Since the ground motion on the tunnel surface is relatively strong,the interval between two adjacent receivers in Set#B is 0.25 m.The Set#C receivers are arranged in a 100 m 100 m square region surrounding the Set#A receivers,with an interval of 1 m.Since the simulated tunnels in different schemes have varied shapes and sizes,the location and number of the three-receiver sets in different schemes are slightly modified.

3.2.Simulation schemes

This paper focuses on the control effect of four factors,i.e.the main seismic source wavelength (λs),the tunnel span (D),the tunnel shape,and the damage zone (including the EDZ and EFZ)range on the amplification effect of seismic waves around deep tunnels.A total of five simulation schemes were designed in the study,as shown in Table 2.First,we performed the simulations of seismic wave propagation excited by the near-field finite seismic source.A model without tunnel excavation(the background model)is adopted.Its simulation results are considered background data for the subsequent ground motion parameter analysis,as shown in Scheme #1 in Table 2.Except for non-excavation scheme,each background model’s parameter settings are consistent with the subsequent Schemes #2-#5.Second,the excavation models are simulated for seismic wave propagation,as shown in Schemes#2-#5 in Table 2.Finally,a comparative study of the seismic wave field characteristics of the excavation model with its corresponding background model is carried out to explore the amplification effect of seismic wave field around deep tunnels in consideration of the amplification factor (α).

Fig.2.Calculation flowchart of model parameters for the FFSSM.Ae represents the area of the fracture surface of the seismic sub-source (Ae=Le2 in the FFSSM); MW and MeW represent the moment magnitudes of the main seismic source and a seismic sub-source,respectively;Vf denotes the seismic source fracture velocity,which is usually taken as 0.6-0.85 times of the shear wave velocity (Vf=0.8Vs in the FFSSM); NR represents the number of seismic sub-sources that have ruptured; Di denotes the distance between the ith seismic sub-source (SS #i) and the initially activated sub-source (SS #1),and Di=niLe; s(t) is the standard definition of the second derivative of Gaussian of Ricker wavelet;m(Xs,t) and mij(Xs,t) represents the moment density tensor of the main seismic source and a seismic sub-source,respectively; Mij denotes moment tensor component of the seismic sub-source; δ(X-Xs) denotes the Dirac Delta distribution function of the seismic sub-source; w,u,and σ represent an arbitrary test vector,the displacement vector the stress tensor,respectively; S represents the fracturing surface of seismic source; n represents the unit outward normal to the boundary (Γ); Ω represents an elastic medium; ρ represents the density.These formulas are derived from Berkhout (1987),Faccioli et al.(1997),Beresnev and Atkinson (1999),Somerville et al.(1999),Aki and Richards (2002),Chapman (2004),Motazedian and Atkinson (2005),Carcione (2007),Chaljub et al.(2007),Potvin et al.(2010),and Fichtner (2011).

The seismic source is located in the surrounding rocks at the excavation’s lower-left corner and is simulated with the FFSSM.The seismic sub-sources are set to 3,marked as 1,2,and 3,respectively.Seismic sub-source 1 is located between 2 and 3,as shown in Fig.4.The fracturing direction of the seismic sub-source is specified as a bidirectional fracturing mode,i.e.the fracture propagates from seismic sub-source 1 to 2 and 3 simultaneously.The seismic subsource 1 is activated firstly,and then the follow-up seismic subsources 2 and 3 are activated with a delay time of ti.The simulated main seismic source is a normal fault.The strike φs,rake γrand dip δdof the seismic sub-sources are-90,-90,and 45,respectively.The Ricker time wavelet is used as the time function of seismic sub-sources.The stress drop of the seismic sub-source is set to 5 MPa.Table 2 lists the study content and model configurations for each numerical scheme and the specific schemes are described below.

Fig.3.Numerical model configuration and layout of receivers:(a)Numerical model configuration with spectral element meshes and PML absorbing layers;and(b)Layout of three sets of receivers in simulation region or around the tunnel.

Scheme #2: The tunnel is fixed as a gate arch-shaped tunnel with a tunnel span D=12 m.The main seismic source wavelength is changed by controlling the length of the seismic sub-source(Fig.4a).The seismic source is located in the rock masses below the lower-left corner of the tunnel.Seismic sub-source 1 is located at the point of X=30 m,Z=-70 m and 28.28 m from the center of the tunnel.

The coordinates of seismic sub-sources 2 and 3 can be determined by the seismic sub-source length(Le).Scheme#2 includes 11 cases with different Levalues (Le=0.5 m,0.75 m,1 m,1.5 m,2 m,4 m,6 m,9 m,12 m,15 m,and 25 m).Table 3 lists the source parameters relating to each case.Each parameter is determined according to Path#1 (Fig.2).

Scheme #3: Keeping the source parameters constant(Le=1.5 m,same as Case 4 in Scheme #2 in Table 3).The wavelength of the S-wave is therefore fixed at 3 m.The cross-sectional shape of the tunnel remains circular,and its spans in 7 cases in Scheme#3 are set at 1.5 m,3 m,6 m,10 m,12 m,16 m,and 20 m,respectively (see Fig.4b,Table 2).The change in tunnel span (D),implies a change in the interaction range between seismic waves and the tunnel section,leading to changes in reflection,diffraction,and refraction of seismic waves,ultimately affecting the site effect of ground motion.

Scheme#4:The seismic source parameters and the tunnel span remain unchanged,and only the tunnel shape is changed.Scheme#4 includes two sets of cases (Set 1 and 2),and each set includes four cases(see Fig.4c,Table 2).The cases in the same set have the same seismic source parameters(4 m and 6 m for Leof Set 1 and Set 2,respectively;see Table 3 for cases Nos.6 and 7)and tunnel span(D=12 m),but have different tunnel shapes (circular,gate archshaped,rectangular,and arched tunnels for cases Nos.1-4,respectively).Different tunnel shapes lead to different directions of reflection and diffraction of seismic waves,changing the characteristics of the seismic wave field and thus resulting in variations in the magnification and distribution of the amplification factor (α).

Scheme#5:The seismic source and tunnel parameters are fixed in this scheme,and only the range of damage zones is changed.In the simulation,the surrounding rocks are subdivided into the EFZ,EDZ,excavation influence zone(EIZ),and original rock zone(ORZ).Numerous studies show that the excavation damage in overstressed rock masses decreases sequentially from EFZ to ORZ until it disappears (Siren et al.,2015; Bao et al.,2020).The corresponding rock wave velocity gradually increases to the original rock wave velocity.The elastic wave velocities are assigned in each zone to simplify the simulation.In the EFZ and EDZ,the elastic wave velocities of rock masses are taken as 0.4 and 0.7 times of the elastic wave velocity of original rock masses in the ORZ,respectively.In particular,the elastic wave velocity of rock masses in the EIZ matches with the ORZ.According to the acoustic test results conducted in headrace tunnel#2 at Jinping II Hydropower Station,the rock masses in the ORZ,EIZ,EDZ,and EFZ are given different elastic wave velocities listed in Table 4.

Moreover,as the degree of damage increases,the range of the damage zones increases,as shown in Fig.4d.Scheme #5 includes two sets (Set 1 and Set 2),each including four cases (as shown in Fig.4d and listed in Table 2).The cases in the same set have the same seismic source parameters(4 m and 6 m for Leof Set 1 and Set 2,respectively),tunnel span of D=11.67 m and different EDZ and EFZ ranges.The ratio parameter k is used to control the regions of EDZ and EFZ.In this scheme,k is set to 0,0.5,0.8,and 1,respectively,as shown in Fig.4d.

Fig.4.Layouts of seismic sub-sources in the proposed FFSSM for different parametric analysis schemes:(a)Scheme#2(Cases 5-9);(b)Scheme#3(Cases 2-6);(c)Scheme#4;and(d)Scheme#5.Not all cases are drawn in the figure due to the limited size.However,only some typical cases in Table 2 are selected to illustrate the changing law of the parameters of different schemes.

4.Model results

4.1.Influence of the main seismic source wavelength

Schemes#1 and#2 in Table 2 were implemented to explore the effect of the main seismic source wavelength on the characteristics of the seismic wave field and the distribution of strong ground motions,which can be formed by the dynamic interaction of seismic waves with tunnels.This paper focuses only on the interaction between S-waves with tunnels,mainly because S-waves transmit more energy and can cause a more significant ground motion response in the surrounding rocks of tunnels,while the Pwave induced ground motion is much smaller(Hashash et al.,2001;Aki and Richards,2002; Wang and Cai,2015).

Fig.5 shows snapshots of the vertical velocity components in the seismic wave field for seismic sub-source lengths of 0.5 m,1.5 m,and 4 m at three-time points.The moment tArepresents the time when the leading edge of the S-wave first reaches the left foot of the tunnel.tBis the moment when the S-wave fully interacts with the tunnel.tCis the time when the S-wave passes through the right arch of the tunnel.It should be noted that in Fig.5,the specific values of the above three moments are different,which is caused by the propagation change of the induced seismic waves due to the difference in the wavelength of the seismic sub-source (Le) in the FFSSM.It can be seen from Fig.5 that (1) As the λsin the FFSSM increases(Table 3),the interaction times(tC-tA)between S-waves and tunnel increase accordingly.When Leis 0.5 m,1.5 m,and 4 m,the interaction times are 5.2 ms,8.8 ms,and 14.4 ms,respectively.(2)The increase in wavelength increases the delay time of seismic waves passing through the tunnel.However,more significantly,the seismic wave field formed by the interaction of seismic waves with the tunnel is more complex,as shown in Fig.5c.The range of their interaction is large.(3) When the maximum wave front of seismic S-wave propagates to the tunnel,the reflected wave of S-wave is generated on the left arch shoulder side and the right bottom corner side.It acts jointly with the original incident S-wave and its diffracted wave to form a seismic wave superimposed area (white dashed ovals in Fig.5),where the ground motion is superimposed and enhanced.This is quite different from the background model’s distribution of strong ground motions.In the background model,the maximum seismic wave response region is mainly along the two directions shown by the white dashed lines in Fig.5,i.e.the two directions with counterclockwise intersection angles of 45 and 135 with the X-axis.The propagation direction of the FFSSMinduced seismic waves changed due to the obstruction from the tunnel,which results in the fact that the area of the maximum wave front in the background model is weakened in the excavation model to form a wave weakened area(green dashed oval in Fig.5).An enhanced area is formed on either side of the weakened wave area.

Table 2 Simulation schemes for exploring influencing factors of the near-field ground motion amplification effect around deep tunnels based on the FFSSM.

Table 3 List of seismic sub-source parameters of different seismic sub-source lengths used in Scheme #2.

Table 4 Parameters of different surrounding rock zones for Scheme #5.

Fig.6 shows the PPV contours and amplification factor distributions for Le=0.5 m,1.5 m,4 m,6 m,and 9 m.Fig.6a and b shows that PPVs increase with increasing Lefor both the background model and the excavation model.The zones with high-PPV values in the background model are mainly distributed along the line connecting the center of the tunnel and the center of the FFSSM seismic source.The PPVs in the background model also follow the pattern of decay with the epicenter distance (i.e.a distance to the center of the FFSSM seismic source),which show the maximum at the FFSSM seismic source and decrease gradually with the increase of the epicenter distance.In the excavation model,high-PPV zones(Zone-A and Zone-B) are formed due to the generation and propagation of reflected S-waves and diffracted S-waves.Zone-A is located within a certain depth range of the surrounding rock on the tunnel’s left sidewall and left floor.It is formed by the superposition of the reflected S-waves from the initial fracturing seismic subsource 1 and the incident S-waves of the follow-up activated seismic sub-sources 2 and 3.As shown in Fig.6b,the extent of Zone-A increases with the increase of Leand gradually expands from the left sidewall and floor area to the left arch shoulder and right floor area.Zone-B appears far away from the left wall and floor of the tunnel,which is formed by the joint superposition of reflected S-wave,incident S-wave and diffracted S-wave from the FFSSM seismic sub-sources.Additionally,the shielding effect of the tunnel and diffraction loss of seismic waves lead to the formation of a low-PPV zone,Zone-C,which locates on the right arch shoulder side of the tunnel,and a striped high-PPV zone,Zone-D,which is formed on both edges of Zone-C due to diffraction waves.As the FFSSM seismic sub-source length(Le)increases,the extent of Zone-C increases gradually,and Zone-B intermingles gradually with Zone-D.Moreover,Fig.6c and d shows the effect of the main seismic source wavelength on the distribution of PPV amplification factor.Two types of PPV characteristic zones are formed: the PPV amplified zone(α ＞1,as shown in Fig.6c and d for A1,A2 and A3,where A1,A2,and A3 are PPV amplified areas),and the PPV weakened zone (α ＜1),which reflects the blocking and shielding effect of the tunnel on ground motion.These two PPV characteristic zones correspond well to the high-and low-PPV zones in Fig.6b.As the seismic sub-source length (Le) increases,zone A1 shifts gradually from the left sidewall and left floor corner areas of the tunnel to the left arch shoulder and right floor corner areas.Meanwhile,zones A2 and A3 also approach gradually each other and finally merge into a PPV amplification zone.This is caused by the change of the maximum wavefront coverage area of seismic waves and their interaction region with the tunnel,accompanied with the change of the interaction of multiple waves(i.e.incident S-wave,reflected Swave,and diffracted S-wave of the FFSSM seismic sub-sources).The former is due to the increase in the seismic sub-source length,which increases the length of the main seismic source.The latter is caused by the increase in the wavelength of the main seismic source.The increase in the main seismic source wavelength changes the excitation time of the FFSSM seismic sub-source.It causes the phase difference between the seismic sub-source waves,thus increasing the interaction time between the seismic waves and the tunnel (Table 3),and changing the multi-wave superposition effect.

Fig.5.Vertical velocity fields for FFSSM with different seismic sub-source lengths recorded at three-time points:(a)Le=0.5 m;(b)Le=1.5 m;and(c)Le=4 m.The colors of the seismic waves represent the positive or negative of the vertical velocity components,red represents positive while blue represents negative; the white dashed ovals and green dashed ovals represent wave superimposed areas and wave weakened areas,respectively.

The maximum value of PPV,PPVm,is one of the most important parameters.For this,the relationship between the PPVm,the maximum amplification factor (αm) and the main seismic source wavelength(λs)in the near-tunnel area are given in Fig.7a and b.It is worth noting that the main focus here is on the ground motion characteristics of the rock masses in a certain depth range around the tunnels.Firstly,most seismic hazards often occur within the shallow rock masses near the tunnel boundary (Feng et al.,2012;He et al.,2018).Secondly,as revealed from Fig.6,the PPV and amplification factor of the near-tunnel zone A1 are high.Compared with zones A2 and A3,the ground motion in zone A1 is much stronger.As shown in Fig.7a,as λsincreases from 1.05 m to 45.12 m,PPVmincreases from 0.16 cm/s to 30.72 cm/s in the excavation model and from 0.14 cm/s to 21.15 cm/s in the corresponding background model.In Fig.7b,it can be seen that αmcan reach around 2.13.Regularly,as λsincreases,the α shows a trend of first increasing and then decreasing.Especially,the PPV amplification effect is the strongest at a specific S-wave wavelength,such as λs=20 m in Fig.7b.

Fig.6.PPV contours and amplification factor distributions for Le=0.5 m,1.5 m,4 m,6 m and 9 m:(a) PPV contours in the background model;(b)PPV contours in the excavation model;(c)Amplification factor distributions around the tunnel(25 m 25 m);and(d)Amplification factor distributions for the entire computation domain(100 m 100 m).Zone-A,Zone-B,and Zone-D are the high-value zones of PPV.Zone-C is the low-value zone of PPV.A1,A2,and A3 are PPV amplified zones.

Fig.7.The relationship between the near-field ground motion in the near-tunnel area and the main seismic source wavelength: (a)The variation of the maximum PPV with the main seismic source wavelength; and (b) The variation of the maximum amplification factors with the main seismic source wavelength.

4.2.Influence of the tunnel span

This section focuses on the effect of the tunnel span (D) on the ground motion at a fixed wavelength(Le=1.5 m,λs=3 m).Seven tunnel spans(D=1.5 m,3 m,6 m,10 m,12 m,16 m,and 20 m)are chosen to perform the simulations.Details of the scheme are discussed in Scheme #3 in Section 3.2.

For tunnel spans of 6 m,10 m,and 16 m,Fig.8 illustrates the distribution of the vertical velocity components in the seismic wave field at four-time points.It can be seen that the change in the tunnel span can induce significantly different seismic wave fields despite the constant seismic source parameters.As D increases,the time for the seismic wave to reach the lower-left boundary of the tunnel is earlier,and the time to leave the upper-right boundary of the tunnel is delayed,which means that the time difference (tCtA)for the S-waves to interact with the tunnel structure increases.The larger D,the shorter the time for the seismic S-wave to reach the lower-left boundary of the tunnel,and the later it leaves the upper-right boundary of the tunnel.It means that the time difference of S-wave interacts with the tunnel (denoted by tC-tAin Fig.8)increases with the increase of tunnel span.From Fig.8,when D is 6 m,10 m,and 16 m,tC-tAare 4.8 ms,6.4 ms,and 8.4 ms,respectively.Moreover,as D increases,the seismic wave field around the tunnel also has become more complex,such as the appearance of stronger scattered and diffracted waves.

For D=1.5 m,3 m,6 m,10 m,and 16 m,the distributions of PPV and amplification factor of ground motion around the tunnel are illustrated in Fig.9.The corresponding PPV in the background model can be found in Fig.6a(Le=1.5 m)because the FFSSM source parameters used in Figs.6a and 9 are the same.In Fig.9a,three high-PPV zones,shown as Zones-A,B,and D,and one low-PPV zone,Zone-C,could be found.Moreover,the PPV amplification zones are marked as A1,A2,and A3,as demonstrated in Fig.9b and c.Similarly,a PPV weakened zone can be found in the tunnel upperright area.For the case of Le=1.5 shown in Fig.9a,analysis of the PPV values for different tunnel spans shows that as D increases,the PPV values also increase accordingly.For the PPVmillustrated in Fig.10a,D increases from 6 m to 16 m,and the corresponding PPVmincreases from 0.65 m/s to 0.95 m/s.The PPVs in the excavation model increase significantly relative to those in the background model of Fig.6a,as indicated in Fig.10a,and the site effect of the tunnel is the major reason for the seismic wave amplification effect.In Fig.10a,when D increases from 2 m to 20 m,the PPVmin the background model increases by approximately 20%,while the PPVmin the excavation model increases by approximately 50%.The amplification factor in Fig.10b gives a more visual indication of the amplitude,positively correlated with the increase of D.The αmincreases from 1.61 to 1.97 when D increases from 2 m to 16 m.However,as D increases,the αmshows a trend of first increasing rapidly and then leveling off,as shown in Fig.10b.Additionally,for the PPV weakened zone shown in Fig.9b and c,it can be found that the weakened zone increases significantly with increasing tunnel span because a larger tunnel span can lead to a longer tunnel boundary to impede seismic wave propagation.

4.3.Influence of tunnel shape

The circular tunnel and the seismic source parameters used in Section 4.2 result in symmetry of the PPV distribution field.The symmetry axis of PPV distribution is the white dashed line shown in Fig.8,which is a special scenario.However,by comparing Figs.6 and 9,it can be observed that the symmetry of the PPV distribution disappears for the gate arch-shaped tunnel.Thus,the tunnel shape controls the seismic wave field and the near-field ground motion amplification effect.To explore the influence of tunnel shape on the near-field ground motion amplification effect,Scheme #1 and Scheme #4 in Table 2 are implemented.

Fig.11 shows snapshots of the vertical velocity component at 15.2 ms of four cases with different tunnel shapes in Set 1(Le=4 m)of Scheme #4.Under the conditions of the same FFSSM seismic source location,the variability of the four tunnel shapes slightly changes the distribution pattern of the seismic wave field.This is mainly reflected in the slight difference in the interaction processes between the seismic waves and the tunnel boundary,making the composite seismic wave field among the incident wave,diffracted wave,and reflected wave of the seismic wave different.However,Fig.11 shows that the four tunnel shapes have a limited degree of alteration to the seismic wave field.Therefore,we further analyze the PPV and its amplification factor distribution under four tunnel shapes,as shown in Fig.12.Different tunnel shapes possess roughly the same PPV range (0.6-3.4 cm/s) and similar distribution patterns (Fig.12a).Under the conditions of the FFSSM seismic source location and parameters set in Scheme#4,the PPV distribution has a good symmetry in circular and rectangular tunnels,approximate symmetry in gate arch-shaped tunnels,and the asymmetry in arched tunnels.For example,the PPVs at the floor of arched tunnels are higher compared to those at the left wall.Compared with the other three tunnels,the PPVmof the circular tunnel is slightly larger,approaching 3.4 cm/s; whilst the other three values are closer,about 3 cm/s.This difference is more sound in the amplification factor distribution shown in Fig.12b and c.The amplification factor distributions of the circular and rectangular tunnels are similarly symmetric about the direction of seismic wave incidence.Those of the gate arch-shaped tunnels are also approximately symmetric.Conversely,arched tunnels have some variability in zones A1 and A2+A3.The αmof zone A1 on the left floor of the arched tunnel is 1.7,while that of the left wall is 1.6.However,in high-PPV regions,αmis reached to 4.2 for the right-lower region and 2.4 for the leftupper region of the arched tunnel.

Fig.8.Vertical velocity fields for the circular tunnel with different tunnel spans recorded at three-time points:(a)D=6 m;(b)D=10 m;and(c)D=16 m.The colors of the seismic waves represent the positive or negative of the vertical velocity components,red represents positive while blue represents negative;the white dashed ovals,and green dashed ovals represent wave superimposed areas and wave weakened areas,respectively.

The above evidence shows that the tunnel shape influences the site effect of seismic waves,which can dominate the location of ground motion amplification and weakened zones and affect the degree of ground motion amplification to some extent.However,in Fig.12,the effect of tunnel shape on the ground motions in the near-tunnel rock masses is relatively minor in terms of the change of PPVs and amplification factors.Fig.13 shows the maximum amplification factors of tunnel-induced ground motions for four tunnel shapes with FFSSM seismic sub-source lengths of 4 m(Set 1)and 6 m(Set 2).For four tunnel shapes,the maximum amplification factors,αm,are 1.87 and 2.04 for Sets 1 and 2 of circular tunnels,and 1.84 and 1.99 for rectangular tunnels,1.77 and 1.94 for gate archshaped tunnels,and 1.75 and 1.79 for arched tunnels,respectively.Taking the amplification factor of circular tunnels as the base,the amplification factor of the rectangular tunnel,gate arch-shaped tunnel,and the arched tunnel is 0.98,0.96,and 0.9 times of the αmof the circular tunnel,respectively,which is calculated by the ratios of the average amplification factor of Sets 1 and 2.

Fig.9.PPV contours and amplification factor distributions for D=1.5 m,3 m,6 m,10 m,16 m: (a) PPV contours for the excavation model; (b) Amplification factor distributions around the tunnel(25 m 25 m);and(c)Amplification factor distributions for the entire computation domain(100 m 100 m),in which Zone-A,Zone-B,and Zone-D are the highvalue zones of PPV.Zone-C is the low-value zone of PPV.A1,A2,and A3 are PPV amplified zones.

Fig.10.The relationship curves between the near-field ground motion in the neartunnel area and the tunnel span: (a) The variation of maximum PPVs with the tunnel span;and(b)The variation of maximum amplification factors with the tunnel span.

4.4.Influence of the range of damage zones

To understand the influence of the range of damage zones,the simulations of Schemes #1 and #5 in Table 2 are conducted.Four ranges of damage zones are selected as shown in cases 1-4 of Fig.4,set by a damage range parameter,k.

For four cases with k=0,0.5,0.8,and 1,the vertical velocity distributions at 16.8 ms in Set 1 of Scheme#5 are shown in Fig.14.As the range of EDZ and EFZ of the rock masses around the tunnel increases,that is,the k value increases,the site effect caused by the interaction of seismic waves is enhanced accordingly.In Fig.14,various diffractions,refractions,and reflections of the incident waves at the EDZ,EFZ,and tunnel borders make this site effect more noticeable on both sides of the left arch shoulder and the right floor of the tunnel.The interaction of these multiple seismic sub-waves changes the location and extent of the appearance of strong ground motion and the degree of its amplification.Additionally,it can be found from Fig.14 that the change of k can result in different tunnel shielding effects.With the increase of k,the ground motion on the right arch shoulder side of the tunnel tends to weaken,implying the gradual enhancement of the shielding effect.

Fig.15 illustrates the distribution of PPV contours and amplification factors corresponding to the conditions of Fig.14.As the k increases,PPVmin the Zone-A increases from 3 cm/s to 6.2 cm/s,and the extent of PPVmgradually extends from the left floor to the left sidewall of the tunnel.Additionally,with the increase of k,several small high-PPV regions also appear gradually in the right arch shoulder and right wall areas of the tunnel,where their PPVmvalues also increase progressively and can reach 5 cm/s when k is 1.In contrast to Zone-A,the range of Zone-B and Zone-D decreases as k increases.For the PPV amplification factors as shown in Fig.15b and c,with the increase of k from 0 to 1,the αmin the zone A1 increases gradually from 1.7 to 3.7.The range of A1 extends from the left floor region to the left sidewall region.Meanwhile,the αmincreases progressively to 3.3 in both the tunnel’s right shoulder and sidewall regions.However,in the other two PPV amplification zones,A2 and A3,the αmdecreases from 4.2 to 2.1 as k increases,whose range also decreases.Moreover,with the increase of k,the αmis reduced from 0.4 to 0.2 in the PPV weakened zone.In short,combining the results of Figs.14 and 15,the increase of the extents of EDZ and EFZ can increase ground motions in the near-tunnel areas,which is manifested by the increase in PPVmand αm.On the contrary,it also decreases ground motions in the far-tunnel and tunnel shielding areas,where the shielding effect of the tunnel increases significantly,as illustrated in the PPV weakened zones in Fig.15b and c.

Considering four ranges of damage zones and two conditions of the FFSSM seismic sub-source length in Scheme#5,Fig.16 shows the αminthe near-tunnelareas.Taking the excavationmodelwithk=0as a reference,when k is 0.5,0.8,and 1,the αmvalues are 1.4,1.8,and 2 times that of k=0,respectively.These relative times are determined basedontheaverageamplificationtimesfortheFFSSMseismicsource lengths of 4 m(Set 1 in Fig.16)and 6 m(Set 2 in Fig.16).

5.Empirical formula of the amplification factor in the neartunnel area

Based on the above results of the control effects of the four factors (e.g.the main seismic source wavelength,tunnel span,tunnel shape,and range of damage zones) on the amplification of near-field ground motion in the near-tunnel area,an empirical formula for αmcan be proposed,as a product similar to that used in the Q-system for classifying rock masses:

5.1.The wavelength control factor,Fλ

According to the relation between αmand λsillustrated in Fig.7b in Section 4.1,a best-fit curve of a Gaussian function between αmand λsis proposed using the least square method,as shown in Fig.17a.The wavelength control factor,Fλ,is established through a formula expressed by Eq.(2).The coefficients in Eq.(2) are finally obtained by the least-squares iteration,and the minimum fit converges error optimization.

5.2.The tunnel span factor,FD

The tunnel span factor (FD) is introduced to quantify the contribution of the tunnel span to the increase in near-tunnel ground motion.It can be found that Eq.(3) has a typical logistic functional form,and the best-fit curve between αmand D given by Eq.(3)is illustrated in Fig.17b.According to the analysis in Section 4.2,αmshows a tendency to increase with increase of D,and the increasing trend decreases gradually.

5.3.The tunnel shape factor,Fs

Fig.11.Snapshots of the vertical velocity component at 15.2 ms of four cases with different tunnel shapes in Set 1 of Scheme#4:(a)Circular tunnel;(b)Gate arch-shaped tunnel;(c)Rectangular tunnel;and(d)Arched tunnel.The colors of the seismic waves represent the positive or negative of the vertical velocity components,red represents positive while blue represents negative; the white dashed ovals and green dashed ovals represent wave superimposed areas and wave weakened areas.

Similarly,the tunnel shape factor,Fs,is introduced to quantify the contribution of the tunnel shape to the increase in near-tunnel ground motion.However,unlike the Fλ and FDthat are established by using the best-fit relationship between the maximum amplification factor and its control factors (λsand D),the Fsis quantified using a rating scale,as listed in Table 5.In Table 5,only the four most used tunnel shapes are considered and given separate scores,with circular tunnels selected as a baseline for scoring,i.e.Fs=1 for a circular tunnel.The scoring values for the other three tunnel shapes are given based on the laws obtained in Section 4.3.As listed in Table 5,the values of Fsfor four tunnel shapes ranged from 0.9 to 1,which illustrates that the tunnel shape has a relatively small effect on the near-tunnel ground motion amplification effect.

5.4.The excavation damage factor,Fd

The excavation damage factor,Fdis shown in Table 6.From the simulation results in Section 4.4,the range of the damage zones has a significant impact on αm.As the ratio parameter k increases,that is,the range of EDZ and EFZ increase,the amplification effect of ground motion in the near-tunnel area is significantly enhanced.For the convenience of the Fdvalues,the tunnel damage is classified into five levels according to the maximum excavation damage depth (the maximum depth of the EDZ,dm): no,weak,medium,strong,and extremely strong.The αmvalues for different damage levels are interpolated based on the simulation results in Section 4.4.In Table 6,for very strong excavation damage tunnels,the recommended value of Fdis 2.5,but the value can be higher than 2.5 in special cases.However,for deep hard rock tunnels,the engineering conditions where the depth of EDZ is greater than 5 m are relatively uncommon.The Fdvalue of 2.5 has given sufficient consideration to the ground motion amplification effect.

Fig.12.PPV contours and amplification factor distributions for the circular tunnel,gate arch-shaped tunnel,rectangular tunnel,and arched tunnel in Set 1:(a)PPV contours for the excavation model; (b) Amplification factor distributions around the tunnel (25 m 25 m); and (c) Amplification factor distributions for the entire computation domain(100 100 m2).Zone-A,Zone-B,and Zone-D are the high-value zones of PPV.Zone-C is the low-value zone of PPV.A1,A2,and A3 are PPV amplified zones.

Fig.13.The maximum amplification factor of the near-tunnel area around deep tunnel with different tunnel shapes in Scheme #4.Cases 1 to 4 are for circular tunnel,gate arch-shaped tunnel,rectangular tunnel,and arched tunnel,respectively.

5.5.Applications

To verify the reliability of the empirical formula,the αmfitting value obtained through Eq.(1)is compared with the recorded value obtained by the receiver in Schemes #2-5 (see Table A1).The αmfitting values of the empirical formula are very close to the αmrecorded values of the seismic wave propagation simulation based on the FFSSM.Among the 34 cases in the four schemes,the relative error (δ) is less than 5% in 21 cases,the δ value is between 5% and 10%in 11 cases,and the δ value exceeds 10%(14.14%and 17.53%)in only two cases.Among these 34 cases,the proportion of cases with relative errors within 5%is as high as 61.8%.The proportion of cases with relative errors within 10% is nearly 94.1%.Therefore,for the cases in this paper,the proposed empirical formula can estimate the αmin the near-tunnel area very well.

Fig.14.Snapshots of the vertical velocity component at 16.8 ms of four cases with different k values in Set 1 of Scheme#5:(a)k=0;(b)k=0.5;(c)k=0.8;and(d)k=1.The colors of the seismic waves represent the positive or negative of the vertical velocity components,red represents positive while blue represents negative; the white dashed ovals,and green dashed ovals represent wave superimposed areas and wave weakened areas,respectively;k refers to the ratio of the modeled EDZ and EFZ of the tunnel to the actual EDZ and EFZ.

Furthermore,for six published typical deep tunnel projects,the estimated αmvalues from Eq.(1) in this paper and the measured values from the field monitoring under nine cases were further compared,as shown in Table 7.Those estimated αmvalues in nine cases match the measured values,suggesting that our proposed empirical formula can be applied to practical engineering with high accuracy.Nevertheless,due to few studies on the ground motion amplification effect in the near field of the tunnel,some parameters in Eq.(1) are adopted from engineering experiences when calculating the maximum amplification factor in Table 7.This may affect the accuracy of the proposed empirical formula.However,in any case,the empirical formula proposed in this paper is very practical and simple for estimating the amplification effect of near-field tunnel ground motion.

Table 5 List of the tunnel shape factor,Fs.

Table 6 List of the excavation damage factor,Fd.

Table 7 Analysis table of relative error between the fitting values and measured values of αm in six underground engineering projects.

Fig.15.PPV contours and amplification factor distributions of four cases with different k values in Set 1 of Scheme#5:(a)PPV contours for the excavation model;(b)Amplification factor distributions around the tunnel(25 m 25 m);and(c)Amplification factor distributions for the entire computation domain(100 m 100 m).Zone-A,Zone-B,and Zone-D are the high-value zones of PPV.Zone-C is the low-value zone of PPV.A1,A2,and A3 are PPV amplified zones.

6.Conclusions

Fig.16.The maximum amplification factor in the near-tunnel area with different EDZ and EFZ ranges in Scheme #5.

Fig.17.The fitting relationship between the amplification effect in the near-tunnel area and λs and D (a) Fitting curve of αm and λs (b) Fitting curve of αm and D.

This paper studied the characteristics and distributions of strong ground motions caused by the interaction between near-field seismic waves and deep tunnels based on the proposed FFSSM.The control effects of the main seismic source wavelength,the tunnel span,the tunnel shape,and the range of damage zones on the near-field ground motion amplification behaviors are revealed.Then,an empirical formula in the form of αm=FλFDFsFdthat can be used to estimate the maximum amplification factor of near-field ground motions in the near-tunnel area is proposed.The conclusions are made as follows.

(1) The main seismic wavelength (λs) plays a crucial role in controlling the ground motion behavior around the tunnel,which affects the interaction time of seismic waves with the tunnel and changes the distribution of PPVs and their amplification factors.The interaction time and PPVs increase with increasing λs,while the maximum amplification factor αmin the near-tunnel area shows a trend of first increasing and then decreasing.Additionally,several PPV amplification zones are formed in the interaction regions of the incident waves and the secondary reflected waves.Conversely,the PPV weakened zones are formed on the back side of the tunnel in the direction of the incident seismic waves.

(2) The tunnel span (D) also affects ground motion around the tunnels.When the tunnel span is changed,even if the source parameters of the FFSSM are constant,the interaction time of seismic wave with the tunnels and the induced strong ground motions in the near-tunnel area is changed.As D increases,the PPVm,and its corresponding αmincrease.

(3) The tunnel shape slightly affects the propagation of seismic waves radiated by the FFSSM seismic source and its strong ground motion distribution in the near-tunnel area.A study of four common tunnel shapes (circular,gate arch-shaped,rectangular,and arch) shows that circular tunnels have the largest maximum PPV amplification factors.

(4) The ranges of the EDZ and EFZ of the surrounding rocks induced by the high-stress unloading of the deep tunnels have a significant impact on the near-field seismic wave propagation and its strong ground motion distributions around deep tunnels.As the EDZ and EFZ ranges increase,the PPV and its amplification factor distribution change,while the PPVmand αmalso increase significantly.

The applicability and accuracy of the empirical formula are verified.This empirical formula provides an easy approach for accurately estimating the ground motion parameters in seismic hazard risk evaluation (e.g.fault-slip rockburst and seismicinduced collapse) and the rock support design of the deep tunnels under dynamic load conditions.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was jointly supported by the National Natural Science Foundation of China (Grant No.41877256),the Natural Science Foundation of Hubei Province (Grant No.ZRQT2020000114),the Key Research Program of the Chinese Academy of Sciences (Grant No.KFZD-SW-423).

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.12.024.