# System Identification of Tunnel Response to Ground Motion Considering a Simplified Model

^{1}School of Civil Engineering, University of Leeds, Leeds, United Kingdom^{2}School of Science and Engineering, University of Dundee, Dundee, United Kingdom^{3}Department of Engineering Science, University of Greenwich, London, United Kingdom^{4}Chair of Geotechnical Engineering, ETH Zürich, Zurich, Switzerland

This paper presents a simplified model using coupled Bernoulli and shear beams supported by Winkler-type springs for system identification of tunnels subjected to earthquake-induced ground motions. The model formulation is introduced and closed-form solutions to the modal characteristic equation and mode shapes are derived. The latter are verified against the results of numerical models. Then, the system identification algorithm is presented, showcasing its ability to recover model parameters when recorded acceleration time histories along the tunnel are compromised by noise and sensor locations are varied. The presented framework can be used to initially and, simply, recover tunnel response when monitoring data is available, or to plan for monitoring campaigns in newly constructed or existing tunnels.

## Introduction

Underground structures are less susceptible to earthquake-induced damage as displacements and total accelerations increase per the free surface (Towhata, 2008). However, due to their significant construction cost, tunnels may become bottlenecks in transportation routes in case their functionality is even slightly impaired by these events, leading to large downtime, direct and indirect losses. This motivates careful consideration of the effects of seismic motions for design, construction and management of tunnels.

Deployment of structural health monitoring systems (SHMS) can be beneficial in terms of enhancing the seismic resilience of tunnels for the following reasons (Ogie et al., 2017): (a) they allow for assessment of accelerations, displacements and strains in tunnel linings in real time, ensuring that the desired performance levels under serviceability conditions and in the aftermath of extreme seismic events are met; (b) they facilitate implementation of adaptive management by learning from the response observed in previous occurrences of disruptive seismic excitations; and (c) they can help to take timely action to protect roadway users from life threatening conditions (e.g., preventive closure).

However, deployment of SHMS in tunnels has not been as widespread as for other infrastructure assets, such as bridges and buildings (Fraser et al., 2010; Dolce et al., 2017), mainly due to the lack of a standardized approach resulting from the special nature of each project and the harsh environment faced by sensors (Bhalla et al., 2005). Fortunately, such issues have recently been addressed (Bennett et al., 2010; Tondini et al., 2015; Li et al., 2016) resulting to an increase in their deployment (Dikmen, 2016; Wang et al., 2017).

There are two basic approaches for assessing the output of SHMS with regard to system identification: non-parametric and parametric. In the first case, the fundamental behavior is evaluated by involving solely collected data (Celebi, 2007; Mikael et al., 2013; Saaed Alqado et al., 2015; Loh and Chen, 2017). This can be seen as a “black-box” approach, where the supporting framework does not allow clear explanation of the observed results. In parametric system identification, a model of the structure is formulated and its parameters are adjusted until the measured response is replicated. Detailed numerical models (Hatzigeorgiou and Beskos, 2010; Wu and Gao, 2013; Yu et al., 2016a) can be used for this purpose, but convergence to finding a proper solution can be time consuming (if not impossible) due to the large number of parameters involved.

Another approach to parametric system identification is the use of simplified models. For tunnels, such models typically consist of beams on independent Winkler-type springs (Winkler, 1867), involving both numerical (e.g., Anastasopoulos et al., 2007; Yu et al., 2017; Oh and Moon, 2018) and analytical schemes (St. John and Zarah, 1987; Karadeniz, 2001; Sanchez-Merino et al., 2009; Yu et al., 2016b). However, the seismic response of continua (i.e., soil) can be modeled better by including a shear beam over the Winkler base, allowing for interaction among individual springs (Girija-Vallabhan and Das, 1991), rather than using a single layer of independent springs (Nogami and Lam, 1987). In geotechnical earthquake engineering, a similar approach has been employed for assessing the seismic response of piled foundations (Valsangkar and Pradhanang, 1988) and retaining walls (Younan et al., 1997).

This paper develops a tunnel system identification method that uses recorded acceleration time histories. The analysis is performed by modal superposition of the response of a simplified model comprising coupled flexural and shear beams supported by Winkler-type springs. The inspiration comes from Osawa (1965) and Miranda (1999), who successfully developed a similar approach for the assessment of the seismic response of buildings. Its scope has recently been widened to system identification of large sets of buildings (more than 60 cases) in California (Alimoradi et al., 2006).

## Model Description and Assumptions

