Scaling laws for shaking table testing of reinforced concrete tunnels accounting for post-cracking lining response

Abstract This paper proposes a new set of scaling laws for the study of the post-cracking behaviour of lightly reinforced concrete tunnel linings during 1g shaking table testing. The post-cracking behaviour scaling laws are formulated using two non-dimensional parameters: the brittleness number s, which governs the fracturing phenomenon for unreinforced concrete elements and N P , which plays a primary role for the stability of the process of concrete fracture and steel plastic flow in reinforced concrete elements. The proposed laws allow for the development of an “adequate” experimental model and are validated using numerical analyses of a reinforced tunnel in rock, in both prototype and 1:30 model scale. The adopted experimental set-up is inspired by an existing 1g physical testing campaign on the seismic response of a concrete tunnel in rock and the postulated laws are shown to grant satisfactory similitude between the cracking behaviour of the model and prototype tunnel under two examined earthquake records. The potential of using the proposed laws in 1g tests for Class A predictions of evolving crack patterns in reinforced concrete tunnels is highlighted. The proposed laws are examined under three possible boundary conditions, indicating that both rigid and laminar boxes can still change the behaviour significantly compared to an envisaged free field boundary model. The analysis shows though that for larger soil to lining stiffness ratios, boundary artefacts could be greatly reduced. The present study provides useful recommendations for future 1g tests that did not exist to date, while the proposed scaling laws allow for versatility in the design of novel tunnel lining model test materials.


