Modelling tunnel behaviour under seismic actions: an integrated approach

: This paper intends to describe the integration of physical and numerical modelling, focusing on tunnels under seismic actions. It shows how numerical calculations can be used in association with centrifuge testing to model different aspects of tunnel behaviour during earthquakes. The scope of the paper has been limited to a few aspects, mainly concerning the change of internal forces in the tunnel lining during shaking and the effect of soil liquefaction. The interaction between a tunnel and a building in a soil layer undergoing liquefaction has also been taken into account.

focusing on broad and impacting research streams only. Less appealing problems, that receive lower attention, might be excluded from the benefit of such a combined approach. Repositories of the experimental data that are produced by different facilities for different purposes all over the world, play in this case a fundamental role. This paper intends to describe the integration of physical and numerical modelling from the point of view of numerical modellers. Focusing on the dynamic behaviour of tunnels, and in particular on the internal forces in the tunnel lining, the use of the centrifuge results to calibrate a numerical model and extend the scope of application is shown in section 2 and section 3. In the former the results of centrifuge testing are boosted by including numerically the effect of a construction process for tunnelling and a more complex structural behaviour of the lining. In section 3 the back-analysis of the centrifuge test in dry sand is used in association with the results of cyclic simple shear testing in undrained conditions, for the calibration of a constitutive model suitable for modelling pore-water pressure build-up in undrained conditions. The effect on the lining of the excess pore-pressure arising during shaking is analysed. Finally, in section 4 the process is reversed from numerical analysis to centrifuge modelling. The calibration carried out in the previous sections is used to perform a preliminary analysis of tunnel-building interaction in liquefiable soil, in order to design a series of centrifuge tests.

Background
Internal forces in the tunnel lining change during earthquakes. They can be calculated following several approaches (Hashash et al., 2001;Pitilakis & Tsinidis, 2014). Pseudo-static or uncoupled dynamic analyses are usually carried out in routine design. Full dynamic analysis, that is including dynamic soil-structure interaction, must be performed however, if the influence of the existing stress state around the tunnel has to be considered. Moreover, the latter may include the irreversible behaviour of soil that is likely to produce permanent ground deformation during shaking. Since the tunnel construction process may affect the static conditions before shaking, numerical analyses can include this aspect. Compared to plane strain, three-dimensional models permit the construction phases to be simulated in a more accurate fashion, including geometrical details of the lining that may affect its structural behaviour (for instance the segmental layout of precast lining). The effect of seismic waves propagating in any direction can be also analysed in a three-dimensional numerical model.
On the other hand, direct measurements of the effect of the complex interaction between a tunnel lining model and the surrounding soil during ground shaking can be achieved in centrifuge tests. This enables a large amount of experimental data to be collected and used for validation of numerical analyses. As part of a research within the ReLUIS project funded by the Italian Civil Protection Department, a series of centrifuge tests were carried out at the Schofield Centre of the University of Cambridge on circular tunnel in dry sand, undergoing dynamic excitation (Lanzano et al., 2012). Internal forces (bending moments and hoop forces) in the tunnel lining were measured during shakings. Hence, experimental evidence was gained on a problem that had been previously explored via analytical solutions (mainly based on the elastic theory) and numerical modelling only.
Such tests, which are briefly recalled in the next section 2.2, were later used as an experimental benchmark for numerical modelling, aimed at extending the scope of the study. In fact, three-dimensional finite element analyses were performed that take into account the non-linear and irreversible soil behaviour. The tunnel excavation process, that is neglected in the centrifuge tests, was modelled to achieve a realistic state of stress effect before shaking. Moreover, the segmental layout of a precast tunnel lining was modelled, although with a few simplifying assumptions (Fabozzi, 2017).

Experimental benchmark
The experimental benchmark used for validating the numerical model of a rather shallow tunnel (C/D=2) in dense sand is the centrifuge model T3 (Figure 1), described in details by Lanzano et al. (2012).  Lanzano et al., 2012).
In the model (Fig. 1a), an aluminium tube (diameter D = 75 mm, thickness t = 0.5 mm, cover C = 150 mm) representing a circular tunnel is embedded in a layer of dry Leighton Buzzard sand (fraction E), pluviated in the container at a relative density of 75% ( Figure 1). The tube is instrumented with strain gauges in four positions along its transverse section (indicated as NE, NW, SW and SE in Fig.1a). After spin up at 80g, a series of pseudoharmonic signals of increasing amplitude and frequency was applied at the base of the model. The time histories of bending moment, M, and hoop forces, N measured during four subsequent dynamic events are shown in Fig. 1b at the model scale (Lanzano et al., 2012). It is worth noticing that permanent increments of internal forces arose in the tunnel lining after each events. These seem well correlated to the progressive densification of the sand layer that was observed in the experiments.