As previously mentioned, a simplified model is developed, comprising Winkler-springs supporting a Bernoulli and a shear beam coupled with axial links, in such way that they are subjected to the same transverse displacements in the *y* direction (vertical) as shown in Figure 1. The model is based on the following assumptions and simplifications:

a) The response is purely elastic, allowing for superposition;

b) No axial loads along the tunnel axis are considered, the only loading stems from the ground motion;

c) The far-field ground motion along the subgrade is the same along the tunnel length;

d) There is no interaction between axial and transverse displacements, or among out of plane displacements, and therefore the response in this direction can be superimposed directly;

e) The effects of rotational inertia are not accounted for;

f) The deformation of the cross section of the flexural beam follows Bernoulli's hypothesis; and

g) The tunnel and subgrade materials are uniform.

The equation of motion describing the vibration of the model under the above assumptions can be described as follows (Yu et al., 2016b):

where: *z* is the ordinate along the tunnel length, *t* is time, *u* is the vertical deflection of the tunnel, *EI* is the flexural stiffness of the Bernoulli beam, *AG* is the shear stiffness of the shear beam, *k* is the Winkler subgrade reaction modulus, ρ is the unitary mass of the tunnel per length, and *u*_{g} denotes the far–field uniform vertical ground motion under the foundation.

Equation (1) can be solved through modal analysis. Firstly, its free vibration is governed by the following differential equations:

where: *q* determines the temporal nature of the response, φ is the mode shape and ψ [1/s] is the coupling parameter between Equations (2) and (3) that will define the system response; and take a set of discrete infinite values, each one related to a different mode shape. Furthermore, the variable *z* has been normalized by the tunnel length *H* (*x* = *z/H*) in Equation (3), introducing the following parameters:

where: α is the non-dimensional ratio of shear to flexural stiffness, effectively controlling the role of coupled shear deformation among the supporting springs (Pasternak, 1954; Miranda and Taghavi, 2005), and τ is the tunnel mass parameter normalized by length and flexural stiffness, allowing for definition of a non-dimensional frequency parameter, as ωτ.

From Equation (2) the modal, non-dimensional, circular vibration frequency of the *i*-th mode is:

where *k*_{b} is the non-dimensional ratio between subgrade and flexural stiffness defined below,

The general solution of Equation (3) is given by the following expression:

where: *A*_{1}*, A*_{2}*, A*_{3}*, A*_{4} are coefficients that depend on the end boundary conditions. The variables γ_{i} and β_{i}, which describe the *i*-th mode shape, are related to model parameters as follows:

Consequently, if γ_{i} and β_{i} are found, they will define the non-dimensional product ψ_{i}τ according to Equation (9). With this in hand, it is possible to use Equation (5) to calculate the normalized, non-dimensional, modal circular vibration frequencies. Furthermore, if modal vibration periods are scaled by the first fundamental period, *T*_{1}, the following is obtained:

The above equations fully define the simplified model. Its response to ground motion can be obtained by setting these four parameters: α, *k*_{b}*, T*_{1} and the fraction of critical damping for each mode ξ_{i}.

Fixed–fixed end boundary conditions at the ends are considered, being representative of a tunnel spanning two close-by underground stations. Under such fixed supports at both beam extremities, the following boundary conditions apply for each mode:

Applying the above boundary conditions to Equation (7) defines the following system of equations in matrix form:

In order to avoid a null solution to the above system, the determinant in Equation (12) must be zero, leading to the following modal characteristic equation:

The roots of Equation (13) define the β_{i} and γ_{i} parameters needed for evaluating the period ratios in Equation (10), and the following mode shapes:

where η_{i} is:

It must be noted that the mode shapes are independent of the value of the subgrade reaction modulus, *k*_{b}, as they only depend on γ_{i} and β_{i}. This is valid because a uniform support is used as a simplification.

## Response to Ground Motion and Modal Properties Verification

Equation (1) can be solved through modal decomposition, and so the displacement is given by (Paz, 2006):

Equation (16) constitutes an exact solution only if an infinite number of modes *N* is considered, but usually up to 5 modes are sufficient (Miranda and Taghavi, 2005). The values of *r*_{i} can be found after solving the following ordinary differential equation:

where *Ug*_{i} is the equivalent ground motion for each mode:

However, evaluation of Equation (18) requires knowledge of the ground motion field at all times and positions under the Winkler base. This information cannot be expected to be available under normal circumstances, as it would require deployment of dense sensor networks within the support medium around the tunnel, and would still not provide the exact *u*_{g}(*x, t*) function. Therefore, a more practical and realistic source of excitation data is out-of-site strong motion recording stations. Hence, a more convenient approach is to assume uniform base excitation, in which case Equation (17) becomes:

where parameter Γ_{i} is given by:

Equation (19) is the same with the one that governs the equation of motion of a Single Degree of Freedom (SDOF) oscillator, except that it is scaled by Γ_{i}. For the stated model hypotheses, closed form expressions for this parameter are given in Equation (A.1) in Appendix A. In this study, the algorithm developed by Nigam and Jennings (1969) is employed to solve Equation (19) for each mode. Total acceleration is found by taking the second derivative of Equation (16), considering Equation (20) and adding the base ground motion:

### Verification of Modal Parameters

The developed simplified model is verified against the results of a dense finite element (FE) model, which is used as a benchmark. Its correctness will be assessed in the future against field measurements and laboratory tests. The closed–form expressions proposed for both modal period ratios and mode shapes considering the fixed-fixed boundary conditions are used for the comparison.

The FE model is implemented in SAP 2000 (Wilson, 2002). Two beams, comprising 1,000 beam elements each (with a unitary length) are placed next to each other (separated by an unitary distance) and they are subjected to equal vertical displacements. For one of the two beams, the rotation is also constrained (set to zero), while its flexural stiffness is adjusted to represent the shear stiffness associated with the α value being considered. The flexural stiffness of the other beam is set to a unitary value, while its nodes are coupled with springs and unitary (translational) masses. The stiffness of the Winkler springs is adjusted to represent the *k*_{b} values.

For buildings, values of α range between 0 and 30 (Alimoradi et al., 2006). As tunnels usually have large length-to-diameter ratios, they are usually idealized as thin-walled cylinders (Kouretzis et al., 2006, 2011). It is expected that their deformation patterns resemble more what is observed in flexural (Bernoulli) beams, rather than shear beams. Consequently, as a first approximation, values between 0 and 20 were considered in this study. Previous research regarding the dynamic behavior of beams on Pasternak foundations considered values up to α = 3 (Valsangkar and Pradhanang, 1988; Filipich and Rosales, 2002).

The parameter *k*_{b} takes values between 4 × 10^{−3} *H*^{2} and 3 × 10^{−2} *H*^{2}. The derivation of these values is presented in Appendix B and was achieved by considering the tunnel design recommendations of the Japanese Railway Research Institute, which are based on data inferred from plate tests (Zhang et al., 2014), while adopting a tunnel cross-section similar to the one of the T3 Santiago de Chile Metro Line, which is idealized as a hollow circle of 10 m outer diameter and 0.5 m wall thickness, and concrete modulus of elasticity of 30 GPa. Furthermore, tunnel lengths (*H*) between 100 and 250 m were considered, which constitute a reasonable spacing between successive metro stations in dense urban areas.

The above values result in *k*_{b} limit values of 40 and 1,000. Beyond this upper value, tunnel flexibility is so large that its displacement follows the ground motion, rendering pseudo-static analysis valid (Valsangkar and Pradhanang, 1988)—less interesting in the context of this paper. Hence, based on the above, the following cases are considered: α = 0 and *k*_{b} = 5; α = 0 and *k*_{b} = 25; α = 7 and *k*_{b} = 1; α = 7 and *k*_{b} = 1,000; α = 20 and *k*_{b} = 10; α = 20 and *k*_{b} = 100, representing extreme values. As shown in Figure 2, for all cases examined the comparison proves the correctness of analytical solution, with overall differences between the closed-form solution for period ratios, mode shapes and the benchmark numerical simulations being < 0.1%.

## System Identification Algorithm and Results

The simplified model, defined by four parameters, can be used to undertake system identification under the assumption that sensor networks are available to measure actual tunnel response. Setting proper values for parameters describing the model can be done by involving system identification procedures on both frequency (McVerry, 1980) and time domain (Beck and Jennings, 1980) following the outline described in Figure 3.

The procedure starts by defining an initial set of parameters, which is iteratively updated until the difference between modeled and measured response at the desired locations becomes less than a threshold. In principle, gradient dependent methods (Arora, 2015) could be employed. However, they are prone to yield local instead of global optimal values. An alternative is to use heuristic methods (Bozorg-Haddad et al., 2017), which mimic biological and physical phenomena to systematically search a dominion of plausible values, allowing for aleatory search beyond regions where local minima are observed. It must be stressed that heuristic algorithms cannot guarantee that an absolute minimum is achieved; consequently, it is common practice that the final results based on multiple runs are critically reviewed. In this paper, the square difference between measured and predicted total acceleration time histories is set as the objective function (*OF*) to be minimized:

where: *n* is the number of total acceleration time histories recorded (ü_{TM}) at locations *x*_{j}.

### Genetic Algorithms

Among Heuristic methods, Genetic Algorithms (GA) and GA-based procedures have been employed widely to solve problems in geomechanics, including the estimation of lateral load capacity of piles, undrained side resistance of drilled shafts, tunneling-induced ground motion settlements, and soil liquefaction (Gandomi and Alavi, 2012).

GA are programming schemes for searching a domain space involving a process that mimics natural selection of biological species (Affenzeller et al., 2009); firstly, initial values (“first generation”) of parameters are generated randomly. Each full set of parameters (α, *k*_{b}, *T*_{1}, and ξ) constitutes an individual. Once the first generation is produced, the *OF* for all individuals is calculated and the best performing become “parents”. These parents are then combined through a process called crossover to define new individuals, called “children”. In this study, crossover is done through roulette selection, in such way that the chances of being selected for crossover are proportional to the inverse value of the *OF*, i.e., the lower the *OF* it is, the higher the possibility of being combined with another high-valued individual.

Crossover is made through a weighted average of parameters that comprise each parent. Weights are defined according to the value of *OF* for each one of them. Therefore, parameters in children will resemble mostly the parent with the lowest *OF* value. Also, the algorithm introduces aleatory modifications randomly in children parameters, following a prescribed probability of occurrence. This process is called “mutation”. Finally, all samples are ranked. The best suited individuals are included in the next generation, while the worst performing ones are discarded. The process continues until no meaningful improvement in minimizing the *OF* is found. In this study, the algorithm is implemented using Matlab Optimization Toolbox function GA (The Mathworks, 2016), with the following parameters: population size of 50 individuals, 400 maximum allowable generations, a stop criterion if changes in the *OF* are < 10^{−6} and 5% chance of mutation.

### Investigated Cases

Ten discrete samples were selected from a wide range of possible values for each one of the four parameters describing the model. A range of indicative values is considered for each variable: α ranges from 1 to 19; *k*_{b} between 1 and 400; *T*_{1} between 0.3 and 1.2 s (Wang et al., 2015); and ξ (which was considered uniform for all modes) between 0.5 and 5%. Then, 10 random combinations of these parameters were created, leading to the test cases outlined in Table 1.

The vertical component of two ground motions were selected to excite the model at its base: the Hulverstone Drive Pumping Station (HDPS) record, obtained within the city of Christchurch, during its 2011 M_{w} 6.2 earthquake, and the CHY080 record obtained during the 1999 M_{w} 7.6 Chi-Chi Taiwan earthquake mainshock. Both records are characterized by high peak ground acceleration (PGA) values in the vertical direction, i.e., 10.5 m/s^{2} for HDPS and 7.50 m/s^{2} for CHY80. The different acceleration spectral response of each motion is presented in Figure 4 (Pacific Earthquake Engineering Research Center, 2018).

**Figure 4**. Ground motions considered in this study. Time histories on left, acceleration response spectra (ξ = 2%) on right.

The model response was recorded at 9 locations, considering 3 configurations: uniformly spaced along the full length; clustered in the middle between *x* = 0.3 and *x* = 0.7 considering spacing half the one employed in the uniform case; and deployed just in the first half (*x* ≤ 0.5) again with half the spacing of the uniform layout. Additionally, the records of all sensors were compromised by adding independent white Gaussian noise with a variance equal to the one observed for the ground motion. Finally, a reduced network of only 5 sensors, deployed at *x* = 0.2, 0.4, 0.5, 0.6, 0.8 was assessed fora subset made by the best performing tests cases, i.e., 1, 6 and 7 (see Table 1 for naming convention).

### Results

As summarized in Table 2, for values of *T*_{1} < 0.5 s the algorithm is unable to recover a useful solution. Beyond this limit, the algorithm successfully recovers model parameters, with divergences < 5% in more than 75% of all runs as long the fundamental period is < 1 s and in more than 50% for values beyond that threshold. Among successfully obtained parameters, differences between estimates and target values are < 1% in more than 90% of successful runs, and in all cases < 2%. Indicative target against estimated responses are shown in Figure 5.

