A numerical model for ground-borne vibrations from underground railway traffic based on a periodic finite element–boundary element formulation

A numerical model is presented to predict vibrations in the free field from excitation due to metro trains in tunnels. The three-dimensional dynamic tunnel–soil interaction problem is solved with a subdomain formulation, using a finite element formulation for the tunnel and a boundary element method for the soil. The periodicity of the geometry in the longitudinal direction of the tunnel is exploited using the Floquet transform, limiting the discretization to a single-bounded reference cell. The responses of two different types of tunnel due to a harmonic load on the tunnel invert are compared, both in the frequency–wavenumber and spatial domains. The first tunnel is a shallow cut-and-cover masonry tunnel on the Paris metro network, embedded in layers of sand, while the second tunnel is a deep bored tunnel of London Underground, with a cast iron lining and embedded in the London clay.


Introduction
Underground trains create vibrations which are transmitted through the soil and interact with the foundations of adjacent buildings, resulting in disturbance from vibrations (1-80 Hz) and re-radiated noise .
Numerical models with varying degrees of sophistication are under development to predict vibrations from trains running in tunnels. This work has benefited to a large extent from recently developed models for vibrations due to trains running at grade on discretely supported and continuous tracks.
Dynamic axle loads are mainly caused by wheel and rail roughness, impact excitation due to discontinuities of the wheels and the rails and parametric excitation due to the periodic support of the rails on sleepers. Stateof-the-art modelling of vibrations due to trains in tunnels assumes that the tunnel invert and the surrounding soil are relatively stiff with respect to the track, resulting in a two-step approach where the dynamic axle forces on the rails are derived from an uncoupled track model and subsequently applied on a track-tunnel-soil 1 model to compute the radiation of waves from the tunnel into the soil. The objective of the present paper is to focus on the second problem.
Two-dimensional (plane strain) finite element (FE) models have been used to solve the dynamic (track-)tunnel-soil interaction model, where the infinite extension of the soil domain is represented by local absorbing boundary conditions [1][2][3]. Alternatively, two-dimensional coupled finite element-boundary element (FE-BE) formulations have been used [4,5], where the discretization effort is limited to the finite element mesh of the tunnel and the boundary element (BE) mesh, which can be restricted to the interface between the tunnel and the soil if the Green's functions of a layered half-space are used as an alternative to the full space solutions. However, two-dimensional models necessitate important simplifications to translate the three-dimensional (moving) load into an equivalent line load and underestimate radiation damping into the soil.
Three-dimensional FE models with local absorbing boundary conditions become prohibitively large for applications in the wide frequency range under consideration, as they require an extended finite element mesh together with a small element size. BE formulations partially solve this problem as they reduce the discretization effort. Stamos and Beskos [6,7] have used a three-dimensional BE formulation for the seismic analysis of large underground structures embedded in a homogeneous half-space. For applications to railwayinduced vibrations in tunnels, however, involving the analysis at higher frequencies, three-dimensional coupled FE-BE formulations would also become very expensive from the computational point of view due to the large dimension of the fully populated BE matrices resulting from the coupling of all degrees of freedom on the interface between the tunnel and the soil. This has been the main motivation for many researchers to exploit the longitudinal invariance or periodicity of the track-tunnel-soil system for trains running in tunnels or the track-soil system for trains running at grade.
Forrest [8], Hussein and Hunt [9,10] and Hussein [11] have developed a semi-analytical insertion loss model, that accounts for the essential three-dimensional dynamic interaction between the track, the tunnel and the soil, as well as for wave propagation in the soil. The tunnel is modelled as a thin cylindrical shell surrounded by a homogeneous elastic medium. The model allows for the analysis of deep-bored tunnels embedded in a homogeneous soil. The equilibrium equations are formulated in the wavenumber-frequency domain, exploiting the invariance of the geometry along the longitudinal axis of the tunnel; models of this type are therefore referred to as two-and-a-half dimensional (2.5D), although it must be understood that they capture the three-dimensional behaviour of the invariant tunnel-track-soil system. Hussein [11] has also studied the vibration isolation efficiency of a floating slab in the tunnel, accounting for three different arrangements of the slab bearings.
For trains running at the surface, Sheng et al. [12,13], Madshus and Kaynia [14] and Lombaert et al. [15][16][17] have developed 2.5D models in the wavenumber-frequency domain that account for the invariance of the track in the longitudinal direction. These models are used to compute train-induced vibrations due to moving axle loads; the model of Lombaert et al. also accounts for the interaction between the train and the track.
In the context of the seismic analysis of sheet pile walls of infinite length, considerable progress has been made regarding the solution of dynamic soil-structure interaction problems defined on periodic domains, using coupled periodic FE-BE formulations [18][19][20]. Within the frame of the EC-Growth project CONVURT [21], this methodology has been used to develop a modular numerical model to predict vibrations from excitation due to metro trains in tunnels for both newly built and existing situations [22]. It is assumed that the computation of the dynamic axle loads on the rail can be decoupled from the solution of the dynamic track-tunnel-soil interaction problem. The present paper concentrates on the three-dimensional dynamic track-tunnel-soil interaction problem that is solved with a subdomain formulation, using a finite element formulation for the tunnel and a boundary element method for the soil. The periodicity of the tunnel and the soil is exploited using the Floquet transform, limiting the discretization to a single bounded reference cell of the tunnel. It has been demonstrated that the proposed approach can reuse the existing three-dimensional BE technology for layered media as the periodic Green's kernels have the same singularities as the threedimensional Green's kernels [20]. This is a major advantage compared to the so-called 2.5D approaches, where a translation invariant model is assumed and for which the analysis of all singularities has to be repeated after the Fourier transform along the tunnel axis, except for surface tracks where the singular stress kernel vanishes [23].
The model will be validated by means of in situ experiments that have been performed at a site in the CiteÚ niversitaire on the line Re´seau Express Re´gional (RER B) of Re´gie Autonome des Transports Parisiens (RATP) in Paris [24][25][26][27] and a site in Regent's Park on the Bakerloo line of London Underground [28,29]. The tunnel in Paris is a shallow cut-and-cover masonry tunnel with two tracks, embedded in layers of sand, gravel and marl, while the tunnel in London is a deep cut-and-cover tunnel with a cast iron lining and a single track, embedded in a homogeneous layer of London clay.
After a brief review of the governing system of equations, details on the geometry and construction of both tunnels are presented in the present paper (the reader is referred to other papers [18,20,22] for a complete account of the formulation and solution of periodic dynamic soil-structure interaction problems). The response of both tunnels (and the surrounding soil) due to a harmonic load on the tunnel invert is considered, allowing to draw conclusions on the dynamic behaviour of both tunnel-soil systems. It is also demonstrated how the vibration isolation efficiency of a floating slab track can efficiently be computed using a Craig-Bampton substructuring technique [30]. Validation of the numerical model using the experimental results will be considered in a sequel publication.
The response for a harmonic load at different positions along the longitudinal axis in the reference cell (or the transfer functions) is needed to compute the response due to moving loads. This methodology has been outlined in general terms by Clouteau et al. [31], and has been applied to compute the response of both invariant [17] and periodic [32] surface tracks. The application of this methodology to the response of periodic tunnel-soil systems due to moving loads will also be discussed in a future publication.