Numerical modelling
Numerical analyses were performed at prototype scale, using a scaling factor N=80. Hence, the corresponding tunnel diameter is assumed 6 m, the tunnel axis depth is 15 m and the lining thickness is comparable to that of a concrete lining about 0.06 m thick. The numerical model has been implemented in the finite element code Plaxis 3D (Brinkgreve et al., 2016). The mesh is shown in Figure 2.
While the height of the model is 23.2 m, that is 80 times the relevant size at model scale, its width is larger than that and equal to 200 m, to minimise the influence of lateral boundaries. A reference section at the midspan of the tunnel was assumed to be compared to the experimental results. Hence, in order to guarantee plane strain conditions in the reference section, the size of the model along the axis of the tunnel was assumed as long as 150 m. The vertical sides of the mesh were fixed in the horizontal direction in static condition; viscous dashpots were applied during shaking.
The time history of acceleration recorded by the accelerometer ACC13 at the base of the centrifuge model (see Figure 1) was scaled up to prototype scale and band-pass filtered  in order to reduce the its high-frequency content. This signal (with nominal frequency 0.375 Hz and nominal amplitude 0.05 g) was applied as dynamic input at the base of the model. The lining is an elastic plate (EA=2.810 6 kN/m; EI=3.710 2 kNm 2 /m) with a very smooth interface (the interface factor was assumed as Rint = 0.05).
The sand has been modelled using the Hardening Soil with small strain overlay constitutive model (Schanz et al., 1999;Benz et al., 2009), with the parameters shown in Table 1, derived by Lanzano et al. (2016). Table 1. HS-small model parameters (Lanzano et al., 2016).
This elastic-plastic with isotropic hardening soil model is able to reproduce the decay of shear stiffness with strain level from very small strain and the increase of hysteretic damping. The initial damping ratio at very small strain was modelled through a Rayleigh formulation (αR = 0.0668; βR = 0.704 10 -3 ). Figure 3a compares the time history of acceleration measured in the test by ACC9 with the corresponding computed results. In Figure 3b the corresponding response spectra at 5% of damping are shown. The dynamic response computed for the soil layer is close to the measurements, although there is evidence of a slight overamplification of the signal at high frequencies, as observed also by Amorosi et al. (2014) in similar analyses.  Once validated against centrifuge results, the same 3D model was used to analyse the behaviour in the same sand of a different tunnel lining. This is a reinforced concrete lining with thickness t = 0.3 m (EA = 10.5E 6 kN/m; EI = 78.75E 3 kNm 2 /m; Rint = 0.7) and diameter D = 6 m. A set of natural input signals was applied as time histories of acceleration at the base of the mesh. A few results of the study (Fabozzi, 2017) are presented in sections 2.4 and 2.5: the influence of the construction process on the seismic response of the tunnel is discussed in the former while the latter analyses the influence of the presence of joints in the segmental lining.

Pre-seismic conditions induced by tunnel construction
The influence of the construction process has been taken into account with reference to typical mechanized tunnelling with an earth pressure balance machine. Details of the procedure are described by Fabozzi & Bilotta (2016) and will not be discussed here. The seismic excitation was applied to the numerical model at the state of stress corresponding to the end of construction. Table 2 shows the main characteristics of the input signals applied as time history of acceleration at the base of the model. They are natural time histories of acceleration recorded on a rigid outcropping bedrock (soil type A according to EC8). Their mean response spectrum matches the Eurocode EC8-1 spectrum for ground type A (rock). As an example, Figure 4a shows one of the time histories of acceleration applied at the base of the model. It is the record of the Norcia earthquake in Central Italy on 30/10/2016 (M w =6.5). In Figure 4b the corresponding normalized Fourier spectrum is shown. In all the analyses that have been carried out, permanent changes of internal forces in the lining at the end of shaking were calculated. In some cases, they reach values as high as 30% of the maximum transient change during shaking. As an example, in Figure 5 a pair of time histories calculated for the input signal of the Norcia Earthquake (see Figure 4) are shown. They are the time histories of the increment of bending moment (a) and hoop force (b) calculated at the point NE of the reference central section of the tunnel lining. The experimental evidence obtained by Lanzano et al. (2012) and shown in Fig. 1b are therefore confirmed in part by numerical modelling on a different lining and for different characteristics of ground shaking: permanent changes of internal forces are calculated at the end of shaking, as observed in the experiments, although they do not exceed the transient changes calculated during the event. It should also be remarked that, in order to capture such an effect a suitable elastic-plastic constitutive model for soil must be adopted, as in this case.     6c) were calculated in the transverse reference section, both after simulation of the tunnel construction process (black lines) and for an ideal "wished-in-place" tunnel (grey lines). As one would expect, the stress change due to the excavation produces lower bending moments ( Fig. 6a) and normal forces (Fig. 6b) in the tunnel lining, than in a wished-in-place tunnel. Furthermore, the latter is almost not loaded in longitudinal direction (Fig. 6c).
It is worth noting that, although the maximum values of pre-shaking internal forces (continuous lines) are quite different, such differences reduce after shaking (dashed lines). This implicitly means that the calculated permanent changes of internal forces depend on the pre-seismic conditions: when the excavation process is modelled they are larger than in the case of the wished-in-place tunnel. The effect of the construction stages on the seismic behaviour of the tunnel lining is therefore evidenced by such numerical results.