Introduction
Tunnels are considered to be less vulnerable to seismic damage compared to over-ground structures, mainly because of ground confinement. However, recent seismic events have proven that seismic loading may actually drive the tunnel lining past its elastic regime, inducing non-negligible plastic deformations that may in some cases result to severe cracking. Several cracked tunnels were identified after the 1995 Hyogoken-Nambu and the 2016 Kumamoto earthquakes in Japan (Asakura, 1997;Suzuki, 1996;, the 1999 Chi-Chi earthquake in Taiwan (Wang et al., 2001), the 1999 Duzce earthquake in Turkey (O'Rourke et al, 2001) and the 2008 Wenchuan earthquake in China (Wang & Zhang, 2013;Shen et al., 2014).
Despite the obvious importance of such well-documented case histories, numerical and physical modelling is necessary to derive deeper understanding of the failure mechanisms and to explore the key factors affecting the seismic performance of the soil-tunnel system. Especially when considering complicated nonlinear phenomena, such as cracking of the tunnel lining, numerical modelling can have certain limitations. Physical modelling can be employed to better understand the prevailing failure mechanisms, also serving as a benchmark for validation of numerical analysis techniques. In this context, centrifuge testing has been extensively used to study the seismic performance of tunnel structures (Onoue et al., 1998;Yamada et al., 2002;Billota et al., 2009;Cao & Huang, 2010;Shibayama et al., 2010;Cilingir & Madabushi, 2010;Cilingir & Madabushi, 2011a;Cilingir & Madabushi, 2011b;Lanzano et al., 2012;Tsinidis et al., 2014;Chen et al., 2014;Tsinidis et al., 2016). These studies mainly focused on kinematic soil-tunnel interaction, the effect of strong motion and soil deposit characteristics on the seismic response of the tunnel, or the effect of the tunnel shape. Valuable insights were offered on key aspects, such as the dynamic modes of deformation, the developing soil deformations around the tunnel, the earth pressures exerted on the lining, and the lining internal forces. Less attention has been paid to the tunnel lining replication, which is modelled using aluminium alloy, for complying with an implicit assumption of linear elastic response. To the best of our knowledge, the post-cracking behaviour of the reinforced concrete lining and its effects on tunnel response have not been experimentally addressed, despite the fact that their significance has been proven through numerical analysis (Kampas et al, 2019).
The previously quoted seismic damage to tunnels has motivated an attempt to investigate the development of seismically-induced crack patterns on tunnel linings, by incorporating the post-cracking concrete material response into physical models. Wang et al. (2015) conducted 1 g shaking table tests to identify progressive damage states in unreinforced concrete linings under cumulative earthquake loading. In order to quantify the effect of cracking on the seismic performance of the tunnel, they measured the degradation of the first natural frequency of the soil-tunnel system after each seismic event. In their experimental campaign, the model tunnels were made of "emulation concrete", which has been conjectured (as discussed later on) to be efficient for the simulation of the pre-and post-cracking behaviour of concrete in small scale. Sun et al. (2011), had earlier used the same material to conduct 1 g shaking table tests on the nonlinear seismic response of the portals of two parallel tunnels in rock. Chen et al. (2013), conducted a validation study based on results from standard element tests (uniaxial compression/tension, 4-point bending, etc.) indicating similarity in the stress-strain response and the observed damage (crack) patterns. However, these works neither fully explain the necessary scaling laws used to achieve similitude in reinforced concrete in the post-cracking regime, nor explain the effect of the whole composite system similitude.
To the best of our knowledge, no study has developed appropriate scaling laws for the nonlinear structural behaviour of reinforced concrete tunnel structures subjected to seismic loading, under 1 g conditions. Aiming to bridge this gap, this paper introduces a new set of scaling laws, extending previous endeavours, applicable to lightly reinforced concrete tunnels, accounting for cracking. It is well-known that in order for a small-scale model to accurately replicate the behaviour of the corresponding prototype, certain similitude laws need to be applied. These are a series of relationships between dimensionless ratios formed by corresponding parameters of the model and the prototype, determined by dimensional analysis of the studied problem (Iai, 1989;Meymand, 1998;Muir Wood, 2004). Naturally, depending on the objectives, different similitude rules may be applied. In the case examined herein, where the seismic response of the tunnel is studied through 1 g shaking table testing, the objective is to achieve "dynamic similarity", which as defined by Meymand (1998), implies that homologous parts of the model and prototype systems experience homologous net forces. Moreover, modelling the crack patterns of the tunnel lining through reduced-scale shaking table testing necessitates ensuring similarity between model and prototype in the post-cracking regime, and hence, the behaviour of the model concrete and of the steel reinforcement comprises another important scale-modelling criterion.
The developed scaling laws are validated numerically with the aid of 2D plane strain finite element (FE) analyses using Abaqus (2017), employing a damaged plasticity model for the concrete lining, which is validated against published experimental results of 3-point bending tests on lightly-reinforced micro-concrete beams. The potential of using the proposed scaling laws in 1 g tests for Class A predictions of evolving crack patterns in reinforced concrete tunnels is highlighted, while the paper extends its applicability to future 1 g tests on tunnels, by examining boundary artefacts that may hinder proper simulation of the soil-tunnel behaviour in the tests and possible remedies.

Scaling laws for post-cracking behaviour
For the small scale modelling of concrete structures subjected to inelastic deformations, concrete-like materials have been the default choice to date (Moncarz & Krawinkler, 1981). Various such materials have been used in shaking table experiments of past studies, namely micro-concrete, gypsum-based mortars, cement-based mortars, barite powder or bentonite-concrete mixtures (Knappett et al., 2011;Phansri et al., 2015). Gypsum has been favoured for many years as a representative concrete model material, especially for analyses in the elastic regime, due to its high workability, low Poisson's ratio and linear stress-strain relationship. However, its stress-strain response is substantially different to that of prototype concrete when entering the nonlinear post-cracking regime, while it also exhibits a dissimilarly larger tensile strength.
In order to investigate the applicability of the proposed similitude laws for reinforced concrete linings to reduced-scale 1 g shake table testing, the experimental campaign of Wang et al. (2015) is used as a benchmark. In their experiments, a real tunnel in China served as inspiration for the prototype structure. The tunnel has a horseshoe shape and is founded in rock, while the 1 g experiments were performed at a 1:30 geometric scale down. In the present study, the tunnel configuration is modified by considering a minimal amount of reinforcement (e.g. for non-structural cracking control). Everything else (model scale, material properties, etc.) is kept the same to allow some form of comparison. To this end, emulation concrete is adopted as the model material to simulate the nonlinear behaviour of the concrete tunnel lining.
Developed by researchers in the Dalian University of Technology, emulation concrete is a low-strength, fine aggregate material, made by a mixture of cement, water, barite powder, mineral powder and barite sand in suitable proportions to meet the case by case similitude requirements. The material has been increasingly used in 1 g shaking table tests, aimed at detecting the cracking patterns of concrete structures subjected to dynamic loading, including tunnels (Sun et al., 2011;Wang et al., 2015) and gravity dams (Chen et al., 2013). Its efficacy to model the post-cracking unreinforced concrete response is demonstrated experimentally by Chen et al. (2013). Emulation concrete contains only small quantities of cement to achieve strength similitude, and the sizes of its aggregates are calculated according to the geometric scale. The particle sizes of barite, barite sand, and ore powder range from 0.05 mm to 2 mm (Chen, et al., 2013).
Having assumed the use of a model material of proper stress-strain behaviour for concrete, our focus lies on the definition of an applicable set of 1 g scaling laws for modelling lightly reinforced concrete tunnels. To this end, this study makes use of two dimensionless parameters, s and N P . The brittleness number s, governs the fracturing phenomenon for brittle unreinforced concrete elements and is given by the following equation (Carpinteri, 1982): where: K IC is the material fracture toughness; f t the ultimate tensile strength of concrete; and h a characteristic linear size of the specimen. Carpinteri (1982) used dimensional analysis and the parameter s to demonstrate that for the accurate reproduction of the mechanical response of a concrete structure and its collapse mechanisms during fracture tests on scaled specimens, the model size should not be too small, or more specifically, the brittleness number s should not exceed the limit values of 0.5 and 0.54 for 3-point bending tests and tension tests, respectively.
The dimensionless number N P applies to lightly reinforced concrete elements and, according to the study of Carpinteri (1984), it plays a primary role for the stability of the process of concrete fracture and steel plastic flow. It is defined as follows: where: f y and ρ t are the yield strength and the percentage of longitudinal steel reinforcement within the cross section, respectively. According to Corrado et al. (2011), N P values close to zero correspond to plain concrete (e.g. N P ≤ 0.12) that will lead to an unstable reinforced M. Antoniou, et al. Tunnelling and Underground Space Technology 101 (2020) 103353 concrete fracture process. Corrado et al (2011) applied dimensional analysis using the parameters s and N P to study the flexural behaviour of lightly reinforced concrete beams, taking into account the nonlinear behaviour of concrete in tension and considering an interface that allows slippage between concrete and steel reinforcement. They concluded that physical similitude for the structural response of reinforced concrete beams (expressed in terms of dimensionless flexural moment versus normalised cross section rotation), can be achieved only when both s and N P are kept constant for different specimen sizes. Based on their observation and according to Eqs. (1) and (2), the following similarity scales are derived for parameters s and N P : where the subscripts m and p denote properties of the model and of the prototype, respectively. Equations (3) and (4) yield: where: S f t , S f y and S ρ t are the similitude ratios for the concrete tensile strength, the steel yield strength, and the reinforcement percentage, respectively. Equation (7) shows that for a prototype of known reinforcement yield strength, reinforcement ratio and tensile strength, an appropriate model-scale material can be found that fulfills the equation. It should be noted that although steel is normally referred to as the reinforcement material, there is no reason why a different material could not be used as long as it complies with Eq. (7), which offers great possibilities to be imaginative in developing future novel model-scale materials. Finally, one can observe that Eq. (5) results in the following similitude relationship for concrete fracture toughness: where S KIC and S L are the similitude ratios for concrete fracture toughness and length respectively. The postulated post-cracking similitude laws refer to lightly reinforced concrete sections, where the maximum moment capacity owing to the concrete tensile strength is comparable to the moment capacity as calculated based on the steel (or other) reinforcement contribution. In practical terms, this translates to elements reinforced with close to the minimum codified reinforcement percentage. The latter, in most cases, is required in order to stabilize the crack formation and avoid phenomena of brittle failure.
Without being able to explicitly define the crossover point from lightly to heavily reinforced, as the mechanical response of reinforced concrete during failure depends on both the steel percentage and the concrete section properties, we trust that the applicability of our methodology extends in all cases where the reinforcement lies close to widely-used minimum required reinforcement predictions/criteria. On the other hand, a lower limit of reinforcement can be found using the formula N P = The proposed scaling provides a viable approach (Meymand, 1998), meaning that we basically identify and successfully model the primary forces and processes in the system, while suppressing secondary effects (discussed later on). It therefore constitutes an "adequate" model, rather than applying similitude theory to achieve a "true" (1 to 1) model similarity. Table 1 in conjunction with Eq. (7) summarize the full set of scaling relationships, i.e. the applied elastic laws and the post-cracking relationships proposed herein. The following sections are dedicated to the application of the proposed scaling laws to the previously discussed experimental model, aiming to demonstrate their validity and to propose a robust methodology for the assessment of tunnel seismic response employing 1 g shaking table testing.

Numerical analysis methodology
The seismic response of the previously discussed tunnel is analysed employing the finite element (FE) analysis environment of Abaqus 6.17. The main focus of the analysis is to explore the validity of the developed scaling laws.

Constitutive modelling of concrete
The Concrete Damaged Plasticity (CDP) constitutive model is used to simulate the nonlinear response of reinforced concrete (Fig. 1). The CDP model is a continuum, plasticity-based, damage model, which is available in the Abaqus FE environment. It follows a yield condition with a non-associated plastic flow rule, based on the yield functions proposed by Lubliner et al. (1989) and Lee & Fenves (1998). The Drucker-Prager hyperbolic function is used as the flow potential. The model was proposed for monotonic, cyclic and/or dynamic loading of concrete structures. It has subsequently been demonstrated to adequately capture dynamic material response (e.g., Krahl et al., 2018;Benham et al., 2018;Zoubek et al., 2013;Nagy et al., 2010), also including buried structures (Nagy et al, 2010).
It aims at capturing the effects of irreversible damage associated with the two main failure mechanisms that occur in concrete: tensile cracking and compression crushing. The stress-strain relations are governed by scalar damaged elasticity (SIMULIA, 2014): where f is the Cauchy stress tensor; ε is the strain tensor; D el 0 is the initial (undamaged) elastic stiffness of the material; −  Antoniou, et al. Tunnelling and Underground Space Technology 101 (2020) 103353 degraded elastic stiffness; ε pl is the plastic strain tensor; and d is a scalar stiffness damage variable, which can take values in the range from zero (undamaged material) to one (fully damaged material). For a given cross-section, the factor − d (1 ) represents the ratio of the effective load-carrying area (i.e., the overall area minus the damaged area) to the overall section area. To model the different damage states for tensile and compressive loading, the damage variable d in Eq. (9) is defined using d t and d c , which are two independent state variables: d d ( t , d c ). Model response to uniaxial loading in tension and compression is schematically shown in Fig. 1. Damaged states in tension and compression are characterized independently by two hardening variables, ε t pl and ε , c pl which are referred to as equivalent plastic strains in tension and compression, respectively. The cracking and crushing in the concrete are represented by increasing values of the hardening variables, while they also control the evolution of the yield surface and the degradation of the elastic stiffness. The damage variables in uniaxial tension and compression (d t and d c ) are also functions of the concrete tensile/compressive strengths ( f tm and f cm ), fracture and crushing energies per unit area (G F and G ch ) and a characteristic element length (l eq ) of the mesh, namely: The description of nonlinear behaviour of a concrete specimen requires definition of the stress-strain relationship under uniaxial tension and compression, as well as the evolution of the damage variables d t and d c .
The CDP model requires the following four material constants: K c , ψ, f bo /f co and ∊. Parameter K c is the ratio of the second stress invariant on the tensile meridian q(TM) to that of the compressive meridian q(CM) at the initial yield (0.5 < K c ≤ 1.0), ψ is the dilation angle, f bo /f co the ratio of initial biaxial compressive yield stress to the initial uniaxial compressive yield stress, and ∊ the flow potential eccentricity. The four parameters used in this study are adopted from closed-form expressions of Alfarah et al. (2017), listed in Table 2. The dilation angle ψ = 13°is a lower bound for concrete (it may as well be as high as 50°), denoting that a less stiff response is taken into account, while all other parameters equal the default values suggested by Abaqus (SIMULIA, 2014). It should be noted that strain rate and long-term effects were not accounted for.
The use of damaged-plasticity laws predicting strain localization with finite elements suffers from well recognized mesh dependency effects, in such a way that the dissipated energy decreases with mesh refinement. To alleviate this limitation, Alfarah et al. (2017) incorporate the characteristic element length l eq in the calculation of the damage variables and the yielding parts of the stress-strain curves, to allow for proper simulation of the dissipated energy with respect to the utilized FE mesh. Before proceeding to the analysis of the soil-tunnel system, the numerical model is validated against well-documented three-point bending tests of reinforced concrete model beams. Ruiz et al. (1998) conducted a series of 3-point bending tests of lightly reinforced micro-concrete beams to examine fracturing phenomena with respect to size effects. Beams with a reinforcement percentageρ t = 0.13% are selected for the validation conducted herein, considering two different values of section thickness: D = 300 mm, case (a), and D = 150 mm, case (b). Both beams are of width b = 50 mm, while ribbed steel wires of 2.5 mm nominal diameter are used as reinforcement [two wires for case (a) and four wires for case (b)],for ensuring better bond between micro-concrete and steel. 2D FE analyses are conducted, considering the out-of-plane dimension equal to the beam width. Concrete is modelled with nonlinear quadrilateral continuum elements, while the steel reinforcement is simulated with a single row of truss elements with zero torsional stiffness, whose area corresponds to the entire reinforcement area. The properties of microconcrete and steel wires used for the validation are summarized in Table 3. The FE models are discretized in square elements of uniform length l eq = 10 mm for case (a), and l eq = 5 mm for case (b). The   Table 3 Material properties of the micro-concrete and reinforcement steel bars used in the three-point bending tests FE simulation (Ruiz et al,1998 (Ruiz et al., 1998). Loaddisplacement (P-δ) curves of the reinforced beams are presented for the cases of: (a) D = 300 mm, ρ t = 0.13%, (b) D = 150 mm, ρ t = 0.13%. In (c), ABAQUS snapshots illustrating the tensile crack pattern of the RC beams of D = 300 mm (on the left) and D = 150 mm (on the right) are displayed. derivation of stressstrain curves under uniaxial tension and compression, and the evolution of the damage variables d t and d c is conducted with respect to the characteristic lengths for each case. The FE model involves simulation of the cylinder supports, by considering contact elements that allow sliding. A Coulomb friction coefficient μ = 0.05 is assumed, representing the expected low friction between the rigid cylinder supports and the model beams. The loading protocol includes increments of displacement at a rate of 0.3 mm/min, applied at the top centre of the beams until failure. Representative FE analysis results are presented in Fig. 2, in terms of load-displacement curves, revealing good agreement with the experimental measurements for both beam configurations, not only in terms of initial stiffness and maximum load capacity, but most importantly, in terms of post-peak behaviour. Note that after reaching the initiation of cracking and the peak value, the beam resistance decreases rapidly until the fracture zone reaches the reinforcement. Then, it increases again as fracturing is retained by the steel wires, and thus hardening due to composite action takes over. Fig. 2c displays snapshots of the tensile crack pattern developed on the beams by the end of loading.

Modelling of an reinforced concrete tunnel in rock
As previously discussed, the geometry of the examined reinforced concrete tunnel is inspired by a real tunnel in China. Located in an area of high seismicity between the towns of Mianning and Chengxiang, the Lebuguo Tunnel, has a horseshoe shape with excavated dimensions of 12.64 m width and 10.14 m height (Fig. 3). Focusing on the ovalingracking response of the tunnel due to shear waves propagating normal (or nearly normal) to its axis, the analysis is conducted assuming planestrain conditions. This is typically the case for several previous studies that focused on the same mode of deformation, involving either numerical or centrifuge modeling (e.g., Billota et al., 2009;Tsinidis et al., u Fig. 4. General configurations of the employed FE models: (a) Model A replicates setup within a rigid box, with the tunnel located at 3D distance from the rigid boundary, (b) Model B replicates setup within a laminar box, with the tunnel located at the same distance and (c) Model C intends to eliminate the lateral boundaries effect, placing them as far as 16D from the tunnel centerline.  2016). Naturally, there are certain limitations as the following mechanisms cannot be captured: (a) axial compression/extension stemming from the seismic wave components that propagate along the longitudinal tunnel axis of the tunnel (wave passage effects), and can result to significant damage, especially in the case of tunnels with joined segments (Anastasopoulos et al., 2007); (b) the longitudinal bending due to incoherence effectsalthough these are considered minor compared to wave passage effects (Anastasopoulos et al., 2007); (c) the damage associated with the response of tunnel portals. Such mechanisms require a more detailed 3D analysis. An initial sensitivity study is performed with respect to the lateral boundaries of the model. Fig. 4 shows snapshots of the FE meshes. Model A (Fig. 4a) corresponds to the geometrical replication of the shaking table test experimental setup within a considered rigid box utilized in the study of Wang et al. (2015). Models B and C refer to different simulations of the lateral model boundaries, aiming to examine their effect on soil-tunnel response. More details about the boundary conditions are offered in the following sections. The characteristic element length is set equal to l eq,p = 0.05 m for the prototype tunnel and subsequently scaled to l eq,m = 0.0017 m for the 1 : 30 model tunnel (according to the proposed length similitude law).

Material response
The properties of the prototype concrete lining and steel reinforcement are listed in Tables 4 and 5, respectively. Typical values for a Class 25/30 concrete mix are adopted, assuming damping ratio ξ = 5%. A lightly reinforced prototype is considered with ρ t = 0.169%. We assume one array of steel bars at the outer surface of the lining, where preliminary analyses showed that initiation of cracking is expected. More specifically, the prototype section embodies Ø22 (dp p = 22 mm) rebars placed at 0.1 m from the tunnel outer face (minimum cover depth of 5 cm, according to ITA Guidelines for the Design of Tunnels, 1988) and spaced at Sp a,p = 0.5 m along the tunnel longitudinal axis (Fig. 3). When reinforcement is provided for crack control rather than covering inner stresses, as for the case examined herein (the Lebuguo Tunnel is described as unreinforced in Wang et al., 2015), ITA guidelines propose 1.5 cm 3 /m of steel reinforcement at the outer lining surface. The assumed reinforcement of Ø22 per 0.5 m corresponds to 7.6 cm 3 /m of steel reinforcement at the transverse lining direction.
Concrete is modelled with nonlinear continuum quadrilateral elements, while truss elements are embedded within the concrete to model the steel rebars. The emulation concrete considered for the scaled-down model has compressive and tensile strengths f c = 0.67 MPa and f t = 0.067 MPa, respectively, Young's modulus E = 660 MPa, and density ρ = 1.67 t/m 3 (Table 4). These values lie within the limits found in the literature regarding the properties of emulation concrete (Chen et al., 2013;Wang et al., 2015).
An elastic-perfectly plastic constitutive law is adopted for the steel rebars. Likewise, the steel wire used as reinforcement in the scaleddown models is reproduced by setting the steel yield strength (f y ) and stiffness (E) to the typical values of Table 5. The required percentage of steel reinforcement for the small-scale model may be calculated with respect to Eq. (7). Moreover, selecting the desired wire size, here dp m = 0.2 mm, allows calculation of rebar spacing, according to the following equation: A perfect-bond interface is introduced between the concrete and rebar element, not allowing sliding or separation. In order to follow the similitude law of Eq. (7) for post-cracking behaviour, the dimensions, density, elastic modulus, and yield stress of the wire material do not follow the elastic scaling laws of Table 1. For example, the reinforcement yield stressalthough being stress in magnitudeis considered as an additional, different state-variable that follows different scaling. The scaling derives from non-dimensionalising the post-cracking similarity governing equation (i.e. Eq. (7)), which defines additional internal variables that are expected, at some level, to contradict the elastic scaling laws. However, the new reinforcementfocused similitude conditions become dominant only within the cracking/plastic regime, where one would expect the reinforcement contribution to become the governing factor of the overall structural behaviour. The difference in scaling values is readily observed in Table 5, where the actual scaling values for the reinforcement material (stemming from the needs of the post-cracking similitude laws) are compared to the respective scaling values applied to concrete and soil materials (which are based on the elastic laws).
The required stress-strain curves for the concrete behaviour are deduced for the prototype concrete and then scaled-down according to the relevant similitude laws to describe the response of the emulation concrete. This is demonstrated in Fig. 5 by comparing the dimensionless compression stress-strain curves adopted in this study and those  derived from element tests on emulation and prototype concrete by Chen et al. (2013). The soil employed in the FE analyses is modeled with quadrilateral plane-strain continuum elements. Its stress-strain response displays an elastichardening plastic behaviour, defined by Mohr-Coulomb plasticity and a constant Young's modulus E. The hysteretic soil behaviour obeys the Masing rule, so that the soil stiffness of the unloading/reloading branch equals the initial stiffness value of the backbone curve. The constitutive model works with a single value of small-strain elastic modulus, and thus, is not able to realistically reproduce the complex degradation of soil stiffness due to increase in plastic strain. However, its simplicity serves well the purpose of this study, as the emphasis here lies on the verification of the proposed laws with respect to the structural behaviour of the tunnel lining and not the replication of an actual soil -tunnel interaction problem. The properties of the soil in the smallscale model are adopted from Wang et al. (2015), in order to have a reference to a tested laboratory material. In prototype scale, the soil properties correspond to a tunnel founded in rock. Both model and prototype soil properties are listed in Table 6. Rayleigh damping is introduced in the model according to the classical relationship: where ξ is the Rayleigh damping ratio, ω is the system's angular frequency, a R and b R are the Rayleigh damping coefficients proportional to the mass and stiffness respectively. The coefficients a R and b R are derived so that ξ corresponds to 5% around the first and second undamped natural frequencies of the system; therefore different coefficients where calculated for model and prototype. We hereby select this value as rather typical to match the modal damping value at the first and second coupled modes of our soil-structure system, whose structural component consists of a cracking concrete structure. Moreover, the adopted ξ = 5% value allows for the viscous behaviour of soil in the elastic regime, thereby compensating for the fact that the adopted soil model assumes linear isotropic elasticity prior to yield, and, therefore, does not allow for any hysteretic damping before reaching the Mohr-Coulomb plastic state. A similar practice is adopted in the numerical study of Bilotta et al. (2009), who investigated the seismic performance of shallow tunnels: in their case, an elastic soil model was considered in the FE analyses, and increased Rayleigh damping was incorporated in order to simulate soil hysteretic behaviour. Special zero thickness contact elements are introduced between the soil and the lining to simulate a tensionless interface, allowing sliding and separation. The tangential interface follows the Coulomb friction law, with a friction coefficient μ. According to relevant studies (Sederat et al., 2009;Tsinidis et al. 2016), interface conditions significantly affect the tunnel response during seismic shaking, with a full-slip interface (μ = 0) resulting to reduced developing shear stresses on the tunnel sidewalls, compared to no-slip conditions. A value of μ = 0.5 was selected herein as representative for the concretesoil interface; it is slightly lower than the tan(φ) value (φ: internal soil friction angle) recommended by ASCE-ALA guidelines (2005) for concretesoil interfaces, when considering axial static pipeline movement relative to the surrounding soil.
Additionally, several authors (e.g., Hashash et al., 2005;Hashash et al., 2001;Kontoe et al., 2014;Penzien, 2000;Pitilakis and Tsinidis, 2014;Tsinidis, 2017;Zhang and Liu, 2018) have shown the importance of the flexibility ratio F in characterising the interaction between structure and soil. In particular, F as defined by Wang (1993) and Hashash et al. (2001) should be the same for the prototype and the model to guarantee similarity. Using values from Tables 4 and 5 confirms that the flexibility ratio F is 96.3 for both the prototype and the 1 g models.

Boundary conditions
With the intention to apply the postulated scaling laws in a range of plausible boundary conditions used in shaking table tests, the three cases of Fig. 4 are examined: two more realistic (in terms of size) experimental setups, referring to the use of a rigid and a laminar box, and a very wide model (box) to reduced boundary effects.
As previously quoted, Model A (Fig. 4a) is representative of shaking table tests conducted within a rigid box. The experimental campaign of Wang et al. (2015) is used as an example. The rigid box utilized in their 1 g tests, conducted at the Southwest Jiaotong University (China), is 2.5 m × 2.5 m × 3 m (length × width × height). In this model, the tunnel centerline is located at approximately 3D distance (D = tunnel width) from the rigid box boundary. To minimize boundary effects (reflection and scattering of waves on the boundaries), EPS geofoam material of 10 cm thickness was placed at both side walls of the model box. EPS geofoam is a low stiffness, lightweight material, aiming to minimize boundary effects. It has a density ρ = 0.016 t/m 3 , elastic modulus E = 4.7 MPa and Poisson's ratio ν = 0.09 . The respective numerical simulations have accounted for both the material response of the geofoam (assumed elastic) and the soil-geofoam interface response, characterized by μ = 0.27 in accordance to the range proposed by .
Model B (Fig. 4b) is the equivalent of a shaking table test conducted within a laminar box. For the sake of comparison, the box dimensions remain the same as Model A, while the geofoam material is substituted with soil. Proper kinematic constraints ('node-to node' multipoint constraint, MPC, pins that are forcing two nodes to have identical displacements) are applied along the model boundaries to simulate the response of a soil column subjected to in-plane vertically propagating waves. The same method has been previously applied by Bilotta et al. (2009) and Tsinidis et al. (2016) for the numerical simulation of centrifuge tests on tunnels, conducted in a laminar box and an Equivalent Shear Beam (ESB) container, respectively. M. Antoniou, et al. Tunnelling and Underground Space Technology 101 (2020) 103353 Model C (Fig. 4c) is employed in an attempt to diminish boundary effects. The lateral boundaries are placed remotely enough, i.e. at a distance of 16D from the tunnel centreline, while using the same flexible boundary kinematic constraints as in Model B. In reality, beyond laborious (see e.g. box filling considerations) such a configuration is difficult to achieve rigid box like stiffness.

White noise excitation for eigenfrequency identification
Prior to studying their response to real earthquakes, the FE models were subjected to low-amplitude, band-limited white noise vibrations in order to estimate their dynamic characteristics. To this end, a zero mean random Gaussian signal within the 0-50 Hz frequency band and standard deviation of 0.01 g is employed. Fig. 6a portrays the acceleration time-history and the power spectral density of this signal in prototype scale.
Six different cases are modelled: three different boundary conditions (i.e., Models A, B and C) in two scales (prototype and 1:30 model). Geometric nonlinearities are also taken into account by activating the NLGEOM option available in ABAQUSthis option enables the FE analysis to distinguish between reference/original (undeformed) and current (deformed) configurations, so that stresses and strains are calculated per current element area, in view of the plastic deformations expected to develop in the examined problem. The acceleration time histories and respective response power spectra are derived with respect to the tunnel invert (Point A in Fig. 6). Estimation of the systems eigenfrequencies is based on the transfer function S, i.e. the power spectrum calculated for the tunnel invert (Point A) divided by the M. Antoniou, et al. Tunnelling and Underground Space Technology 101 (2020) 103353 power spectrum of the base excitation (Point I). Results are plotted in Fig. 6b for the prototype and the scaled-down models, indicating that similitude is achieved in terms of system initial dynamic properties. With respect to Model A, the first and second natural frequencies of the prototype are f p 1, = 4.0 Hz and f p 2, = 6.4 Hz, while the respective values for the 1:30 model are f m 1, = 21.9 Hz and f m 2, = 35.1 Hz. Evidently, the 1: 5.5 similarity requirement (Table 1) is satisfied. The same conclusion applies also for Models B and C.
The first natural frequency of the soil-laminar box system (Model B) without the presence of the tunnel structure is also shown on the respective graph as f NT 1, = 1.85 Hz. Interestingly, this value is slightly higher than the f p 1, = 1.81 Hz value of the Model B, illustrating that the presence of the tunnel structure slightly decreases the stiffness of the soil system, due to the excavated soil, rather than increasing it, as one might expect due to the presence of the stiffer lining. This is, of course, a matter of the relatively high stiffness of the considered soil profile.
The most interesting conclusion stemming from the white noise excitation analyses is the marked deviation of the Model A vibrational characteristics in comparison to the other two systems, which highlights the significance of boundary conditions even though the existence of foam. For Models B and C (Fig. 6b), where the soil movement is not restrained by a rigid box, response is dominated by the natural frequency of the soil column. Note that in both cases f 1 is practically equal to = ≈ f V H /4 s s 1.85 Hz, which is the frequency at first resonance of a soil column subjected to vertically propagating shear waves (Kramer, 1996). On the other hand, Model A displays a substantially stiffer response. The exact same observation holds true also for the second natural frequencies of the systems (i.e., frequency decreases from f p 2, = 6.4 Hz in Model A to f p 2, = 5.5 Hz and f p 2, = 5.4 Hz in Models B and C, respectively). Naturally, this difference is bound to play a significant role in the subsequently discussed response to real earthquakes.

Similitude
The objective of this section is to investigate the effectiveness of the proposed scaling laws, by subjecting the 1:30 models and the equivalent prototypes to recorded seismic motions. Ground acceleration records from two seminal earthquake events are selected for this purpose. Both have been recorded in rock strata: (a) the North-South acceleration component (Component 0) from Kobe University in Japan, recorded during the M6.9 Kobe earthquake in 1995; and (b) the East-West acceleration component (Wolong EW) from Wolong station in the Sichuan province, China, recorded during the M8.0 Wenchuan Earthquake in 2008. Their acceleration time-histories, along with the respective elastic (ξ = 5%) response spectra are plotted in Fig. 7. The Kobe University record, being rich in spectral peaks and having a peak ground acceleration (PGA) equal to 0.29 g, may be regarded as a medium intensity excitation. The station's distance from the epicentre of the earthquake was registered as 24 km, according to the Center for Engineering Strong Motion Data website (https://strongmotioncenter. org/). On the other hand, the Wolong recordrecorded at an epicentral distance of 23 km (Chen & Booth, 2011) -is considerably stronger, with a PGA value of 0.96 g (Yu et al., 2012). It is selected due to the proximity of the examined Lebuguo Tunnel to the epicenter of the 2008 Wenchuan earthquake, but also due to the extensive damage that this earthquake caused to tunnels in the area (Wang & Zhang, 2013).
According to Cilingir & Madabhushi (2011a,b), the magnitude of the maximum input acceleration plays a crucial role on the maximum and residual forces acting on a tunnel section during shaking. Hence, the different intensities of the selected earthquake motions are expected to excite different system responses, allowing assessment of the postcracking tunnel response for both low-to-medium and intense seismic shaking. The proposed methodology is believed to be more suitable for the assessment of plastic response of tunnel linings subjected to earthquake loading, compared to existing studies that can only capture residual moments and forces on the lining, as the experimental setup is based on elastic scaling laws. The resulting deformation and stress of the tunnel lining due the applied seismic loading shall be investigated thoroughly in the following sections.
The FE analyses are conducted in two consecutive steps: (a) a static gravity step, where loading is applied simultaneously to the entire soiltunnel system; and (b) a dynamic step, where the earthquake time history is applied to the model. It is important to note that the gravity step simulates the cast-in-place construction of a model tunnel in the laboratory and, hence, it does not account for staged tunnel construction including excavation. According to the proposed scaling laws of Table 1, a dynamic time step ofdt p = 0.05 s for the prototype translates to dt m = 0.0091 s for the 1:30 model, to maintain similarity.
Figs. 8 and 9 compare the acceleration time-histories at the tunnel invert and the soil surface for the prototype (α p ) and the 1:30 model However, within the range of high amplitude accelerations (i.e., between 3 and 13 s for the Kobe Record; and 5 -20 s for the Wolong Record), an apparent deviation in response between the model and the prototype is observed for the high amplitude acceleration peaks. For the Kobe record, this difference reaches a maximum value of α Δ = 0.15 g at the tunnel invert and α Δ = 0.18 g at the surface (Fig. 8a), while for the stronger Wolong record it momentarily reaches α Δ = 0.37 g at the tunnel invert and α Δ = 0.65 g at the surface (Fig. 9a). The difference should be attributed to the nonlinear behaviour that has commenced during high amplitude response. Namely, although both model and prototype systems commence with close to identical transfer functions (i.e. linear behaviour), they progressively deviate due to even minor scaling differences in their inelastic regime. Deviations accumulate with time and this could probably explain why the longer and stronger Wolong record shows greater differences when compared to Kobe. Moreover, it is interesting to notice how the aforementioned deviations become much less significant when performance is compared in terms of model and prototype response spectra at the same positions (Figs. 8b & 9b). This is particularly clear for the Wolong record (Fig. 9b), as the highest Δα values appear for accelerations other than the maximum values utilised for the construction of the respective response spectrum. Having substantially lower inertia in comparison to the surrounding soil, the tunnel is predominantly subject to kinematic loads during shaking, as repeatedly reported in several studies (e.g., Paolucci & Pitilakis, 2007;Billota et al., 2009;. As a result, its response can be directly associated to soil deformations. Fig. 10a b compare the response of the prototype and the 1:30 Model A systems, in terms of dynamic ground settlements (w) at two characteristic locations, namely the invert (w 1 ) and the crown (w 2 ), for both of the aforementioned seismic motions. They indicate that the 1:30 scaled-down physical model reproduces the transient ground displacements experienced in the prototype with sufficient accuracy, irrespective of the severity of the seismic motion. The model also effectively reproduces the kinematic distress imposed upon the tunnel in terms of drift ( u Δ 1,2 ) (i.e., the difference of horizontal movement between crown (u 2 ) and invert (u 1 ), divided by the height of the tunnel) as shown in Fig. 10c. The permanent vertical and horizontal deformations of Fig. 10 are partially caused by the development of plastic points within the soil, which for the case of Kobe University record mainly occurred around the tunnel invert for all three models, as revealed by the soil plastic deformation contours of Fig. 11b. For more accurate soil movement predictions, centrifuge modelling or a more complex soil constitutive model would be required to allow for appropriate scaling of soil stresses and plastic behaviour respectively.

Boundaries
Although ensuring adequate reproduction of the nonlinear response of both the soil and the structure is a major concern, it is not the only prerequisite for realistic physical modelling of such soil-structure Fig. 13. Effect of the soil stiffness relative to the lining material on the tunnel convergence with respect to boundary effects. Note that results refer to the Kobe University excitation, for a linear elastic tunnel with E lining = 29.7 GPa (prototype scale). M. Antoniou, et al. Tunnelling and Underground Space Technology 101 (2020) 103353 interaction problems. Boundary effects also play an appreciably critical role, as highlighted by the snapshots of FE meshes of Fig. 11a. As revealed by the contour plots of the settlements accumulated after excitation of models A, B and C with the Kobe University record, there is a twofold problem: (i) the distance between the tunnel and the boundaries (i.e., the model size); and (ii) the "flexibility" of the boundaries. Both affect significantly the distribution of permanent ground deformations and hence the kinematic distress of the tunnel. It should be noted here, that in the experimental study of Meguid & Mattar (2009), who examined the effect of soil-pile-tunnel interaction in clayey soils through 1 g testing, the distance between the tunnel circumference and the lateral boundaries of the utilized rigid box was set to 4D, in order to minimize their effect on the measured lining stresses.
Boundary effects are quantified in Fig. 12 in terms of stresses and strains developed in the soil, at critical locations near the tunnel. Model C, where free field (flexible) boundaries are placed far away from the tunnel at a distance equal to 16 times its diameter (D) from the centreline, can be considered as the benchmark that most accurately approximates reality. However, it is admittedly impractical, and often impossible (see previously quoted stiffness issue), to use such a large physical model in an experimental campaign. It is, therefore, worth investigating the potential of using models that are smaller in size, if reasonably accurate results can be obtained.
Comparison of the seismically induced normal and shear stresses at the boundaries (Fig. 12ab) indicates that Model A deviates substantially from other responses, irrespective of the earthquake intensity. M. Antoniou, et al. Tunnelling and Underground Space Technology 101 (2020) 103353 As would have been expected for a model with boundaries that do not reproduce free field conditions, the presence of the rigid box in Model A leads to a stiffer system response (and subsequently to lower plastic strain levels within soil), which in turn yields significant under prediction of tunnel response (Fig. 12c). Results dispute the effectiveness of the 10 cm thick geofoam material used at the sides in minimising boundary effects, according to Wang et al. (2015). Conversely, thanks to its flexible boundaries, Model B seems to provide reasonably close estimates of the stress and strain field around the tunnel during the medium intensity earthquake (Kobe University), yet its accuracy reduces when the stronger Wolong Station motion is considered (observe the difference in developed strains in Fig. 12c). However, looking closely at the horizontal displacement profiles of the lateral boundaries in Model B during the Kobe University record (Fig. 11a), a different conclusion appears. A distorted soil response is revealed for Model B, indicating that the 3D lateral distance between the boundaries and the tunnel is insufficient. This results in the formation of a shear band at the lower left region of the model (Fig. 11b) that does not inhibit the model from adequate prediction of the normal stresses at the tunnel crown and the maximum values of shear stress/strain in the vicinity of the tunnel. Such soil shearing alters the overall dynamic characteristics of the soiltunnel system, and as illustrated in the next section, it actually leads to a less similar response to Model C than initially assumed. Fig. 13 displays comparative results on the convergence at the middle tunnel section (i.e. variation of the tunnel height = − δ w w w crown invert ) for Models A, B and C with respect to the ratio of soil stiffness to the lining stiffness (E E / soil lining ). Results correspond to a set of numerical analyses where the soil stiffness is the only parameter varied, to result in different E E / soil lining ratios with respect to the initially assumed. The figure provides encouraging evidence that boundary effects can be probably reduced for a range of E E / soil lining ratios; in the examined case when the latter is larger than 6. Fig. 7a and b showed the initiation of cracking behaviour for Model A (rigid box configuration), which is marked on the acceleration time histories of the Kobe University record and the Wolong Station record. Α limited elastic region is observed before the reinforced concrete lining enters its nonlinear regime due to shaking. The first eigenperiods of all three models are then marked in the elastic spectra of Fig. 7c and d. As observed, they all lie at a period range associated with high acceleration amplification: around 0.5 g for the Kobe University record, and higher than 1 g for the Wolong Station record (approximately equal to 2.7 g for Model A and 1.3 g for Models B, C).

Crack patterns
The evolution of damage is investigated in Fig. 14 Fig. 14) conforms to the previously discussed observations on boundary effects: Models B and C display greater permanent damage in comparison to Model A for both excitations, while the difference in response (in terms of tensile damage) for Models B and C under the Kobe University excitation is clear. Yet, it is worth noting that the mean value of tensile damage may be either slightly over-predicted or under-predicted by the model in comparison to the prototype, depending on the characteristics (e.g., the frequency content) of the applied motion.
Numerical results for the moderate intensity Kobe motion show slight to moderate damage for the tunnel lining, with few cracks detected at the tunnel invert and sidewalls at the end of shaking (Fig. 15). On the contrary, the tunnel suffers extensive damage under the very strong Wolong record, where the tunnel is subjected to accelerations in excess of 1 g (Fig. 16). A plethora of deep tensile cracks develop along the entire lining section, while the previously unaffected tunnel crown and shoulders are also now experiencing severe cracking. It is interesting to observe that the most vulnerable part of the tunnel crosssection lies at the region between the invert and the sidewalls. This part is consistently the first to crack during both records, something that can be attributed to the increased residual stresses it experiences prior to seismic shaking (maximum radial stresses equal to σ r,max ≈ 2500 kPa, i.e. very close to the concrete tensile stress, appear at these points due to gravity loading alone).
The evolution of cracking is best viewed in Fig. 17, which portrays the cracked tunnel section of Model A, at different time frames during the Wolong excitation. As already mentioned, the first cracks appear at the sidewallsinvert connection around t = 6 s, when the tunnel invert is experiencing approximately 0.5 g of acceleration. Between t = 9-10 s, where the acceleration amplitude reaches a = 1.3 -1.5 g, the interior part of the invert has also been damaged, with a pattern similar to the one observed at the end of the Kobe excitation. Tensile cracks appear in the interior of the left and right tunnel shoulders by t = 11.5 s, while the exterior part of the crown also displays extensive damage after the last great acceleration pulse (a = 1.6 g) at t = 13 s.
In qualitative agreement with the empirical observations of Power et al. (1998), a very strong intensity motion is required to observe M. Antoniou, et al. Tunnelling and Underground Space Technology 101 (2020) 103353 extensive damage in a bored tunnel in rock. However, the fact that cracks develop on the tunnel exterior and, most importantly, that is where they start from, is an observation of particular interest. In a real tunnel under earthquake loading (even the low intensity one), these cracks would not be detectable during post-earthquake inspections, rendering a simple visual tunnel assessment problematic. A more sophisticated monitoring system would be needed, or an alternative means to indicate the existence of such damage. Fig. 18 displays results from white noise shaking in Model A prior and after the Wolong Station record. Degradation of the system's first natural frequency from = f m 1, 21.9 Hz to = f ' m 1, 21.0 Hz is vividly illustrated in the transfer functions plots, and is explicitly attributed to the tensile damage developed at the lining during Wolong station excitation. The use of Mohr-Coulomb constitutive law for the soil (which does not account for the degradation of shear stiffness G due to increase in shear strain γ during dynamic loading) allows us to isolate the change of the system's dynamic properties solely due to the tensile damage on the lining. Future studies should also focus on including the soil non-linear behaviour more exhaustively. The soil model used, alongside the reinforcement influence, also explains the small decrease we observe after such a strong seismic event, compared to the results of Wang et al., (2015) who reported a significant drop in the system's first natural frequency after seismic motions that exceed PGA = 0.5 g. However, this indicates that acceleration monitoring in tunnels accompanied by an appropriate analysis, such as the one of Fig. 18, can help to identify changes in the dynamic properties of a tunnel and offer indirect warning for structural degradation.

Summary & conclusions
The paper summarizes and tests a new set of elastic and postcracking scaling laws for 1 g shaking table tests, simulating the earthquake performance of lightly reinforced concrete tunnels founded in rock. Adopted from Corrado et al. (2011), the proposed post-cracking similitude relations are validated through numerical analyses using Abaqus. The applicability of similitude laws for a range of plausible boundary conditions used in typical shaking table tests is shown for three cases: Models A and B correspond to experiments conducted in a rigid and a laminar box respectively, while Model C refers to a very wide model, aiming to exclude the presence of boundary effects.
The nonlinear concrete response is simulated with the Concrete Damaged Plasticity model. The results confirm a good similitude between model and prototype behavior for both medium and high intensity seismic motions. The calculated mean tensile damage of the tunnel lining at the end of the records confirms the previous observation for all models A, B and C. The proposed similitude laws open the doors for the development of novel model materials for the tunnel lining to study the post-cracking behavior of tunnels.
For the examined tunnel structure, the observed crack patterns on the lining indicate that the initiation of cracking systematically occurs at the exterior side of the invertsidewalls connection (mainly due to the high initial stresses at this point prior to seismic loading), revealing that the most vulnerable part of the tunnel section is a non-accessible point for visual inspections. This means that alternative means to assess the state of a tunnel lining post-earthquake is necessary to capture changes due to cracking.
The significant role of boundary effects in the encountered soil-structure interaction problem is highlighted and gives guidance on future shaking table tests. Models A and B, examined and compared here as possible alternatives to the, impractical very wide benchmark Model C, exhibit clear boundary effects. However, it is shown that in cases where the model width is a prescribed, or constrained, parameter, boundary effects may be limited if appropriate ratios of soil stiffness to stiffness of the lining material are used.

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. Fig. 18. Degradation of the system's natural frequency after the Wolong station excitation, explicitly due to the tensile damage developed at the RC lining. Results of white noise shaking before and after the EQ are presented in small scale for Model A. The ratio S of power response spectra between the tunnel invert (Location A) and the model base (Location I) is considered.