Case 2 is particularly interesting because the α value is not well estimated, ranging between 0 and 1.3. This is because for values of α below 3.0 the influence of this parameter in the acceleration response is small: this was observed by the small variations in the mode shapes and period ratios for α < 3 as shown in Figure 6. In other words, for this range of values, considering the Bernoulli Beam on a spring base is a reasonable approximation, as changes in a value for this range produce slight changes on acceleration response. Beyond 3.0, the effects of α become more pronounced, noticeably influencing acceleration response and therefore, allowing the algorithm to capture them more accurately.

Effects of changes in sensor arrangement on finally predicted results are not significant. Relevant predicted parameters can be recovered from any of the assessed sensor layouts. Also, results for the reduced array of just 5 sensors were proven satisfactory (Table 3). Furthermore, recovered parameters match closely the target values, showing results that are comparable to what was observed for the denser 10 sensor arrangement.

## Conclusions

The use of a simplified model for system identification of tunnels considering earthquake–induced ground motion is presented in this paper. The model includes coupled flexural and shear beams, supported by Winkler-type springs, in such way that they are subjected to the same vertical accelerations. The proposed model effectively extends a previous model, which was developed for buildings. Closed form solutions for modal assessment of response were obtained in terms of trigonometric and hyperbolic functions. These solutions were verified against the results of a dense FE model, used as a benchmark.

It was shown that the proposed simplified model can be properly defined by 4 parameters only; the first mode period *T*_{1}; the ratio α of shear to flexural stiffness; the ratio *k*_{b} between subgrade stiffness and flexural beam stiffness; and the overall damping ratio ξ for each mode.

Using a simple genetic algorithm and the proposed simplified model, it was shown that a system identification can be performed for variable conditions. The model was tested for two different earthquakes of varied frequency content, and the influence of sensor distribution along the tunnel length was parametrically investigated. The output of tests was intentionally compromised by adding white Gaussian noise with a variance equal to the one observed in the ground motion. It was observed that for α values >3.0 (i.e., when the effect of the shear beam is of significance), model parameters can be successfully recovered. Hence, for such conditions, the proposed approach can be used to recover important dynamic properties of the tunnel (e.g., *T*_{1} or α), but also soil-structure interaction by tracking changes of *k*_{b} following earthquakes. For values < 3, the Bernoulli beam with a Winkler Base is a valid representation, observing differences with a model including a Pasternak foundation to be marginal.

The different sensor distributions and the reduced number of sensors did not lead to significant differences in the final results. This initially shows that sensor position is not the dominating parameter for the case of 5 sensors along a tunnel length fixed at both ends, and the uniform ground conditions and other assumptions made in the paper. Further tests should be conducted on real data and sensors to study the influence of other aspects such as sensor noise.

## List of Variables

*A*: coefficients of the general modal solution.

*AG*: shear stiffness. [N]

*EI*: flexural Stiffness. [Nm^{2}]

*H*: tunnel length. [m]

*k*: Winkler grade reaction modulus. [N/m/m]

*k*_{b}: non-dimensional, normalized Winkler reaction coefficient by flexural stiffness.

*i*: mode index.

*OF*: Objective Function.

*q(t)*: time variation of modal ordinate during free vibration.

*r*: time variation of modal ordinates when subjected to ground motion.

*T*: mode vibration period.

*u*: vertical deflection of the tunnel. [m]

*u*_{g}: displacement at the base of the Bernoulli beam. [m]

ü_{T}: total tunnel acceleration.

ü_{TM}: measured tunnel acceleration.

*U*_{g}: generalized displacement at the base. [m]

*x*: normalized length ordinate by tunnel length.

*x*_{j}: normalized location of j-th sensor.

*z*: ordinate along tunnel length. [m]

α: non-dimensional ratio between shear and flexural stiffness.

β_{i}: parameter of the i-th mode shape.

γ_{i}: parameter of the i-th mode shape.

Γ_{i}: i-th mode Participation Factor.

η_{i}: parameter of the i-th mode shape.

φ_{i}*(x)*, i-th mode shape.

ψ: coupling parameter between Equation (2) and Equation (3).

ξ_{i}: damping coefficient to the critical (for mode i).

ρ: mass per unit length. [kg/m]

τ: normalized mass parameter normalized by flexural stiffness. [1/s]

ω: circular vibration frequency of mode i. [1/s].

## Author Contributions

AA-R: Development of the model and implementation of the Test Algorithm; RF: Development of the model and revision of paper draft; NN: Development of tests cases, definition of the Test Algorithm; JK and GK: Definition of Range of test parameters, Review of Paper and Validation of the Model; IA: Definition of the System Identification Algorithm, Final review of the document.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