Problem outline
The analysis of the three-dimensional dynamic soil-tunnel interaction problem is restricted to straight tunnels that are periodic with a period L in the longitudinal direction e y along the tunnel axis and embedded in a horizontally layered soil medium. The resulting problem is periodic and can be restricted to periodic fields of the second kind defined on a reference cellÕ (Fig. 1) by application of Floquet's theorem. The boundary qÕ of this domain is decomposed into the free surfaceG s and the boundaries S 0 and S L on which periodic conditions are imposed. The generic cellÕ is decomposed into two subdomains: the soilÕ s and the tunnelÕ t . The interface between these subdomains is denoted byS ts . The boundaryG ss is the free surface of the soil, while a surface forcef t is applied onG ts (Fig. 1).
The position vector x of any point in the problem domain O is decomposed as x ¼x þ nLe y , wherex is the position vector in the reference cellÕ and n is the cell number. The Floquet transformationf ðx; kÞ of a nonperiodic function f ðxÞ ¼ f ðx þ nLe y Þ defined on a three-dimensional domain O, that is periodic in the direction e y with period L, transforms the distance nL between the nth cell and the reference cellÕ to the wavenumber k and is defined as [20] The complex-valued functionf ðx; kÞ is a periodic function of k with a period 2p=L (or a periodic function of the first kind with respect to k), as the following periodicity condition holds: This complex functionf ðx; kÞ is also periodic of the second kind with a period L in the direction e y , since the following condition holds for allx:f ðx þ Le y ; kÞ ¼ expðÀikLÞf ðx; kÞ.
The function f ðx þ nLe y Þ can be reconstructed for any x ¼x þ nLe y using the inverse Floquet transform:

Navier equations
Using the Floquet transformation, all displacement and traction fields uðx; oÞ and tðx; oÞ defined on the periodic domain O are transformed to the fieldsũðx; k; oÞ andtðx; k; oÞ defined on the generic cellÕ. Subscripts s and t will be used in the following to denote displacements, tractions and other field variables related to the soil and tunnel subdomain, respectively.
The Navier equations in the soil domainÕ s and the boundary conditions onG ss are written as follows for every frequency o 2 R and wavenumber k 2 À p=L; þp=L½: where Eq. (5) expresses the conservation of momentum (with r s the density of the soil material) of the soil in the reference cell, Eq. (6) represents zero stresses on the partG ss of the free surface and Eq. (7) imposes the periodicity condition of the second kind on the soil displacement vectorũ s ðxÞ. This displacement field should also satisfy Sommerfeld's radiation conditions. The following Navier equations and boundary conditions hold in the tunnel domainÕ t with boundaryG ts : where r t is the density of the tunnel material andf t is the force per unit area applied on the boundaryG ts . Similar remarks can be made as for the Navier equations and boundary equations for the soil subdomain. Continuity of displacements and equilibrium of stresses must hold on the tunnel-soil interfaceS ts : t t ðũ t Þ þt s ðũ s Þ ¼ 0 onS ts .

Weak variational formulation
The Navier equations (9) onÕ t and the boundary equations on the boundaryG ts are multiplied with the complex conjugateṽ t ðx; kÞ of any virtual fieldṽ t ðx; kÞ and integrated over the problem domainÕ t and the boundaryG ts . Integration by parts on the volume integral containing the divergence of the stress tensor allows to derive the following weak variational form for the tunnel [18,20,22]: ZÕ whereẽðṽ t Þ :rðũ t Þ denotes the double contraction of the complex conjugate of the virtual strain tensorẽðṽ t Þ and the stress tensorrðũ t Þ, which is equivalent to the trace of the tensor productẽðṽ t Þ Tr ðũ t Þ. In Eq. (13), the boundary qÕ t consists of the tunnel-soil interfaceS ts and the two boundaries S 0 and S L at the two edges y ¼ AEL=2 of the generic cell. As the actual and the virtual displacement fields are periodic of the second kind, the contribution of the integral on the sum of the boundaries S 0 and S L vanishes [18]. Accounting for the stress equilibrium (12) along the tunnel-soil interfaceS ts , the weak variational equation becomes ZÕ whereũ sc ðũ t Þ denotes the wave field that is scattered by the tunnel into the soil and that obeys the displacement continuity relation (11) on the tunnel-soil interfaceS ts .

Coupled periodic FE-BE formulation
As the tunnel is bounded, the displacement fieldũ t ðx; k; oÞ can be decomposed as the following linear combination of functionsw m ðx; kÞ that are periodic of the second kind: where the modal coordinates a m ðk; oÞ are collected in a vector a t . The kinematic basis for the tunnel is determined such that the modesw m ðx; kÞ are periodic of the second kind: The modesw m ðx; kÞ are constructed asw where the vectorsw 0 m ðxÞ represent the eigenmodes of the reference cell, which are enforced to be periodic of the first kind (for The number M of modesw m ðx; kÞ in the expansion (15) of the tunnel displacement vector must be sufficiently high so as to incorporate the contribution of the relevant lowest modes of the tunnel, which depends on the frequency content and spatial distribution of the loading on the tunnel invert. This study is the subject of a convergence analysis. The soil displacementsũ s ðx; k; oÞ can be written as the superposition of waves that are radiated by the tunnel into the soil:ũ s ðx; k; oÞ ¼ũ sc ðũ t Þðx; k; oÞ ¼ The numerical solution of the dynamic tunnel-soil interaction problem is obtained using the classical domain decomposition approach based on the finite element method for the structure and the BE method for the soil.
The displacements in the structureũ t ðx; k; oÞ are interpolated as where N t is a matrix of global finite element shape functions (of dimension 3 Â N in a three-dimensional continuum formulation with N degrees of freedom) andũ t is the N-dimensional vector with nodal degrees of freedom. Global shape functions are preferred to local element-based shape functions so as to obtain a concise mathematical derivation. In practice, however, global finite element matrices are obtained by the assembly of element matrices, based on locally defined shape functions. Employing the same interpolation for the virtual displacementsṽ t ðx; k; oÞ, the following system of equations is finally obtained: where K t ðkÞ is the projection of the stiffness matrix K FE t of the tunnel reference cell on the tunnel modes (which are periodic of the second kind): In a three-dimensional continuum formulation, L is a 3 Â 6 matrix with differential operators that relates the strain vector to the displacement vector, while D is a 6 Â 6 matrix with constitutive coefficients that relates the stress vector to the strain vector. The finite element matrix K FE t of the reference cell is an N-dimensional matrix, whereas the dimension of the stiffness matrix K t ðkÞ is equal to the number of modes M in the decomposition (15).
The matrix M t ðkÞ is the projection of the mass matrix M FE t of the tunnel reference cell on the tunnel modes: F t ðk; oÞ is the generalized force vector applied on the tunnel invert: and K s ðk; oÞ is the dynamic stiffness matrix of the soil: The stressest s ðũ sc ðN tWt ÞÞ on the tunnel-soil interface are calculated with a periodic boundary element formulation with Green-Floquet functions or periodic Green functions defined on the periodic structure with period L along the tunnel. The reader is referred to complementary literature for more details on the periodic boundary element formulation [18][19][20]22].