Segmental layout of the tunnel lining
Mechanised tunnelling in soft ground is usually associated with the use of a pre-cast concrete segmental lining to withstand external loads from interaction with the surrounding soil. Due to such a segmental layout the structural demand of the lining under static conditions is usually lower, because its flexural and axial stiffness is lower compared to a continuous lining of the same thickness. The same numerical model that was described in section 2.4 was improved to introduce a segmental lining. The segments were modelled as elastic volumes of reinforced concrete with the same thickness as the continuous lining (EA = 10.5E 6 kN/m; EI = 78.75E 3 kNm 2 /m). Following Fabozzi (2017), the longitudinal joints between the segments were modelled as elastic-plastic elements (thickness = 0.30 m, width = 0.30 m): the values adopted for their mechanical parameters are shown in Table 3. Interface elements with the same behaviour were assumed to represent the transverse joints between rings. Figure 7 shows details of the structural model. The excavation stage was not modelled. Table 3. Model parameters for the lining (Fabozzi, 2017)

Continuos lining Segmental lining
A lower structural requirement for the segmental lining compared to the continuous lining is evident also at the end of shaking. In Figure 8, the distributions of bending moment (Fig. 8a) and hoop force (Fig. 8b) in the transverse section at the end of shaking are shown.
The lower values of structural forces in the segments correspond to a larger deformability of the lining system at the joints, where relative displacements and rotation might be expected to occur. Figure 9 shows the time histories of relative rotation between segments during shaking, calculated in the joints located at 45°, 135°, 225°, and 315° about the horizontal tunnel axis. At the end of shaking permanent relative rotations remain between segments. The magnitude of such permanent rotations is sometimes rather close to the peak values calculated during shaking. This result indicates a possible weakness of the segmental lining at the joints, where the rubber gaskets that guarantee water-tightness of real linings might be dislocated at the end of an earthquake.
The results in terms of relative rotations for the whole set of input signals shown in Table 2 are plotted in Figure 10. In Figure 10a a linear trend can be observed for the logarithm of the calculated peak relative rotation between segments versus the value of the peak ground acceleration of the input signal (PGA). It is worth noting (Figure 10b) that as far as the peak relative rotation increases (with increasing PGA), the permanent relative rotation increases faster. The ratio between the permanent and the peak rotation increases from as low as 10% until almost one half for the stronger earthquakes. (a) (b) Figure 10. Peak relative rotation vs. peak ground acceleration (a) and permanent relative rotation vs. peak relative rotation (b), all input signals in Table 2.
This further highlights the influence of the non-linear behaviour of the surrounding soil on the value of permanent rotations experienced by the segmental lining at the end of shaking.

Remarks
The numerical analyses presented in this section were calibrated on a single benchmark centrifuge test and then extended to model more complex cases in terms of the geometry of the lining. The numerical model also allowed a straightforward application of natural input signals and a consideration of the influence of construction process on the seismic demand of the tunnel lining. It is worth noting that an earthquake can hit a tunnel several years after construction, hence different "preseismic" conditions can be considered. Moreover, in earthquake-prone regions, the same tunnel may be subjected to sequences of seismic events, with variable intensity and effects. Hence, the influence of the "initial state" should be considered in the assessment of tunnel vulnerability.
Moreover, the numerical results from the segmental layout may create some concerns for tunnel linings in highly permeable soils, where an excessive rotation of joints may produce loss of water-tightness. This aspect may deserve attention in design and, at the same time, requires further experimental and numerical investigation.