This work was sponsored by a Newton Fund Initiative, supported jointly by the Engineering and Physical Sciences Research Council of the UK and the National Commission for Scientific and Technological Research of Chile, grant EP/N03435X/1. Strong motion records employed in this work were obtained from the PEER NGA 2 database, hosted by the Pacific Earthquake Engineering Research Center.

## References

Affenzeller, M., Winkler, S., Wagner, S., and Beham, A. (2009). *Genetic Algorithms and Genetic Programming*. Boca Raton, FL: CRC Press.

Alimoradi, A., Miranda, E., Taghavi, S., and Naeim, F. (2006). Evolutionary modal identification utilizing coupled shear-flexural response-implications for multistory buildings. *Struct. Des. Tall Spec. Build*. 15, 51–65. doi: 10.1002/tal.343

Anastasopoulos, I., Gerolymos, N., Drosos, V., Kourkoulis, R., Georgarakos, T., and Gazetas, G. (2007). Nonlinear response of deep immersed tunnel to strong seismic shaking. *J.Geotech. Geoenviviron Eng.* 133, 1067–1090. doi: 10.1061/(ASCE)1090-0241(2007)133:9(1067)

Beck, J. L., and Jennings, P. C. (1980). Structural Identification using linear models and earthquake records. *Earthquake Eng. Struct. Dyn.* 8, 145–160. doi: 10.1002/eqe.4290080205

Bennett, P. J., Kobayashi, Y., Soga, K., and Wright, P. (2010). Wireless sensor network for monitoring transport tunnels. *ICE Geotech. Eng.* 163, 147–156. doi: 10.1680/geng.2010.163.3.147

Bhalla, S., Yang, Y. W., Zhao, J., and Soh, C. K. (2005). Structural health monitoring of underground facilities – technological issues and challenges. *Tunnel. Underground Space Technol.* 20, 487–500. doi: 10.1016/j.tust.2005.03.003

Bozorg-Haddad, O., Solgi, M., and Loaiciga, H. A. (2017). *Meta-heuristic and Evolutionary Algorithms for Engineering Optimization*. Hooboken, NJ: Wiley-Blackwell.

Celebi, M. (2007). “On the variation of fundamental frequency (period) on an undamaged building – a Continuing Discussion,” in *International Conference on Experimental Vibration Analysis for Civil Engineering Structures*. Porto.

Dikmen, S. U. (2016). Response of marmaray submerged tunnel during the 2014 Northern Aegean Earthquake (Mw = 6.9). *Soil Dyn. Earthquake Eng.* 90, 15–31. doi: 10.1016/j.soildyn.2016.08.006

Dolce, M., Nicoletti, M., De-Sortis, A., Marchesini, S., Spiuna, D., and Talanas, F. (2017). Osservatorio sismico delle strutture: the italian structural seismic monitoring network. *Bull. Earthquake Eng.* 15, 621–641. doi: 10.1007/s10518-015-9738-x

Filipich, C., and Rosales, M. B. (2002). A further study about the behaviour of foundation piles and beams in a Winkler-Pasternak Soil. *Int. J. Mech. Sci.* 44, 21–36. doi: 10.1016/S0020-7403(01)00087-X

Fraser, M., Elgamal, A., He, X., and Conte, J. P. (2010). Sensor network for structural health monitoring of a highway bridge. *J. Comput. Civil Eng. ASCE* 24, 11–24. doi: 10.1061/(ASCE)CP.1943-5487.0000005

Gandomi, A., and Alavi, A. (2012). A new multi-gene genetic programming approach to non-linear system modeling. Part II: geotechnical and earthquake engineering problems. *Neural comput. Appl.* 21, 189–201. doi: 10.1007/s00521-011-0735-y

Girija-Vallabhan, C. V., and Das, Y. C. (1991). A refined model for beams on elastic foundations. *Int J. Solids Struct.* 27, 629–637. doi: 10.1016/0020-7683(91)90217-4

Hatzigeorgiou, G. D., and Beskos, D. E. (2010). Soil-Structure interaction effects on seismic inelastic analysis of 3-D Tunnels. *Soil Dyn. Earthquake Eng.* 30, 851–861. doi: 10.1016/j.soildyn.2010.03.010

Karadeniz, H. (2001). Earthquake analysis of buried structures and pipelines based on Rayleigh wave propagation. *Int. J. offshore Polar Eng*. 11, 133–140. Available online at: http://legacy.isope.org/publications/journals/journalJune01.htm