Craig-Bampton substructuring method
In order to analyse the performance of different track structures with respect to vibrations generated in the free field and in adjacent buildings, it is advantageous to differentiate between the degrees of freedom of the tunnel invert and the track. The motivation for this approach is to avoid the recomputation of the soil impedance, which is only governed by the kinematics of the tunnel invert, when alternative track configurations are considered in the tunnel. Therefore, the displacement vectorũ t ðx; k; oÞ of the tunnel is discretized as follows, using a Craig-Bampton substructuring method [30]:ũ where the subscripts t r and t t refer to the track and the tunnel invert, respectively. The modes in Eq. (26) are periodic of the second kind and are constructed as follows from periodic modes of the first kind: where the diagonal matrices K t r t r and K t t t t are constructed according to Eq. (17). The modesW 0 t r are the eigenmodes of the track clamped at the tunnel invert. The modesW 0 t t are the modes of the free tunnel without track. The displacementsW s0 t r are the quasi-static transmission of the tunnel modes into the track, computed as where K FE t r t r and K FE t r t t are block submatrices of the finite element stiffness matrix K FE of the reference cell of the tunnel. The Craig-Bampton decomposition is frequently used for the seismic analysis of structures on flexible foundations [33]. It is equivalent to a classical modal decomposition (based on the modes of the complete track-tunnel structure) if all modes are included in both formulations. If a reduced modal basis is used, a convergence analysis is needed to determine the number of tunnel and track modes that need to be included in the Craig-Bampton decomposition (26), analogously to a convergence analysis for classical modal analysis.
Introducing the decomposition (26) in Eq. (14) results iñ K s ðk; oÞ is the dynamic stiffness matrix of the soil: F t r ðk; oÞ is the generalized force vector applied on the track: The soil impedance in Eq. (30) only depends on the tunnel modes and does not change when a calculation needs to be made for an alternative track structure in the tunnel, which contributes to the versatility of the proposed approach.

Wave propagation in the soil
When the displacementsũ t ðx; k; oÞ and the stressest t ðx; k; oÞ on the tunnel-soil interface are known, the incident wave fieldũ inc ðñ; k; oÞ is obtained by application of the dynamic representation theorem in the unbounded soil domain corresponding to the reference cell:  [18,20]. The incident wave field in the soil is obtained by evaluating the inverse Floquet transform (4).

Site characteristics
The tunnel on the Bakerloo line of London Underground in Regent's Park is a deep bored tunnel with a cast iron lining and a single track, embedded in London clay at a depth of about 28 m. The tunnel has an internal radius of 1.83 m and a wall thickness of 0.022 m. There are six longitudinal stiffeners and one circumferential stiffener at an interval of 0.508 m, resulting in a periodic structure (Fig. 2). The track is a non-ballasted concrete slab track with Bull head rail supported on hard Jarrah wooden sleepers via cast iron chairs. The sleeper distance is 0.9 m. Both ends of a sleeper are concreted into the invert and the space between the sleepers is filled with shingle. Resilience is mainly provided by the timber sleepers, as the rails are not supported by rail pads.
Geological maps show that the thickness of the London clay layer at the site is 40 m. GeoDelft has performed cone penetration tests (CPT) upto a depth of 21 m [34]. A shallow top layer with a thickness of 4-6 m has inclusions of sand and gravel and varying cone resistance. Bender element tests on undisturbed samples at several confining pressures result in an average shear wave velocity of 124 m/s and a longitudinal wave velocity of 1604 m/s [34]. A material damping ratio of 0.042 in the top layer and 0.039 in the second layer has been determined with free torsion pendulum tests [34]. A seismic cone penetration test (SCPT) upto a depth of 21 m confirmed the presence of a shallow stiff layer with a thickness of 4-6 m and a shear wave velocity C s ¼ 325 m/s on top of a homogeneous half-space with C s ¼ 220 m/s [34]. Spectral analysis of surface wave (SASW) tests confirmed the latter and revealed the presence of a homogeneous clay layer with a shear wave velocity between 200 and 260 m/s [35].
In the numerical calculations, the tunnel is assumed to be embedded in a layered soil consisting of a single shallow layer with a thickness of 5 m on top of a homogeneous half-space consisting of clay. The top layer has