Background
Soil liquefaction may induce buoyancy of underground structures such as tanks, tunnels and pipelines. This is triggered when high excess water pressures develop, as those induced by strong motions. Several cases of uplift of underground tanks and pipelines have been observed in the past.
Although little evidence of liquefaction-induced damage to tunnels exists, physical modelling has shown that the high mobility of liquefied soil near surface would encourage floatation of very shallow or immersed tunnels. As a matter of fact, the uplift behaviour of underground structures caused by liquefaction has often been studied by physical models: for instance 1-g shaking table models of buried box structures, sewers and pipes and relevant possible mitigation measures (Koseki et al., 1997;Otsubo et al., 2014;Watanabe et al., 2016) or centrifuge models of tunnel of different shapes embedded in sand layers of different density, with several overburden and groundwater level (Yang et al. 2004;Chou et al, 2010;Chian and Madabushi, 2011;Chian & Madabhushi, 2012;Chian et al., 2014).
Experimental evidence has indicated that both the width of the underground structure and the depth of the liquefied layer have a large influence on the uplift displacement.
In general, physical modelling has been useful to collect important information and quantitative data on the phenomenon. In fact, although many numerical tools have been developed in the last decades to assess soil liquefaction, prediction of soil behaviour after liquefaction is still a challenging task. Hence, physical modelling has a further important role, that is to validate numerical models that can be used later for sensitivity analysis.
In this section, it is shown how starting from the back-analysis of the results of a centrifuge test on a model tunnel in dry sand undergoing shaking, the behaviour of the same tunnel in sand that has been saturated can be modelled numerically. The centrifuge model T4, described in detail by Lanzano et al. (2012), was used as an experimental benchmark. This centrifuge model has the same layout as model T3 in Figure 1, although the sand layer was looser (D r =40%).
The UBC3D-PLM constitutive model (Beaty & Byrne, 1998;Galavi et al., 2013) was used to represent the sand. It includes hardening plasticity and strain dependency of stiffness and damping. Hence it is able to capture the permanent deformation of the ground and changes in internal forces in the tunnel lining due to dynamic loading. Moreover, in undrained conditions it models the pore pressure build-up that may produce soil liquefaction and tunnel uplift.
The model is available in the 2D finite element code Plaxis (Brinkgreve et al., 2016) that has been used for the analyses. It was calibrated on the results of laboratory tests on the sand used in the centrifuge test along monotonic (Lanzano et al., 2016) and cyclic (Mele et al., 2018) stress paths.

Numerical analyses
A plane strain numerical model was defined in Plaxis 2D (Brinkgreve et al., 2016), at prototype scale. 'Tied degrees of freedom' between vertical sides were used as boundary conditions to simulate the laminar box behaviour during shaking. The nodes at the base of the finite element model were fixed in the vertical direction and a time history of acceleration was applied in the horizontal direction. The input signal applied at the base of the model is a pseudo-harmonic signal with nominal frequency 0.375 Hz and nominal amplitude 0.05 g at prototype scale. It was obtained after scaling up and filtering of the record of the base accelerometer in the centrifuge model ACC13 (see Figure 1).

Model calibration
The UBC3D-PLM is an elastoplastic constitutive model, which is a generalized formulation of the original UBCSAND model proposed for cyclic loading by Beaty & Byrne (1998). The model uses isotropic hardening and a simplified kinematic hardening rule for primary and secondary yield surfaces respectively, in order to take into account the effect of soil densification and transition to the liquefied state during undrained cyclic loading. Table 3. UBC3D-PLM model parameters sand The constitutive model is capable of modelling cyclic liquefaction for different stress paths (Galavi et al., 2013). Table 3 reports the input parameters used in the UBC3D-PLM model. The calibration of the model mechanical parameters was performed by Colamarino et al. (2017) using the results of laboratory tests on the sand used in the centrifuge test along monotonic (Lanzano et al., 2016) and cyclic (Mele et al., 2018) stress paths.

Response of the model with dry sand
The numerical results are compared to the centrifuge test results in terms of time history of acceleration and the relevant response spectrum, 1.6 m below the ground surface (position of ACC9 in Figure 1) at prototype scale ( Figure 11). For the sake of comparison, here and in the following figures, the experimental results of the centrifuge test are scaled up to prototype scale. Other relevant comparisons between recorded data and simulation results are reported in terms of vertical displacements and bending moment in Figure 12a and b, respectively. The main features of the experimental data are well-reproduced by the numerical predictions, both in terms of amplitude and frequency content, although an amplification larger than measured is calculated around 0.6 s, that is close to the natural period of the soil layer. Here, the agreement between the measured and the calculated amplitude achieved by using the UBC3D-PLM model is worse than by using HS-small model in similar conditions (see Figure 3). This might be in part due to the different amount of damping that the two constitutive models generate in stress-strain cycles. Figure 12a shows that the numerical model computes settlement at ground surface since the very beginning of the analysis, before 10 s, that is when the amplitude of the input signal is still negligible. In the same time a slight increase of bending moment is calculated (Fig. 12b). On one hand this confirms the influence of sand densification (hence plastic volumetric deformation) on the permanent change of internal forces in the lining; on the other it also shows that the numerical model tends to overpredict the plastic volumetric strain during shaking. As a consequence, the residual value of bending moment that is calculated at the end of shaking is even larger than the experimental value, although the corresponding transient cyclic changes are very similar (Fig. 12b).