Kouretzis, G. P., Bouckovalas, G. D., and Gantes, C. J. (2006). 3-D shell analysis of cylindrical underground structures under seismic shear (S) wave action. *Soil Dyn. Earthquake Eng.* 26, 909–921. doi: 10.1016/j.soildyn.2006.02.002

Kouretzis, G. P., Bouckvalas, G. D., and Karamitros, D. K. (2011). Seismic verification of long cylindrical Underground structures considering Rayleigh wave effects. *Tunel. Underground Space Technol*. 26, 789–794. doi: 10.1016/j.tust.2011.05.001

Li, H., Ren, L., Jia, Z., Yi, T., and Li, D. (2016). State-of-the-art in structural health monitoring of large and complex civil infrastructures. *J. Civil Struct. Health Monit*. 6, 3–16. doi: 10.1007/s13349-015-0108-9

Loh, C., and Chen, J. (2017). Tracking modal parameters from building seismic response data using recursive subspace identification algorithm. *Earthquake Eng. Struct. Dyn.* 46, 20163–22183. doi: 10.1002/eqe.2900

McVerry, G. H. (1980). Structural identification in the frequency domain from earthquake records. *Earthquake Eng. Struct. Dyn.* 8, 161–180. doi: 10.1002/eqe.4290080206

Mikael, A., Gueguen, P., Bard, P., Roux, P., and Langlais, M. (2013). The analysis of long-term frequency and damping wandering in buildings using the random decrement technique. *Bull. Seismol. Soc. Am.* 113, 236–246. doi: 10.1785/0120120048

Miranda, E. (1999). Approximate seismic lateral deformation demands in multistory buildings. *J. Struct. Eng*. 125, 417–425. doi: 10.1061/(ASCE)0733-9445(1999)125:4(417)

Miranda, E., and Taghavi, S. (2005). Approximate floor acceleration demands in multi-story buildings I: formulation. *J. Struct. Eng. ASCE* 131, 203–211. doi: 10.1061/(ASCE)0733-9445(2005)131:2(203)

Nigam, N. C., and Jennings, P. C. (1969). Calculation of response spectra from strong-motion earthquake records. *Bull. Seismic. Soc. Am.* 59, 909–922.

Nogami, T., and Lam, Y. C. (1987). Two-parameter layer model for analysis of slab on elastic foundation. *J. Eng. Mech ASCE* 113, 1279–1291. doi: 10.1061/(ASCE)0733-9399(1987)113:9(1279)

Ogie, R. I., Perez, P., and Dignum, V. (2017). Smart infrastructure: an emerging frontier for multidisciplinary research. *Smart Infrastruct. Construct. ICE* 170, 8–16. doi: 10.1680/jsmic.16.00002

Oh, J., and Moon, T. (2018). Seismic design of a single bored tunnel: longitudinal deformations and seismic joints. *Rock Mech. Rock Eng.* 51, 891–910. doi: 10.1007/s00603-018-1434-0.

Osawa, Y. (1965). Response analysis of tall buildings to strong earthquake motions. Part 1: linear response of core-wall buildings. *Bull. Earthquake Res. Inst.* 43, 803–817.

Pacific Earthquake Engineering Research Center, P. E. E. R. (2018). *NGA-West 2 Peer Ground Motion Database*. Available online at: https://ngawest2.berkeley.edu/site (Accessed April 18, 2018).

Pasternak, P. (1954). *On a New Method of Analysis of an Elastic Foundation by Means of Two Foundation Constants.* Moscow:Gosudarstvenrwe Izdatelslvo Literaturi po Stroitclstvu i Arkhitekture (In Russian).

Saaed Alqado, T. E., Nikolakopoulos, G., and Jonasson, J. (2015). Comfort level identification for irregular multi-storey buildings. *Automat. Const.* 50, 40–49. doi: 10.1016/j.autcon.2014.10.009

Sanchez-Merino, A., Fernandez-Saez, J., and Navarro, C. (2009). Simplified longitudinal response of tunnel linings subjected to surface waves. *Soil Dyn. Earthquake Eng.* 29, 579–582. doi: 10.1016/j.soildyn.2008.06.003

St. John, C., and Zarah, T. (1987). Aseismic design of underground structures. *Tunnel. Underground Space Technol.* 2, 165–197. doi: 10.1016/0886-7798(87)90011-3

Tondini, N., Bursi, O., Bonelli, A., and Fassin, M. (2015). Capabilities of the Fiber Bragg Sensor System to monitor the inelastic response of concrete sections in new tunnel linings subjected to earthquake loading. *Comput. Aid. Civil Infrastruct. Eng.* 30, 636–656. doi: 10.1111/mice.12106