Kinematics of the tunnel
The tunnel in London is periodic, as its lining is made of cast iron segments with circumferential stiffeners. The length L of the cell is equal to the segment length of 0.508 m. The sleeper distance is equal to 0.9 m and introduces a second periodicity, which is not accounted for as the track is not included in the model.
A FE model of the reference cell is made using the structural dynamics toolbox (Sdt) in Matlab. The exterior shell is modelled with four-node quadrilateral shell elements based on a Q4WT formulation for the membrane behaviour and the MITC4 formulation to represent the plate bending behaviour. In the latter, a mixed interpolation of the transverse displacement, section rotations and transverse shear strains is used to avoid shear locking when applied to thin shells [36]. The longitudinal and circumferential stiffeners are modelled with standard two-node Euler-Bernoulli beam elements, that are rigidly connected to the exterior shell using stiff eccentric connections. Eight-node brick elements are used to model the concrete on the tunnel invert.
The kinematic basis for the tunnel consists of the modesw m ðx; kÞ that are derived from the eigenmodes w 0 m ðxÞ of the tunnel cell with free boundary conditions onS ts and periodicity conditions at both ends S 0 and S L . Due to these constraints and the symmetry of the cell, displacements in the y-direction are decoupled from displacements in the x-and z-direction and only four rigid body modes are found. Fig. 3 shows the first two inplane and the first out-of-plane flexible modes w 0 m ðxÞ of the reference cell. In the following, the response of the tunnel-soil system due to an eccentric harmonic vertical load on the tunnel invert is computed in the frequency range upto 80 Hz. Convergence analysis has shown that minimum 20 tunnel modes of the reference cell need to be included in the analysis. Fig. 4a shows the modulus of the sixth diagonal element Z t ð6; 6Þ of the tunnel impedance matrix Z t ðk; oÞ ¼ K t ðkÞ À o 2 M t ðkÞ as a function of the frequency o and the slowness p ¼ k=o, which is defined as the ratio of the wavenumber k to the circular frequency o. The element Z t ð6; 6Þ corresponds to the sixth mode of the reference cell (Fig. 3b) involving in-plane deformation of the reference cell, equivalent to a mode of azimuthal order m ¼ 2 in an isotropic shell (i.e. with 2 wavelengths along the circumferential direction). Fig. 4a clearly illustrates that the sixth eigenfrequency of the tunnel reference cell corresponds to the cut-on frequency of a dispersive mode that propagates in the longitudinal direction of the tunnel. Fig. 4b shows the variation of the element Z t ð6; 6Þ of the tunnel impedance matrix with the frequency o for a zero value of the slowness p ¼ 0 s/m, corresponding to the two-dimensional plane strain case. At o ¼ 0 rad/s, the tunnel impedance is equal to the square of the sixth eigenfrequency of the reference cell. Due to the inertial term, it decreases with the square of the frequency o and becomes zero at the sixth resonance frequency of the reference cell, which corresponds to the cut-on frequency of a dispersive mode for which the phase velocity is infinite at zero slowness or wavenumber in the undamped case. For non-zero values of p, the tunnel impedance increases with the square of the slowness p. Fig. 4b also shows the variation of the element Z t ð6; 6Þ of the tunnel impedance matrix with the frequency o for two non-zero values of the slowness p ¼ 0:3 Â 10 À3 and 0:6 Â 10 À3 s/m.

Impedance of the tunnel
The response of the tunnel is determined by the contribution of different tunnel modes. Fig. 5 shows, on a logarithmic scale, the modulus of the vertical displacement of the tunnel invert under the load as a function of the circular frequency o and the slowness p for the free tunnel that is not supported by the soil. In all, 20 modes of the reference cell are included in the analysis. The response is governed by dispersive tunnel modes, that start at the eigenfrequencies of the tunnel reference cell; the contribution of the modes 5 till 8 with a cuton frequency at 20.94, 24.32, 62.74 and 68.04 Hz can clearly be observed. part increases for increasing values of p. Due to the sign convention used in the kernels of the Fourier and Floquet transformations, the imaginary part is negative. High (negative) values of the imaginary part reflect dissipation of energy due to radiation damping in the soil. The low imaginary part of the soil impedance for high values of the slowness p (or decreasing phase velocity along the longitudinal tunnel axis) reflects the absence of radiation damping in the soil.