Response of numerical model in saturated sand
The same input motion was applied at the bottom of the mesh modelling the soil as completely saturated. In this condition, significant excess pore pressure developed in the soil layer above the tunnel, although full liquefaction was not triggered, due to the low amplitude of the input signal. The excess pore pressure ratio, ru, defined as the ratio between the generated excess pore pressure and the initial effective vertical stress, did not exceed 0.77 ( Figure 13). Differences in the internal forces between dry and saturated conditions are shown in Figure 14. This figure shows that the hoop force increased at both control points (NE and SE) in saturated sand compared to dry sand (Figure 14a,b), while bending moment increased in the upper part (NE) of the tunnel cross section and decreased in the lower part (SE). This indicates that the pore pressure build-up, associated with changes in effective stresses, affects the distribution of internal forces in the tunnel lining.
In general, a larger change of hoop force is induced in the lining during ground shaking if soil liquefaction approaches. The effect on bending moment depends on the position along the lining. However, for such a lining (the very flexible one used in the experiment) the values of bending moments are very low.
In order to evaluate the preliminary remarks that emerge on the basis of the comparison in Figure 14, the tunnel lining was changed, as in section 2.3, to a thicker reinforced concrete lining (EA = 10.5E 6 kN/m; EI = 78.75E 3 kNm 2 /m).
Moreover, the numerical analyses were performed by applying at the base of the mesh the same signal 'EQ1' recorded in the centrifuge, two more signals obtained simply by scaling up 'EQ1' to twice ('2x') and three times ('3x') its amplitude, and additionally the six natural signals shown in Table 2.
An overview of the analyses is given in Table 4. The peak acceleration of the input signal (amax,b), peak acceleration calculated at the ground surface (amax,s), the maximum change of bending moment (ΔM) and hoop force (ΔN) in the tunnel lining at the end of shaking and the maximum uplift of the tunnel (uv,max) are shown in the table. The last column of Table 4 also reports the average thickness of a continuous layer of soil (if any) where a value of the excess pore pressure ratio ru>0.8 was calculated. In Figure 15 the time history of acceleration at the base of the model (grey line) is compared with that calculated at the surface (black line) for two cases from Table 4: 'Norcia' and 'Northridge' input. The achievement of liquefaction in the soil layer can be noticed in both cases.  Initially the signal is amplified (up to about 0.2 g, that is at about 2.5 s for 'Norcia' and 4 s for 'Northridge') at the surface compared to the base, then liquefaction occurs and the liquefied soil acts as an isolating layer: the amplitude of acceleration at the surface is lower than at the base from this point onwards. Figure 16 shows the distribution of the excess pore pressure ratio ru at the end of shaking in both cases. The shading has been limited to the range 0.8  ru  1. It can be observed that a continuous horizontal layer of soil near to the surface is very close to liquefaction if not liquefied. Moreover, the tunnel itself is partially interacting with liquefied soil, although deeper than the shallow liquefied horizontal layer (C/D=2). In Figure 17 the ratio between peak acceleration at the surface and that at the base is plotted against the peak acceleration at the base, for all the input signals shown in Table 4. It can be noticed that in the cases where the peak acceleration of the input signal is lower than 0.2 g, such a ratio is higher than 1, indicating amplification, while for higher values of peak acceleration the ratio is lower than 1, indicating that deamplification occurred. In all cases de-amplification is caused by liquefaction occurring near the ground surface. The depth of the tunnel in the ground layer does not affect the dynamic response of the soil, as shown in the figure.
The effect of soil liquefaction on the tunnel lining is analysed by looking at the maximum changes of hoop force and bending moment at the end of shaking (Table 4 and Figure 18a, b).   Figure 18 indicates that when the ground amplification prevails (for this ground conditions when amax,b < 0.2 g according to Figure 17) larger changes of internal forces arise for increasing amplitude of shaking. This trend is more evident for deeper tunnels (C/D=2). On the other hand, when liquefaction prevails (when amax,b > 0.2 g), the change of internal forces is independent from the amplitude of the base acceleration.
In terms of permanent displacements induced by soil liquefaction, it is worth noting that the calculations were performed by imposing undrained conditions. Pore-pressure build-up during shaking produces a very limited uplift of the tunnel, unless the soil liquefies and the liquefied ground interact with the tunnel, such as in the cases of Figure 16. In Figure 19 the calculated uplift of the tunnel at the end of shaking is plotted as a function of the average thickness of a liquefied layer. This has been assumed to be a shallow continuous horizontal layer with ru>0.8 (see for instance the shaded areas in Figure 16). The trend in the figure shows that large amounts of liquefaction in the cover soil layer of the tunnel produces significant uplift of a shallow tunnel, although the cover upon diameter ratio is not too low (C/D=2). Similar trends were obtained for shallower tunnels (C/D = 0.5 and 1) using the pseudo-harmonic input signal only (that is 'EQ1', '2x', '3x' in Table 4) and are shown in Figure 20. Although the numerical results are limited, the effect of lower overburden can be read. The shallower the tunnel the larger the uplift associated with the mobility of the surrounding liquefied soil, as observed experimentally by Chian and Madabhushi (2011) in the centrifuge.