Towhata, I. (2008). *Geotechnical Earthquake Engineering*. Berlin, Germany: Springer Science & Business Media.

Valsangkar, A. J., and Pradhanang, R. (1988). Vibrations of beam-columns on two-parameter elastic Foundations. *Earthquake Eng. Struct. Dyn.* 16, 217–225. doi: 10.1002/eqe.4290160205

Wang, B., Zhang, Z., He, C., and Zheng, H. (2017). Implementation of a long-term monitoring approach for the operational safety of highway tunnel structures in a severely seismic area of China. *Struct. Control Health Monitor.* 24, 1–18. doi: 10.1002/stc.1993

Wang, Z.Z., Jiang, Y.-J., Zhu, C.A., and Sun, T.C. (2015). Shaking table tests of tunnel linings in progressive states of damage. *Tunnel. Underground Space Technol.* 50, 109–117. doi: 10.1016/j.tust.2015.07.004

Wilson, E. (2002). *Sap2000 Integrated Software for Structural Analysis and Design, Analysis Reference Manual*. Berkeley, CA: Computers and Structures Inc.

Winkler, E. (1867). *Die Lehre von der Elasticitaet und Festigkeit. Theil, 2*. Prague: H. Dominicus (in German).

Wu, D., and Gao, B. (2013). Dynamic analysis of tunnel entrance under the effect of Rayleigh wave. *Electron. J. Geothech. Eng.* 18, 4231–4246.

Younan, A., Veletsos, A., and Bandyopadhyay, K. (1997). *Dynamic Response of Flexible Retaining Walls.* Report BNL 52519, Department of Advanced Technology, Brookhaven National Laboratory, Upton, NY.

Yu, H., Cai, C., Guan, X., and Yuan, Y. (2016b). Analytical solution for long lined tunnels subjected to travelling loads. *Tunnel. Underground Space Technol.* 58, 209–215. doi: 10.1016/j.tust.2016.05.008

Yu, H., Chen, J., Bobet, A., and Yuan, Y. (2016a). Damage observation and assessment of the Longxi tunnel during the Wenchuan earthquake. *Tunnel. Underground Space Technol.* 54, 102–116. doi: 10.1016/j.tust.2016.02.008

Yu, H., Yuan, Y., and Bobet, A. (2017). Seismic analysis of long tunnels, a review of simplified and unified methods. *Underground Space* 2, 73–87. doi: 10.1016/j.undsp.2017.05.003

Zhang, D., Huang, H., Phoon, K. K., and Hu, Q. A. (2014). Modified solution of radial subgrade modulus for a circular tunnel in elastic ground. *Soils Found.* 54, 225–232. doi: 10.1016/j.sandf.2014.02.012

## Appendices

### Appendix A

The solution to Equation (20) is provided below for completeness.

### Appendix B

Reference values are based on the cross section of the Santiago Metro L3 line. Its section can be idealized as circular, with 9m internal and 10m external diameter. Consequently, the second moment of inertia is:

and the outer, perimeter:

Zhang et al. (2014), based on specifications by the Japanese Railway Research Institute and Chinese standards, outline design values for plate reaction modulus of tunnels (*k*_{p}) between 3 and 150 MPa/m. This would lead to the following thresholds for the base spring coefficient

If concrete with 28 day cylinder strength of 35 to 50 MPa is assumed as material with Young's modulus of 35 GPa, the threshold values for the *k*_{b} parameter can be computed as:

Keywords: system identification, simplified models, tunnel seismic response, closed-form solutions, structural monitoring

Citation: Alonso-Rodriguez A, Nikitas N, Knappett J, Kampas G, Anastasopoulos I and Fuentes R (2018) System Identification of Tunnel Response to Ground Motion Considering a Simplified Model. *Front. Built Environ*. 4:39. doi: 10.3389/fbuil.2018.00039

Received: 01 May 2018; Accepted: 29 June 2018;

Published: 20 July 2018.

Edited by:

Luigi Di Sarno, University of Sannio, ItalyReviewed by:

Michele Palermo, Università degli Studi di Bologna, ItalyAristotelis E. Charalampakis, National Technical University of Athens, Greece

Copyright © 2018 Alonso-Rodriguez, Nikitas, Knappett, Kampas, Anastasopoulos and Fuentes. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Raul Fuentes, r.fuentes@leeds.ac.uk