Impedance of the soil
Bostro¨m and Burden [37] have demonstrated that the response of an infinite elastic medium with a cylindrical cavity is governed by dispersive surface waves modes, corresponding to integer values of the azimuthal wavenumber m; these modes are confined to the vicinity of the cavity and have a cut-on frequency depending on the azimuthal wavenumber m, the characteristics of the soil domain and the radius of the tunnel. They propagate with a phase velocity that varies between the shear wave velocity (at the cut-on frequency) and the Rayleigh wave velocity (at limiting high frequency). As the Bakerloo line tunnel is embedded in a deep homogeneous layer of clay, the behaviour of an element of the soil impedance matrix is very similar to the response of a cylindrical cavity in an elastic medium, when a particular mode of this cavity is considered. The difference, however, is that the kinematics along the tunnel-soil interface are governed by the modes of the (orthotropic) reference cell of the tunnel that do not perfectly match the azimuthal modes of a cylindrical cavity or an isotropic shell. Based on the analytical solution of Bostro¨m and Burden, the cut-on frequency of the mode corresponding to the azimuthal mode number m ¼ 2 of a cylindrical cavity with a radius of 1.953 m, embedded in a homogeneous and undamped medium with the characteristics of the London clay, is equal to 65 Hz.

Response due to harmonic loading
The response of the Bakerloo line tunnel is considered for a harmonic loading on the tunnel invert in the frequency range between 1 and 80 Hz. Fig. 7a shows the modulus of the vertical displacement of the tunnel invert as a function of the circular frequency o and the slowness p for the case of a free tunnel that is not supported by the soil. These results correspond to those presented in Fig. 5, but now a wider range of the slowness p is considered to enable the comparison with results that follow in this subsection. The response is governed by the dispersive character of the tunnel modes at low values of the slowness, as explained before. Fig. 7b considers the response of the layered soil medium with a cylindrical cavity, only accounting for the impedance of the soil and neglecting the dynamic impedance of the tunnel. The kinematics of the free surface along the cavity (i.e. the interface between the soil and the tunnel in the dynamic soil-tunnel interaction problem) is described by 20 modes of the tunnel's reference cell. The response is dominated by dispersive modes that propagate with a phase velocity that varies between the shear wave velocity C s ¼ 220 m/s of the soil around the tunnel (at the cut-on frequency of each mode) and the Rayleigh wave velocity C r ¼ 209:75 m/s at limiting high frequency. The corresponding values of the slowness vary between 1=C s ¼ 4:545 Â 10 À3 s/m and 1=C r ¼ 4:768 Â 10 À3 s/m. Fig. 7c finally shows the response of the coupled soil-tunnel system in the frequency-slowness domain. Furthermore, Fig. 8 shows sections of the vertical response of the tunnel invert as a function of the slowness p for four values of the frequency (21, 41, 61 and 81 Hz). Results are shown for the three cases considered earlier: the free tunnel (Fig. 7a), the soil medium with a cylindrical cavity (Fig. 7b) and the coupled soil-tunnel system (Fig. 7c). Fig. 8b, for example, clearly illustrates that the response of the free tunnel at 41 Hz is dominated by two dispersive bending modes of the tunnel, while the response of the soil medium with a cylindrical cavity is maximum for values of the slowness between the reciprocals of the shear wave velocity and Rayleigh wave velocity of the soil around the tunnel. For low values of the slowness, the dynamic stiffness of the soil is larger than the dynamic stiffness of the tunnel, and the peaks in the response of the free tunnel are no longer apparent in the response of the coupled tunnel-soil system. At high values of the slowness, the dynamic stiffness of the tunnel is very high compared to the soil's stiffness, and the response of the coupled soil-tunnel system converges to the response of the free tunnel, while the effect of dispersive modes along the tunnel-soil interface is no longer observable.
The response in the spatial domain is obtained after evaluation of the inverse Floquet transformation. Figs. 9 and 10 show the response of the Bakerloo line tunnel for a harmonic excitation on the tunnel invert at 10 and 40 Hz, respectively. The tunnel's cross section is relatively stiff but the tunnel diameter is small, contributing to the flexibility of the tunnel in the longitudinal direction. As the Bakerloo line tunnel is embedded at a depth of 28 m in a homogeneous layer of London clay, the response is not affected by resonance in the soil. The wave pattern is similar to that caused by a point load in the soil. The propagation of waves into the soil is spherical but not symmetric, due to the eccentricity of the load on the tunnel invert with respect to the centre of the tunnel.

Floating slab track
Next, the vibration isolation efficiency of a floating slab track with a resonance frequency of 11.68 Hz installed in the Bakerloo line tunnel is considered, using the Craig-Bampton substructuring technique. This example should be regarded as illustrative rather than practical, as the small diameter of the tunnel does not allow to install a floating slab track. Figs. 11 and 12 show the displacements of the tunnel and the soil for a harmonic excitation at 10 and 40 Hz after installation of a floating slab in the tunnel. It should be noted that the grey scales are different from that in Figs. 9 and 10, as the displacements are generally lower. The floating slab track distributes vibrations along the tunnel that are radiated into the soil in directions parallel to the tunnel. Fig. 13 compares the displacements at the tunnel apex and in the free field above the tunnel for the cases without and with floating slab track; the isolation effect above 20 Hz can clearly be observed.