Remarks
This section has shown how an advanced effective stress constitutive model, able to capture the cyclic behaviour of sand in both drained and undrained conditions, has been adapted to back-analyse a centrifuge test in dry sand in order to model afterwards a similar problem in saturated sand. The constitutive model has been calibrated using the results of laboratory tests carried out in monotonic and cyclic loading.
After comparing the numerical simulation in dry and saturated conditions, the numerical model has been used to extend the study to different conditions in terms of lining thickness, tunnel cover, input signal. This provided an insight into the behaviour of a shallow tunnel in a liquefiable sand layer. A form of limiting threshold to the change of internal forces induced by ground shaking in the tunnel lining has been observed in the numerical results once soil liquefaction occurs. At the same time, the influence of the overburden cover on the uplift induced at the tunnel in cases of extensive liquefaction has been discussed on the basis of the calculations. Due to the lack of existing measurements for real cases, experimental campaigns using centrifuge modelling would be highly beneficial to corroborate or debate similar results.

Background
Although uplift mechanisms for an underground structure experiencing soil liquefaction have been identified experimentally and numerically by several authors, the interaction of such mechanisms and the associated displacements of the underground structure with those induced in aboveground structures that may be founded nearby have not yet been investigated.
In urban areas shallow tunnels are likely to be close to the foundations of buildings and easily interact with them during earthquakes (i.e. Soil-Structure-Underground Structure-Interaction, SSUSI). Hence, the reciprocal influence of a tunnel and an adjacent building in the presence of soil liquefaction may be important.
Recent centrifuge testing on the behaviour of buildings founded in liquefiable ground layers has shown that smaller net excess pore pressures are generated within the liquefiable layer under a structure by increasing the contact pressure and height/width ratio of the building (Karimi & Dashti, 2016). Other studies have shown the reciprocal influence of adjacent buildings, affecting non-uniform settlement during liquefaction (Yasuda, 2014).
How the uplift mechanism of an adjacent underground facility is influenced by the presence of the building and how the floating of the underground structure can affect the tilt and settlement of the building are both aspects that deserve attention.
This problem appears rather important considering the rapid extension of the built environment, both aboveand underground, to areas that may be subjected to risk of liquefaction. Hence an insight into such a problem may well contribute to increase the resilience of urban environment to natural hazards.
The project STILUS, within the framework of the European funded network SERA (Seismology and Earthquake Engineering Research Infrastructure Alliance for Europe) intends to investigate this problem through a series of centrifuge tests. In order to plan the centrifuge tests, a preliminary numerical study of tunnel-structure interaction in liquefiable soil was carried out as described in the following sections 4.2 and 4.3. A circular transverse section (modelling a bored tunnel) and a rectangular framed section (modelling a cut-and-cover tunnel) were taken into account at this preliminary stage, since they may be likely to occur in the urban environment.