Site characteristics
The metro tunnel on the line RER B of RATP at Cite´Universitaire in Paris is a masonry cut-and-cover tunnel with two tracks at a shallow depth of about 9.3 m below the surface and a width of 11.9 m (Fig. 14). The  slab thickness is 0.6 m at the top and 0.4 m at the bottom, while the wall thickness is 1.5 m. The masonry is assumed to have a Young's modulus E t ¼ 14; 000 MPa, a Poisson's ratio n t ¼ 0:15, a density r t ¼ 2400 kg/m 3 and a hysteresis material damping ratio b t ¼ 0:02. The track is a classical ballasted track with UIC 60 rails supported every 0.60 m by grooved rubber pads on monobloc concrete sleepers.    [38]. A Poisson's ratio n s ¼ 0:4, a density r s ¼ 1700 kg/m 3 and a material damping ratio b s ¼ 0:05 is assumed in all layers. More recent down-hole measurements have confirmed the variation of the shear wave velocity with depth [39], and have demonstrated that the longitudinal wave velocity increases from 700 m/s at a depth of 3 m to 900 m/s at a depth of 10 m, resulting in a value of Poisson's ratio that decreases from 0.43 to 0.33 (and showing no evidence of the presence of ground water). A higher average density of 1902 kg/m 3 was determined from soil samples, but has not been used in the calculations.

Kinematics of the tunnel
The cut-and-cover masonry tunnel is invariant in the y-direction. The periodicity introduced by the discrete support of the rails is neglected, since the track is not included in the model. The tunnel is invariant in the ydirection, so that the period L can arbitrarily be chosen as 0.3 m; the influence of the assumed period L of the tunnel is presently investigated by a comparison with the results of a 2.5D calculation based on the semianalytical insertion loss model [9][10][11]. The FE model of the tunnel consists of eight-node isoparametric brick elements. As no incompatible modes are included in the formulation of these elements, three elements are used over the height of the walls and ceilings, and two elements for the bottom slab, in order to represent the bending behaviour. The size of the finite elements is governed by the boundary element mesh along the tunnel-soil interface, so that minimum eight elements are used per wavelength at a maximum frequency of 80 Hz; the shear wave velocity C s ¼ 220 m/s in the soil layer around the apex of the tunnel has been used to compute this wavelength. Fig. 15 shows the first two in-plane and the first out-of-plane flexible modes w 0 m ðxÞ of the reference cell. The response of the tunnel-soil system due to a harmonic load on the tunnel invert (in the centre of the left track in Fig. 14) will be computed in the frequency range upto 80 Hz. Convergence analysis has shown that 30 tunnel modes need to be included in the analysis [22]. Fig. 16a shows the modulus of the sixth diagonal element Z t ð6; 6Þ of the tunnel impedance matrix Z t ðk; oÞ ¼ K t ðkÞ À o 2 M t ðkÞ as a function of the frequency and the slowness p. The element Z t ð6; 6Þ corresponds to the  sixth mode of the reference cell (Fig. 15b) involving in-plane bending of the slab of the tunnel. Fig. 16a clearly illustrates that the sixth eigenfrequency of the tunnel reference cell corresponds to the cut-on frequency of a dispersive mode that propagates in the longitudinal direction of the tunnel. The phase velocity is infinite at the cut-on frequency and tends to a constant value of 1=C s ¼ 0:628 Â 10 À3 s/m for limiting high values of the frequency, where C s is the shear wave velocity in the masonry. Fig. 16b shows the variation of the element Z t ð6; 6Þ of the tunnel impedance matrix with the frequency o for three values of the slowness p. Similar figures can be made for other modes of the tunnel's reference cell, but are not shown here. Fig. 17 shows, on a logarithmic scale, the modulus of the vertical displacement of the apex of the tunnel as a function of o and p for the free tunnel that is not supported by the soil, including the contribution of 30 modes of the tunnel reference cell. The response is governed by dispersive tunnel modes, of which the slowness tends to a value 1=C s ¼ 0:628 Â 10 À3 s/m for the in-plane modes and 1=C p ¼ 0:403 Â 10 À3 s/m for the out-of-plane modes, with C s and C p the shear and longitudinal wave velocities in the masonry. Fig. 18 shows the real and imaginary parts of the element K s ð6; 6Þ of the impedance matrix K s ðk; oÞ of the soil, corresponding to the sixth mode of the RER B tunnel. The real part increases for increasing values of p. High negative values of the imaginary part reflect dissipation of energy due to radiation damping or wave propagation in the soil, while the low imaginary part for high values of the slowness p (or low values of the phase velocity along the tunnel axis) represents the absence of this radiation damping. As the RER B tunnel is shallow and embedded in a layered soil medium, the behaviour of the soil impedance as a function of o and p is more complicated than for the Bakerloo line tunnel, which was a deep tunnel embedded in a homogeneous soil layer. The soil impedance is now not only affected by the heterogeneous soil around the tunnel, but also by the proximity of the free surface.    19 shows the transfer functions (vertical displacements) at some points along the free surface in the cross-section where the load is applied, while Fig. 20 shows the transfer functions at points on the tunnel invert, the tunnel apex and the free surface. The peak at 14 Hz is due to resonance of the soil above the tunnel; Fig. 20 confirms that the displacements of the tunnel apex and the soil above the tunnel are equal. Fig. 21 shows the displacements of the tunnel and the soil due to a harmonic loading on the tunnel invert at 14 Hz. The soil above the tunnel can be considered as a mass, which is moving in phase with the tunnel, that behaves as a spring. The peak at 14 Hz corresponds to the eigenfrequency of this equivalent single degree of freedom system (SDOF) system, while damping is due to the radiation of waves away from the tunnel. Fig. 22 shows similar results at a frequency of 80 Hz. On the free surface above the tunnel, higher phase velocities (about 1920 m/s) are observed along the tunnel axis than in the direction perpendicular to the tunnel (about 650 m/s), due to the bending stiffness of the tunnel, resulting in an elliptical wave front. At a frequency of 80 Hz, the displacement of the ground surface and the tunnel apex are almost in phase, as this frequency almost corresponds to the third vertical mode of the soil above the tunnel, as has been demonstrated by analysis of the dispersive characteristics of the waves in the soil above the tunnel [40].

Conclusion
A periodic coupled finite element-boundary element formulation is used to study the dynamic interaction between a tunnel and a layered soil due to a harmonic excitation on the tunnel invert. Two cases have been considered: a deep-bored tunnel on the Bakerloo line of London Underground and a shallow cut-and-cover masonry tunnel on the RER B line of RATP in Paris.
The Bakerloo line tunnel has a complex structure as the cast iron lining is stiffened with longitudinal and circumferential beams, resulting in a periodic structure. The periodic finite element formulation offers full modelling flexibility on a reference cell of the tunnel, which is an important advantage with respect to full three-dimensional or invariant tunnel models, which can be considered as the two extremes of modelling complexity. The soil impedance is formulated with a boundary element formulation on the soil-tunnel interface in the reference cell, allowing to model the soil as a horizontally layered half-space. This is still quite demanding from the computational point of view, however, depending on the size of the boundary element mesh and the evaluation of the Green-Floquet functions. As the Bakerloo line tunnel is an embedded tunnel, the interpretation of the soil's impedance benefits from the analytical study of the response of a homogeneous elastic medium with a cylindrical cavity. Analysis in the frequency-slowness domain demonstrates that, for the case under consideration, the soil dominates the response at low values of the slowness, while the tunnel dominates the response at high values of the slowness. The resulting response of the coupled soil-tunnel system is not dominated by dispersive modes in the tunnel and the soil. This affects the results in the frequency-spatial domain, showing spherical wave patterns in the soil, as would also result from a point load applied to the soil. These preliminary conclusions depend of course on the relative stiffness of the tunnel and the soil, as well as other factors such as the amount of material and radiation damping in the soil. The general  expectation is that, at low frequencies, the influence of in-plane deformations of the tunnel is low, whereas more complicated mode shapes are excited at higher frequencies, affecting the wave field radiated into the soil. However, this could not be clearly observed on the results for the Bakerloo line tunnel.
The case of the RER B tunnel in Paris is more complicated, as a shallow cut-and-cover tunnel with a noncircular cross-section is embedded in a heterogeneous layered soil. The only simplification here is that the tunnel is invariant in the longitudinal direction, allowing more freedom in the choice of a (shorter) period (which has not been studied in this paper and constitutes a subject of further research). It has been demonstrated that the response of this shallow tunnel is importantly influenced by the proximity of the free surface, giving rise to resonance phenomena, as well as to elliptical wave fronts radiated from the tunnel that originate from different phase velocities in the longitudinal and perpendicular direction.
It has also been demonstrated with an example how the Craig-Bampton substructuring technique allows to efficiently incorporate a track in the tunnel and to analyse the vibration isolation efficiency of a floating slab track in the tunnel, without the need to recompute the impedance of the soil. This will substantially facilitate future parametric studies on the efficiency of vibration isolation measures in the tunnel.