Circular tunnel
The numerical calculations were carried out using the same layout as shown in Figure 1, that has been analysed in section 3. The numerical model and the input signal used were the same as described in that section, being the problem modelled using the UBC3D-PLM model for the soil (Galavi et al., 2013). A simple structure was added in the model mimicking a two-storey building as shown in Figure 21 (Colamarino, 2017). The building consists of a two-floors (3 m high each) and a basement (2 m deep). The building rests to one side of the tunnel as shown in the figure.
The building frame was modelled using linear elastic beam elements. Two different material datasets were used, one for the basement (EI = 1.6x10 5 , kNm 2 /m, EA = 1.2x10 7 kN/m) and the other for the rest of the building (EI = 6.75x10 4 kNm 2 /m, EA = 1.6x10 5 kN/m). The mass assigned per unit length to the beam elements takes into account also the presence of the floors and the walls.
The same set of input signals as in the previous section was used in order to compare the results to the "greenfield" conditions considered in that section. Figure 22 shows the highest values of excess pore pressure ratio (ru>0.8) calculated in undrained conditions at the end of the shaking for the three input signals 'EQ1' (Fig.  22a), '2x' (Fig. 22b) and '3x' (Fig. 22c). Insets in the same figure show the corresponding deformed configurations at the end of shaking.
As soon as the amplitude of the signal increases, larger areas of the sand layer are affected by liquefaction or are approaching it (ru>0.8). It is worth noticing that for 'EQ1' (Fig. 22a) the highest values of ru are distributed in the area of maximum shear stresses around the building foundation. Instead, no evidence of liquefaction was observed in the results of the corresponding greenfield analysis in section 3 (ru<0.8, see Figure 13). The larger amplitude of the "2x" input signal (Fig. 22b) produces a continuous layer of shallow soil approaching liquefaction.
The influence of the building is still visible in this distribution but liquefaction does not affect the soil around the tunnel (C/D=2). When subjected to an even stronger shaking, a larger thickness of soil approached liquefaction (Fig. 22c). Compared to the corresponding greenfield analysis, the influence of the stresses induced by the building is evident both in terms of deviator and mean stress: calculated pore-pressure build-up are higher at the corners of the foundation due to initial higher shear stresses and lower towards the centre, where the mean stresses prevail. In this case liquefaction areas reached the tunnel below.
The building settles and tilts. Both settlement and tilt are influenced by the distribution of excess pore pressure around the foundation, that affects the degree of mobilization of shear strength in the foundation ground. However, when liquefaction reaches the tunnel depth, an increased uplift of the tunnel affects the building movements and the building starts to counter-rotate.  Figures 23 and 24 show the time histories of settlement and excess pore pressure calculated at point G and I (see Figure 21) with input signals "EQ1" and "3x". For the weaker "EQ1", the larger positive excess pore pressure that arises around point G (Fig. 23b) produces a larger settlement of the building on that side at the end of shaking (Fig. 23a).
On the other hand, for the stronger "3x", although negative excess pore pressure develop around point I ( fig.  24b), the buildings settles more on the right side ( fig. 24a), indicating an interaction with the uplift of the tunnel on the left side.  Tables 5a and 5b summarize the results achieved in the analyses with C/D=2. Positive tilt is assumed counterclockwise. In Figure 25 the maximum value of uplift of the tunnel axis is plotted for different values of the ratio C/D as a function of the average thickness of a continuous horizontal layer of soil where the excess pore pressure ratio ru is larger than 0.8. As in section 3, such a value has been assumed as a proxy for the effect of liquefaction in the ground layer. In the same figure two curves are shown that represent the trends calculated for C/D=2 and C/D=0.5 in the analyses without buildings (section 3). Although with some scatter, the trends are the same in both sets of analyses (with and without buildings), indicating a minor effect of the presence of a building on the amount of tunnel uplift, providing that similar distributions of pore pressure build-up affect the soil surrounding the tunnel. Trends of increasing building settlement and tilt can be observed in Figure 26 and 27. Very low values of average thickness of the layer with ru>0.8 (close to zero) indicate that liquefaction occurs only in the proximity of the foundation of the building (e.g. Fig. 22a). This corresponds to limited settlement, although non-negligible. Much larger settlement is calculated when liquefaction is approached in larger volumes of soil, as for instance in the cases shown in Fig. 22b-c. Correspondingly, it might be noted that in Figure 27 there is a decreasing trend of tilt towards negative values as the average thickness of the layer with ru>0.8 increases. Hence, a larger pore-pressure build-up generally induces the foundation to rotate clockwise. This indicates the effect of the upheaval associated with the tunnel buoyancy on the left side of the building.  In Figure 28 the change of hoop force (a) and bending moment (b) in the tunnel lining at the end of shaking is plotted as a function of the peak acceleration of the input signal at the base of the model. Trend lines for the case C/D=2 are shown as dashed lines and compared with similar trend lines from Fig. 18 for the 'greenfield' cases, that is without the building (dotted lines). The change of hoop force N induced by pore-pressure buildup is independent of the presence of the building. On the contrary, the change of bending moment is generally larger than in the 'greenfield' case. This finds justification in the less uniform distribution of stresses induced around the tunnel by the presence of the building (compare for instance values of ru: for 'greenfield' conditions in Fig. 16b with 'building" conditions in Fig. 22c).  Table 5 compared to trend lines in Fig. 18 (no building).

Rectangular tunnel
In order to analyse a typical case of a cut-and-cover tunnel in an urban environment, a rectangular section has been assumed, as shown in Figure 29. The same liquefiable sand layer and the same building as in section 4.2 are modelled. Table 6 shows an overview of the analyses that have been carried out. The legend for the input signals has been given in Table 4. A set of analyses with input signals of increasing amplitude was carried out in 'greenfield' conditions (Table 6a), to study the effect on the dynamic response and the pore-pressure build-up of the presence of the tunnel. Similarly, a set of analyses was carried out for models with a building and without a tunnel (Table 6b). Finally, the tunnel-building interaction was analysed with two sets of numerical models with both structures (Table 6c). In the former the building was located on the edge of the tunnel, with a distance to cover ratio, d/C=0 (Fig. 29a). In the latter, the building was located at a distance d=5 m on the right side of the tunnel, corresponding to d/C=1.7 (Fig. 29b).
The tunnel was very shallow, with a cover C=3m, compared to the depth of the basement (2 m).  The dynamic response of the soil layer in 'greenfield' conditions (no building) with and without this tunnel is shown in Figure 30. In the figure the ratio between the peak acceleration at the surface and at the base is plotted against the peak acceleration at the base. It can be noticed that the presence of the larger rectangular tunnel reduces the amplification at the ground surface compared to the circular tunnel. In all cases the dynamic amplification calculated in 'free-field' in the corresponding analyses without a tunnel is much larger. The differences between the curves reduce as the peak ground acceleration increases, when de-amplification occurs due to soil liquefaction, as discussed in section 3.5. Figure 30. Ratio between peak acceleration at surface and at the base vs. peak acceleration at the base: with and without rectangular tunnel and comparison with trend for circular tunnel. The presence of the building affects the distribution of excess pore pressure, as shown in section 4.2. This is confirmed by the analyses without a tunnel as shown in Figure 31a. Here the continuity of the horizontal layer approaching liquefaction is broken below the building due to the higher effective mean stresses. The presence of the tunnel further affects the distribution of excess pore pressure, as can be seen from the comparison between Fig. 31a and Fig. 31b.
It is worth noting that liquefaction is confined most in the free-field areas at the right side of the building and the left side of the tunnel. However isolated liquefied soil volumes are identified below the tunnel, above the tunnel roof and between the tunnel and the building foundations. In the latter area the distribution of excess pore pressure depends on the relative distance d/C.
In Figure 32 the calculated values of building maximum settlement are plotted. Figure 32. Max settlement of the building.  In the figure a trend of increasing settlement of the building with the average thickness of the layer approaching liquefaction is observed, although with some scatter. The presence of the tunnel reduces the calculated settlement of the building. Such a reduction is more evident for the shallower and larger rectangular tunnels than for the deeper and smaller circular tunnels. The calculated tilt and settlement of the building at the end of shaking are plotted one against each other in Figure 33. Despite some scatter, the plot shows a certain degree of correlation among the two quantities.

Remarks
This section has described a preliminary numerical study of tunnel-structure interaction in liquefiable soil. This study has been carried out to identify patterns of deformation that should be expected to occur in centrifuge tests to be carried out in a research project concerning such a problem. The results have clearly shown how the distribution of excess pore pressure induced by shaking in undrained conditions is affected by the presence of the tunnel and of the building.
The relative distance between the two structures (here expressed in terms of ratio d/C of the horizontal distance between the tunnel wall and the building basement upon the tunnel cover) influences the solution both in terms of tunnel lining deformation and of building displacements. Consequence are observed in the distribution of internal forces in the tunnel lining and in the final configuration of the building.
A number of indications for implementing the physical modelling have been suggested by the results of the numerical analyses.
As far as the distribution of internal forces is concerned, since the analyses show that it is influenced by the presence of the building, it would be important to have a large number of measuring points along the tunnel lining, to get an experimental insight into this problem.
Furthermore, since building tilt is expected by the analyses, non-contact laser displacement transducers might be used in the centrifuge to measure such a tilt. They will be then associated to conventional transducers (LVDTs).
Moreover, the distribution of calculated ground movements induced by soil liquefaction may help to define areas where the addition of finer content (down to the nanoscale) may reduce the mobility of the soil. On the experimental side, this show the potential benefit of using digital imaging and particle image velocimetry in centrifuge tests, through a transparent side of the model container.
The numerical results also show that the tunnel uplift is driven by the increase of pore pressure below the tunnel invert and the concurrent reduction of effective stresses (and shear resistance) in the cover. Hence the safety factor against uplift is reduced. This will considered in the layout of the centrifuge tests by deploying transducers to measure pore pressures where the analyses show that a significant build-up may develop.
The calculated distribution of excess pore pressure induced around the tunnel and the building basement may also be useful to identify where mitigation techniques that may locally reduce pore-pressure build-up (for instance: drainage, densification, induced partial saturation) would be most effective against the effects of soil liquefaction and should be implemented in the tests. Nevertheless, it should be pointed out that the results of the numerical predictions should be considered with care. Although the potential of the constitutive model used, and of other models of similar complexity, to predict pore pressure build-up and to identify the occurrence of soil liquefaction has been shown in several studies, their accuracy in predicting the deformation of soil approaching or experiencing liquefaction and large strain is still a matter of study. Hence the need to run tests on physical models, thus achieving an experimental assessment of the behaviour of the soil and interacting structures (tunnel and building) in such conditions. CONCLUSIONS This work has illustrated how numerical calculations can be used in association with centrifuge testing to model different aspects of tunnel behaviour during earthquakes. The scope of the paper has been limited to a few aspects, mainly concerning the change of internal forces in the tunnel lining during shaking and the effect of soil liquefaction. Tunnel-building interaction during shaking in liquefaction prone soil has also been investigated. However, those analysed are only examples of a larger number of applications where an integrated experimental-numerical modelling approach can be followed.
The point of view of this paper is on purpose slightly biased towards the numerical modellers that may benefit of centrifuge tests to calibrate their models. The use of centrifuge testing (and physical models in general) should be considered as complementing laboratory testing on single elements when it comes to study specific aspects of boundary value problems. Indeed the possibility of evaluating numerical models using well-defined and controlled experiments increase the reliability of any numerical study where advanced constitutive models are used.
On the other hand, experimental activities may benefit significantly from a preliminary numerical study that helps to define the scope of testing and the key aspects that the physical model should be able to reproduce. This permits efficient use of resources, possibly reducing the number of experiments, to focus effective efforts on the specified target.
The main achievements of this work are only partial and deserve further investigation. However, they help to show how a combined use of physical and numerical modelling is necessary to analyse earthquake-induced effects on tunnels and other similar subsystems of civil infrastructures