Response of a high-speed train travelling over a long and high-pier viaduct during moderate earthquakes

This contribution presents the results of numerical simulations in which a high-pier viaduct under lateral and vertical non-stationary spatially variable seismic ground motions is crossed by an articulated high-speed train at several speeds. The spectral representation method is employed to generate the required number of acceleration time-history series at each support location of the railway viaduct studied. Wave passage, local soil conditions and loss of coherency effects have been accounted for in the simulated ground-motion time histories, with several peak ground accelerations. Train and bridge responses have been obtained by means of a non-linear dynamic interaction multibody and finite element model. An analysis of train running safety indices has been carried out. The properties of earthquakes and train speeds which may cause problems to traffic safety are detailed, as is the train derailment probability. The numerical simulation outcomes show that it is not safe for the train to travel over the studied viaduct in earthquakes with peak ground acceleration equal to 0.03 g and for train speeds above 280 km/h.


introduction
Strict high-speed lines layout parameters and increasing environmental requirements are known to be leading to the design and building of longer and higher railway bridges, particularly in mountainous areas. Moreover, the frequency of trains' passage over this kind of viaducts is increasing due to the increase in the number of high-speed train passengers. For these reasons, the probability of a high-speed train passing over a high bridge when an earthquake occurs is also rising. The ride comfort, and even more importantly, the traffic safety during seismic ground motions at these bridges, have not been studied thoroughly in the scientific literature.
The behavior of bridges accounting for spatial variability of ground motions has been studied in several works (Shinozuka et al. 2000;Sextos et al. 2003;Lupoi et al. 2005;Zerva 2009). Shinozuka et al. (2000) studied the effects of the spatial variability of ground motion on the seismic response of highway reinforced concrete (RC) bridges to provide guidelines for the seismic design of these structures. For that study they performed inelastic time history dynamic analyses. The results of the analyses carried out showed that the peak of ductility demand at the columns could increase substantially when asynchronous support ground motions were used and different local soil conditions were considered at each bridge support. Therefore, the identical support ground motion assumption seems to be unconservative in concrete bridge design. Sextos et al. (2003) found that bridges subjected to spatially variable input motions excited bridge higher vibration modes which had antisymmetric shapes. The site effects, i.e., the different soil conditions on the bridge supports, also increased the bridge response so those authors introduced different peak ground accelerations and different spectral amplifications. Sextos concluded that in bridges longer than 400 m, the effects of the spatial variation of the seismic action should be considered. Lupoi et al. (2005) also studied the influence of spatial variation of seismic action on the seismic design of RC bridges. Those authors stated that the increase in the support ductility demand when considering asynchronous seismic ground motions depended mainly on the bridge stiffness. Zerva (2009) explained in her book, based on a complete literature review, the influence of the different features of spatially variable seismic ground motions on the response of RC bridges. These features were the wave passage effect, the loss of coherency and the site effects. Zerva affirmed that when the loss of coherency in the seismic motions increased, it overshadowed the weave passage effect. On the other hand, the effects of loss of coherency and the passage of the seismic wave were overshadowed by the effect of site conditions at the bridge supports.
The four works cited above, which studied the effects of the spatial variation of seismic action on bridges, did not consider moderate earthquakes, nor their effect on highpier railway viaducts.
The train response when it is crossing a bridge during an earthquake has also been an objective of some studies where non-linear wheel-contact models were used in order to properly describe the derailment phenomenon (Ju 2013;Montenegro et al. 2016;Zeng and Dimitrakopoulos 2018;Tanabe et al. 2016). Ju (2013) carried out a parametric study in which railway traffic safety on standard multi-span bridges was analyzed. Ju concluded that stiffening the piers of those bridges improved the traffic safety during earthquakes. Therefore, high-pier railway viaducts, such as the one studied in the present paper, could cause greater traffic safety problems under the action of earthquakes. Montenegro et al. (2016) studied the influence of the earthquake intensity, the train running speed, and the track quality on traffic safety. The bridge studied was of several simply supported spans 21 m in length supported by columns between 5 and 15 m high. From the results shown by Montenegro et al. in their study it can be stated that a train traveling at 350 km/h does not travel safely over the bridge studied when an earthquake of Peak Ground Accelerations (PGA) equal to 0.09 g occurs. Zeng and Dimitrakopoulos (2018) studied the behavior of a train traveling on a bridge under the action of different earthquakes using a train-bridge interaction model that included a sophisticated wheel-rail contact model, but did not consider the height of the bridge piers. Tanabe et al. (2016) focused their work on developing the train-wheel-rail-bridge interaction model for earthquakes but did not carry out a parametric study of traffic safety during earthquakes on railway viaducts.
Although in these last four works (Ju 2013;Montenegro et al. 2016;Zeng and Dimitrakopoulos 2018;Tanabe et al. 2016) the authors used non-linear wheel-contact models to suitably obtain wheel-rail interaction forces and therefore, traffic safety indices, they studied short bridges and did not consider the special variation of seismic action.
On the other hand, some contributions have studied the train response over the bridge taking into account the spatial variation of the seismic ground motion, albeit using linear models and the traffic safety was not suitably attended (Zhang et al. 2011;Jia et al. 2013).
It is difficult to find papers about railway viaducts where the running safety of a highspeed train is studied considering the following two issues: non-linear wheel-rail contact models, and the spatially variable ground motions as input excitation. One such study was carried out by Xia et al. (2006). The bridge studied by Xia et al. had seven spans of continuous and variable cross-section prestressed concrete box girders, 48 + 5 × 80 + 48 m. Its piers were up to 22 m high. The train modelled in Xia's study was a China-Star. That study showed that the travel of this train travelling at a speed higher than 300 km/h would not be safe when a transverse earthquake with PGA equal to 0.06 g occurred. Xia et al. (2006) studied the safety of railway traffic on a bridge during earthquakes with appropriate train-bridge interaction models and considering the effects of the spatial variation of seismic action, obtaining interesting results. However, we can wonder, what is the risk to traffic safety on higher viaducts and with lighter trains than that used by Xia? The present study seeks to answer that question.
The non-linear train-bridge interaction model used herein has been developed in previous work by the authors (Olmos and Astiz 2018a, b). This interaction model consists of a multibody model for the train, a finite element model for the bridge, and a non-linear wheel-rail contact model. The non-linearities of this last model are due to the geometry of the wheel and rail profiles and the friction at the contact point. An analysis of train running safety indices has been carried out. In the present contribution, running safety assessment of a high-speed train travelling over a long and high-pier viaduct during moderate earthquakes is studied. The lateral flexibility of these high-pier viaducts may significantly increase the pseudo-static bridge response which is due to the non-uniform seismic ground motion. In this way, the risk of vehicle derailment on this kind of bridges during an earthquake could be higher than in other bridges with short piers. Two important factors that may affect train running safety, namely ground motion intensity and train speed, are considered here in a parametric analysis.
The bridge studied is an idealization of the Arroyo de las Piedras viaduct (Mato et al. 2007), built in Spain in 2007, which is a composite steel-concrete continuous deck bridge. The distribution of span of the bridge is 19 spans which have 50.4 + 17 × 63.5 + 44 m (Fig. 1). Its total length is 1173.9 m. Several of its hollow-section concrete piers are over 93 m high. The links between the deck and the different bridge supports (piers and abutments) consist of pot-type bearings. At least one of the two deck-bearing devices of each support is fixed in the bridge lateral direction. In the bridge longitudinal direction, the bearing devices are sliding in all the bridge supports except in the four central piers (pier 8, 9, 10, and 11), where they are fixed. Although the longitudinal seismic behavior of the bridge is not studied in this work, in order to complete the description of the deck support conditions it can be said that this deck also has two impact transmitters fitted on each abutment, incorporating damping devices against the seismic longitudinal force.
The articulated train type Thalys (or TGV-Alstom) of 10 vehicles has been chosen for this study because it is a light train.
The numerical simulation outcomes show that the train travels over the bridge are not safe for earthquakes with PGA equal to 0.03 g and for train speeds higher than 280 km/h. However, the occurrence probability of this earthquake at this bridge location when a train is passing is lower than the occurrence probability of the design earthquake, considering a train passage frequency of one every hour, from 7 a.m. to 11 p.m. However, in areas of greater seismicity and greater train passage frequency, this probability increases notably. This paper is organized as follows. The equation of motion of coupled train and bridge systems subject to earthquakes is shown and discussed in Sect. 2. Sections 3, 4 and 5 describe the bridge, train, and wheelset-track interaction models, respectively. Section 6 shows the procedure of Deodatis (1996) followed to numerically generate the non-uniform seismic ground motions that have been used as inputs in the traffic safety study performed. In Sect. 7, the implementation and numerical solution of the train-bridge dynamic interaction model used in this study are briefly described. Section 8 shows the outcomes of the railway traffic safety analysis carried out and finally, the conclusions are set out in Sect. 9.

Equations of motion for train-bridge systems subjected to earthquakes
When the finite element model (FEM) and the multibody model are used for describing the dynamic behavior of train and bridge systems subjected to non-uniform seismic ground motions at bridge supports, the dynamic motion equations of the train-bridge system can be expressed as Eqs. (1) and (2). In these equations, the absolute global coordinate system is adopted. v and b subscripts denote train subsystem and superstructure bridge subsystem, respectively, and s subscript denotes bridge supports. vb and bv are non-linear vehicle-bridge interaction forces acting on the train and bridge subsystems respectively. ̈ k (t),̇ k (t) , and k (t) are the acceleration, velocity, and displacement vectors according to the degrees of freedom considered in the train model (k = v) and in the bridge model (k = b) conveniently arranged. k , k , and k are the mass, damping, and stiffness constant matrices of the dynamic vehicle system (k = v), bridge superstructure system (k = b), or bridge supports (k = s). Note that ̈ s (t),̇ s (t) , and s (t) are vectors which contain the seismic ground vertical and horizontal movements at the bridge supports during an earthquake.
When Eq.
(2) is solved using the direct solution method (Xia et al. 2018), the first term on the left side of the equation is expanded, and the motion equations of the bridge in the absolute coordinates can be written as Eq. (3): As in traditional seismic analysis, the term bs in Eq. (3) is neglected. The damping term bs is also neglected because its influence is small and difficult to determine (Xia et al. 2018). Then, the equation of motion of the bridge system is as follows (Eq. (4)): Equations (1) and (4) constitute the motion equation system of train and bridge systems. When the wheel-rail interaction model used is non-linear, as in this study, vb and bv become the non-linear functions of the motion states of the train-bridge coupling system (Olmos and Astiz 2018b). In this case, Eq. (4) cannot be solved by the decomposition method. This decomposition method separates bridge behavior into dynamic and pseudostatic responses, and it is the most used method for analyzing long structures subjected to non-uniform seismic ground motions.
In this study, the mentioned direct solution method (Xia et al. 2018) is used. This method is based on displacement inputs. That means that ground displacements at bridge support points are needed in order to obtain the seismic loads for exciting the train-bridge system. − bs ⋅ s (t) , in Eq. (4) are these seismic loads applied directly on the node close to each support's base in the form of concentrated forces and concentrated bending moments.
The bridge, train, and earthquake models that are involved in Eqs.

Bridge model
To study the dynamic response of this high-pier railway bridge, a finite element model has been built (Figs. 2, 3). The materials of the structural elements in this model have a linear elastic behavior.
Two types of finite elements have been considered in this model: • 1D finite element. This element is a beam finite element. Its cross section may vary along its axis. The shape functions of this element are cubic. Its material behavior is elastic. It has two nodes and six degrees of freedom (d.o.f.) per node. The piers and deck have been modelled with this type of element. • Punctual element with mass. It has six d.o.f.. Its mechanical properties are the mass and inertia of mass according to the element axes. This element has been used at the deck nodes to take into account the ballast mass, sleepers mass, rails mass and the rest of the non-structural mass over the deck. The damping matrix of the bridge dynamic system, b (Eq. 2), is obtained using the modal superposition technique considering the first 336 vibration modes with a damping ratio equal to 1% in all of them. With these vibration modes, the highest frequency considered is Fr 336 = 26.22 Hz, and the mobilized mass percentage in each direction (x, y, and z) is greater than 99%.
The track model used in this study is described by Olmos and Astiz (2018b). The track has been considered infinitely rigid, as in many other train-bridge interaction works available in the literature. The rail displacements, velocities and accelerations at the wheel-rail contact points are obtained through kinematic relationships of rigid section considering the movements of the corresponding deck sections.
In this model, the relative rotation around the vertical axis between the deck, on support sections, and the pier-top is allowed, except for piers no. 8, 9, 10, and 11 (see Fig. 3). Rotation around the vertical axis during a seismic event is also not allowed at deck sections on the bridge abutments. The described movement constraints reflect the real support conditions of the bridge deck.   Table 1 presents the deck mass and stiffness properties of this bridge that has a double composite action in the areas next to the deck supports on piers and abutments.
The modal analysis carried out on the model showed the following bridge frequencies: There were 10 frequencies between 0.34 and 1.270 Hz, which correspond to lateral bending modes. There were another nine frequencies between 1.363 and 1.958 Hz that correspond to deck vertical bending modes. The first and the second shape of the bridge vibration modes are shown in Figs. 4 and 5.

Track model outside the bridge
When the train is travelling outside the bridge, the track, which has been considered as infinitely rigid, moves following the ground seismic motion. Therefore, outside the bridge, the earthquake has been introduced in this study as imposed displacements on the rails (Eq. (5)): In Eq. (5), y L D , y R D are lateral displacements of the left and right rails respectively, and u gy,1 (t) is the ground lateral displacement caused by the earthquake in abutment 1.

Vehicle model
The train used for this study has been described by the authors in a previous work (Olmos and Astiz 2018b) This train is an articulated train type Thalys (or TGV-Alstom) of 10 vehicles. The convoy has two locomotives, two transition-passenger cars and six passenger cars. This train has a length of 193.15 m, from the first wheelset to the last one. The vehicle dimensions of this train are given in Fig. 6. In this study, a multibody model of this train has been built in order to obtain a dynamic train response. The model consists of three types of rigid bodies, (cars, bogies, and wheelsets), and two suspensions systems. The cars and the bogies in the model are joined by means of the primary suspension, whilst the bogies and wheelsets are joined by the secondary suspension. Both suspension systems are made up of linear springs and viscous damper units. The springs and dampers are arranged both horizontally and vertically. There is a third suspension on articulated trains such as this one. This third suspension has horizontal damper units that join car bodies in the train longitudinal direction. Figure 7 shows the multibody model of the train articulated part. All cars except the two locomotives are shown in this Fig. 7  The vertical movement, z W , and the rotation around track longitudinal axis, Rx W , are not independent d.o.f. in a train wheelset with the model used here. That is because these two movements of a wheelset depend on track movements, on the relative displacement between wheels and rails, and on track geometrical irregularities.
The geometry, stiffness, and damping parameters of the springs and dampers in the train suspensions, and the rigid bodies' mass of this train are shown in the scientific literature (Kwark et al. 2004;Song et al. 2003;Lee and Kim 2009).

Wheelset-track interaction model
The validation and the formulation of the wheelset-track interaction model used here can be found in three previous works of the authors (Olmos andAstiz 2013, 2018a, b).
The wheelset-track interaction model used in this study takes into account: the relative displacement between wheel and rail; the geometry of the UIC-60 rail and the S-1002 wheel including its flange; and the adhesion in the contact between the wheel and the rail. This interaction model is therefore a non-linear model and is based on Kalker's linear contact theory (Kalker 1979) and the modification of this theory that Shen et al. (1983) made in their non-linear contact model.
The properties of the track irregularities that have been considered in this study are shown in Table 2.
The quality of the track with the geometric irregularities described corresponds to a track class 5 F.R.A. (Frýba 1996). The standard EN-13848-5 (2009) considers this track class that belongs to a not well-maintained high-speed track. In this way, the study errs on the side of safety since the tracks of high-speed lines are well maintained.

Seismic ground motion at bridge supports
In the present work, a spectral-representation-based algorithm (Shinozuka and Deodatis 1988;Li and Kareem 1991) and an iterative scheme of Deodatis (1996) have been used to generate seismic ground motion (acceleration) time histories at each bridge support that are compatible with prescribed Eurocode 8 (EN-1998-1 2004) response spectra. The acceleration time histories at bridge supports are correlated according to a given coherence function; they include the wave propagation effect, and they have a specified duration of strong ground motion. Ground motion time histories are modelled as a uniformly modulated non-stationary stochastic vector process.
The spectral-representation-based algorithm and the iterative scheme of Deodatis (1996) are explained below.

A spectral-representation-based algorithm
u gy,j (t) and ü gz,j (t) are the ground y-direction horizontal acceleration time history and ground vertical acceleration time history respectively of the j bridge support during an earthquake.
The component of the two non-stationary stochastic vector processes can be expressed as (Eq. (6)): where n is the number of bridge supports, (no. of piers + 2), and a gy,j (t) and a gz,j (t) are the components of two stationary stochastic acceleration vector processes. A j (t) are the modulating functions. It is said that the stochastic vector processes are uniformly modulated because the modulating functions are independent of the angular frequency . The modulating functions adopted in this work are based on the Jennings et al. (1968) model and are detailed in Eq. (7).
In the last expression, 1j is the distance between bridge support no. 1 (first abutment) and bridge support no. j. c s is the velocity of the seismic wave propagation in the ground surface. c s is 2000 m/s in this study. With the formulation of Eq. (7), the wave propagation effect is properly considered.
The acceleration time histories, a gy,j (t) and a gz,j (t) in Eq. (6), have a mean value equal to zero and cross-spectral density matrices given by Eq. (8) and following: These two-vector processes, a gy,j (t) and a gz,j (t) , are stationary because their stochastic properties do not depend on any instant t. The spectral density function S k jm ( ) , in Eq. (7), only depends on the angular frequency . That is (Eq. (11)): In Eq. (9) of the cross-spectral density functions, the factor Γ ji ( ) is the complex coherence function of Eurocode 8, which is defined in Eq. (12): The coherence function penalizes the cross-spectral density corresponding to the points j and m when they are far away. That is to say that when two points are far away the influence of the ground motion generated at one point on the ground motion generated at the other one will be small. d jm is the distance between these points, and d jm,p is the distance between these points projected in the direction of the seismic wave propagation. In the case of this work, both distances, d jm = d jm,p , are the same because the seismic wave propagation direction is assumed to coincide with the bridge longitudinal direction. In Eq. (12), is a constant and i means √ −1 . With this coherency function the loss of the coherence effect is considered in the non-uniform seismic ground motions generated for this study.
For the simulation of each component of the ground motion, in y-horizontal direction and in z-direction, at each bridge support j, the following expressions can be used: In Eqs. (13), (14), and (15), N is the number of frequencies considered, a large enough positive integer number; up is the cut-off frequency or the highest frequency considered; Δ = ( up ∕N) is the frequency increase; Φ 1l , Φ 1l , … , Φ jl they are independent sequences of random phase angles generated from a uniform distribution in the interval [0, 2 ] ; and H k jm ( l ) is a typical element of the triangular matrix obtained from the Cholesky decomposition of the cross-spectral density functions matrix k ( ) (Eq. (8)). That is (Eq. (16)): As Ding et al. (2011) and other authors such as Cao et al. (2000) state, the cost of the digital generation of sample functions, ü gk,j (t) in this case, in a simulated stochastic processes, can be drastically reduced using the Fast Fourier Transform (FFT) technique. Equation (11) can be rewritten as follows (Ding et al. 2011): is a vector of M components for each jm subscript which can be computed as a discrete FFT of the vector B k jm (l ⋅ Δ ) (Eq. (18)): In this study, Eq. (17) and the FFT technique have been used. On the other hand, in order to generate time history sample functions compatible with a Eurocode prescribed response spectrum for each bridge support, the Deodatis (1996) scheme was implemented in the code developed for this work.

Iterative scheme of Deodatis
The target response spectra considered for generating the sample functions are described below. In these spectra, the two types of terrain found at the bridge location have been considered. For piers 8-11, the terrain consists of loose gravel soil in this study. In the rest of the bridge supports, including the abutments, the terrain considered is deep stiff clay soil.
The target elastic (horizontal or vertical) acceleration response spectrum of Eurocode 8 (EN-1998-1 2004) (considering 5% of viscous damping only as a reference), for bridge support j has the following structure: The parameters of Eq. (19) are given in Table 3 for the bridge supports location. The methodology used in this study (Table 4), is explained by Deodatis (1996). It starts by using the spectral-representation-based simulation algorithm described in Sect. 6.1 to generate stationary time histories ( a gy,j (t),a gz,j (t) ), that are compatible with a coherence function Γ jm ( ) and with a velocity of wave propagation but not with the predicted respond spectra (of Eurocode 8). After multiplying these stationary time histories by an appropriate modulating function A j (t) to introduce non-stationarity, the methodology of Deodatis upgrades the spectral density functions S k j ( ) of the vector process components as shown in Table 4, generates new stationary time histories according to the upgraded cross-spectral density matrix, and multiplies them again by the modulating function to introduce non-stationarity. This upgrading or iterative procedure is repeated several times, 10 or less than 10 times, to make the time histories response spectrum similar to that prescribed by Eurocode 8 (EN-1998(EN- -1 2004. Table 3 Target acceleration response spectra parameters for y-horiz. and vertical ground motion to generate a As y-horizontal and vertical seismic ground motions will be used at the same time to excite the bridge, RSA In order to calculate the acceleration response spectra corresponding to the ground acceleration time histories generated ü gk,j (t) at step 5 of Table 4 scheme the following Equations can be used. Table 4 Deodatis iterative scheme to simulate response spectrum compatible acceleration time histories at the n bridge support points on the ground surface

Number of iterations
End equal to 10?
Step 6. Upgrade spectral density functions as: Equation (20) governs the motion of a linear single degree of freedom system, with natural frequency l and damping ratio , subjected to ground acceleration ü gk,j (t) . Fixing the damping ratio in = 0.05 as a reference, Eq. (20) can be solved by the Newmark method for every support (j = 1,…, n) and for every frequency l considered in the study (l = 1,…, N). Equation (21) defines the acceleration response spectra corresponding to a ground acceleration time history. Zerva (2009) states that the convergence criteria for the Deodatis iterative scheme cannot be readily established, but, instead, that the iterations can be stopped when the shape of the target response spectrum is captured by the response spectrum of the simulations. Generally, only a few iterations, 10 or less than 10 iterations, are required in order to achieve this objective. If the mentioned objective is not achieved within the first few iterations, the simulated time history should be discarded and a new one generated (Deodatis 1996).
The iterative scheme of Deodatis (Table 4) to generate the time history ground motions at each bridge support ü gk,j (t) has been applied in the present study, performing 10 iterations. In each new iteration, the power spectral density functions S k j ( ) have been updated as indicated in Eq. (22): In Eq. (22), RSA k,j ( l ) is the target response spectrum of Eurocode 8 corresponding to the bridge support j, and RSÄu g (t) k,j ( l ) is the calculated response spectrum from the time history ü gk,j (t) generated in the previous iteration (i-1). From all 10 iterations carried out, the time history corresponding to the iteration in which the differences between the target and calculated response spectrum were less than 15% (Eq. (23)) have been selected as the most suitable for the present study.
The authors observed that the last iteration, number 10, was not always the one that satisfied the criterion expressed in Eq. (23). Sometimes, the time histories of iteration 8 satisfied the established criteria for obtaining a suitable time history response spectrum (Eq. (23)), and in the following iterations, (9 and 10), this criterion was not satisfied. Therefore, the authors concur with Zerva's claims.

Consistency treatment of earthquake record
As mentioned in Sect. 2, the method used in this study to obtain the train and bridge response under a non-uniform seismic ground motion is based on displacement input. So, displacement time histories of vertical and lateral ground motion at each bridge support are necessary to compute the seismic action. Xia et al. (2018) explained that the direct integration of the generated acceleration time histories may obtain velocity and displacement time histories in error trend lines. Those authors stated that the seismic acceleration time histories need to be modified. The aim of this modification, known as seismic record consistency treatment, is to make the seismic ground velocities and displacements at the initial and the final earthquake moment close to zero. The time domain high-pass filtering method has been used here as the mentioned record consistency treatment. This method consists in modifying the generated acceleration time histories by a high-pass filter with the form of a critical damping oscillator (Xia et al. 2018).
In Eq. (24), ü gk,j o (t) refers to ground acceleration time histories at bridge supports before the modification; whilst ü gk,j (t),̇u gk,j (t), u gk,j (t) are the acceleration, velocity, and displacement time histories at bridge supports after the treatment. c is the angular frequency of the filter.

Seismic ground motion simulations
In order to take into account the uncertainty, with the formulation and procedure described above, as well as the parameters shown in Tables 3 and 5, 20 sets of samples of ground displacement and acceleration time histories were generated at the bridge with the same earthquake intensity PGA y = 0.08 ⋅ g g. Thus, the ground motion of 20 different non-uniform seismic events with the same intensity were simulated. Each one of these events included concomitant movements in y and z directions at each bridge support.
Of all the 20 seismic events generated on the bridge, the one that produced the most adverse effects on traffic safety was selected, considering the derailment factor index and wheelset-track lateral force index that are defined below in this paper. Twenty train travel simulations, with the trains-bridge interaction model described here,were performed over the bridge at 240 km/h in order to determine which of them was the most damaging to traffic safety. Each one of these simulations was carried out with one of the 20 seismic events (24) u gk,j (t) + 2 ⋅ c ⋅̇u gk,j (t) + 2 c ⋅ u gk,j (t) =ü gk,j o (t), k = y, z; j = 1, 2, … , n  At bridge support no.1 (abut. 1) At bridge support no. 12 (pier 11) Fig. 10 a Generated ground vertical accelerations at bridge supports 1 and 12 (concomitant with ground lateral accelerations shown in Fig. 9) and b ground vertical displacements belonging to ground accelerations shown previously in this figure generated, and the aforementioned traffic safety indices were measured. Some results of this analysis carried out are shown in Appendix 1.
The selected seismic event from the 20 samples was the one used in this study because it was the most unfavorable for traffic safety. We shall call this selected seismic event the reference non-uniform seismic event. From this reference seismic event with a PGA y = 0.08 ⋅ g , other events have been obtained with different earthquake intensities ( PGA y = 0.10 ⋅ g, 0.06 ⋅ g, 0.04 ⋅ g, 0.03 ⋅ g, 0.02 ⋅ g, 0.01 ⋅ g ) by adequately scaling its acceleration time histories and subsequently applying the consistency treatment described in Sect. 6.1.3 to obtain the corresponding velocities and displacements time series. Figure 8 shows the agreement obtained between the target and generated horizontal acceleration response spectra for several bridge supports corresponding to the mentioned reference non-uniform seismic event. In this Fig. 8, the horizontal acceleration response spectra RSA y versus the period has been plotted.
In Figs. 9 and 10, ground acceleration and displacement time histories at several bridge supports corresponding to the sample seismic event with reference PGA y = 0.06·g are plotted. Figure 9 shows horizontal movements and Fig. 10 shows vertical movements.
In Figs. 9b and 10b, the amplitudes of the ground displacements at bridge support 12 (pier 11), were greater than those of support 1 (abutment 1). This was due to the terrain in pier 11 being less stiff than on the bridge abutment 1, and therefore the corresponding response spectra at these two locations were different. In these figures, the effect of variable site conditions of this spatially variable seismic ground motion can be noted. On another matter, in these two Figs. 9b and 10b, the wave passage effect can also be observed. The seismic ground motion began earlier at abutment 1 than at pier 11.

Type of earthquakes generated
The spatial variation of the seismic ground motion is characterized by taking several effects not considered in a uniform earthquake (same ground movement at all bridge supports) into account. These effects are the wave passage effect, the loss of coherence between the ground movements at two different points, and the site effect cause by different types of terrain at each bridge support.
In order to study the influence of these effects on the bridge and train responses separately, several types of earthquakes have been generated: type I, II, III, and IV earthquakes.

Type I earthquakes
Type I earthquakes correspond to the fully variable ground seismic movements along the bridge. This kind of seismic movements have been generated following the procedure described above in this Sect. 6 and summarized in Table 3. This type of seismic motion considers all the described effects of the spatial variation of seismic motion along the bridge.
In these earthquakes, the wave passage effect was considered by means of the modulating functions (Eq. (7)). In this Eq. (7) the speed of passage of the seismic wave on the ground surface was c s .
The loss of coherence between the ground movements at two different points was taken into account by means of Γ ji ( ) the complex coherence function of Eurocode 8, which is defined in Eq. (12), and it is a factor of the cross-spectral density functions (Eq. (10).
The site effects considered in this kind of seismic movements are those indicated in Table 2. The terrain at supports j = 9 to j = 12 (piers 8-11) is loose gravel (type D ground according to , and the terrain at the rest of the supports consists of stiff clays (type C ground according to Ec-8).

Type II earthquakes
Type II earthquakes are fully uniform seismic ground movements, i.e., the ground accelerations and displacements time histories are the same at each bridge support in an earthquake of this type (Eqs. (25) and (26)).
In this kind of earthquake, no effects of the spatial variation of seismic motion along the bridge are considered.
These seismic movements have been obtained by copying the generated ground movements at the support j = 1 of the type I earthquakes.

Type III and IV earthquakes
Type III earthquakes are uniform seismic ground movements but considering the wave passage effect. Type IV earthquakes are variable ground seismic movements, although the terrain considered at each bridge support is the same: stiff clays.

Earthquake outside the bridge
In all the earthquake types described above, the ground seismic movement occurs at each bridge support (piers and abutments), and ground seismic movements also occur under the track outside the bridge, across which the train travels before entering the bridge and after leaving the bridge. The seismic movements of the track, located before (25) (ü gy,j ;ü gz,j ) = (ü gy,1 ;ü gz,1 ) j = 2, … , n (26) (u gy,j ;u gz,j ) = (u gy,1 ;u gz,1 ) j = 2, … , n Fig. 11 Forces on bridge deck caused by ground seismic movements ( u gy,1 ;u gz,1 ) of abutment no. 1 and after the bridge, have been modeled as a uniform earthquake equal to the seismic movements of the terrain in abutment 1 and abutment 2, respectively. The mentioned type of earthquake has been considered outside the bridge since it is a reference calculation that is used to evaluate the influence of the bridge on the response of the train to the earthquake.

Seismic forces on the bridge
As mentioned in Sect. 2, the seismic forces acting on the bridge are functions of the ground vertical and horizontal displacements at each support (piers and abutments): The term − bs ⋅ s (t) in Eq. (27) refers to these seismic loads applied directly on the nodes close to the supports' base in forms of concentrated forces and concentrated bending moments. In Figs. 11 and 12, these seismic forces are represented at the mentioned nodes of the finite element model. These seismic forces acting on the nodes of the bridge model close to the bases of the piers and abutments can be expressed as follows (Eqs. (28) to (30)).
(27) s (t) = (… , u gy,j (t), u gz,j (t), …) j = 1, 2, … , n Fig. 12 Forces on a pier caused by ground seismic movement ( u gy,j ;u gz,j ) at this location In the Equations above, l d is the length of the first and last deck finite element (Fig. 11). I zd , I yd are section inertia moments of the first and the last deck finite element. h e,j is the length of the first finite element in the pier j − 1 (support j): the one finite element that is in contact with the support foundation (Fig. 12).A p,j, , I xp,j are the section area and section x-inertia moment of the first finite element in the pier j − 1 (support j). E d and E cp are the deck and piers longitudinal elasticity modulus considered in the model, respectively.

implementation and numerical solution
The implementation and numerical solution of the train-bridge dynamic interaction model used in this study has been described in detail in a previous work by the authors (Olmos 2016).
Several scripts in Matlab language have been used. These scripts build the bridge equation system Eq. (4), using the matrices b , b , b previously obtained with ANSYS software. In addition, the code in Matlab completely builds the train equation system Eq. (1) and obtains the numerical solution of the problem.
The mathematical problem to be solved consists of two systems of non-linear and coupled second-degree differential equations. In this problem the non-linearities are associated with the wheelset-track (or wheel-rail) interaction model. For the systems of equations of the train (Eq. (1)) and of the bridge (Eq. (4)), in the present study, the non-linear functions are located in the interaction force vectors acting on the train and on the bridge ( vb and bv ).
To obtain the numerical solution (dynamic response in the time domain) of a dynamical system from its equations of motion requires the use of a finite differences method to integrate the equations mentioned along time. In this study, the Wilson (2002) time step θ-method was used for that purpose. Since they are non-linear equations, Wilson's method was combined with the Newton-Raphson method.
However, the problem here concerns two coupled dynamical systems, not one. In order to solve the coupling between the two systems (train and bridge), an iterative procedure has been used by which, at every time step of the analysis, and in every iteration of the Newton-Raphson procedure to solve the bridge system of equations, another Newton-Raphson procedure was used to solve the train system of equations.

Running safety analysis performed
In the first part of this work, we have studied the bridge response (accelerations and lateral displacements) and train response (accelerations and traffic safety indices) when the train traveled on the bridge and a moderate spatially variable earthquake occurred. Several simulations or train travels over the bridge with a seismic event were performed. In the simulations the train travelled at the following speeds 180, 200, 220, 240, 260, 280, 300, 320, 340, and 360 km/h. For each train travel speed, earthquakes of peak ground acceleration (PGA) 0.01 g, 0.02 g, 0.03 g and 0.04 g were tested.
In the second part of this work, we studied the influence of several factors that characterized the spatial variation of seismic movement on the train and bridge responses. The purposed of this was to qualitatively determine the errors that can arise from not considering these factors in the characteristics of the earthquake.
In all these simulations, the earthquake began 2.5 s before the train arrived at the bridge and the simulation ended when the last car of the train left the bridge. In this way, the bridge was excited by the earthquake when the train entered it. In addition, for all train speeds tested, the train had wheelsets within the bridge at least between second 2.5 and

Influence of earthquake intensity and train speed on the bridge and train responses
This section shows the outcomes of the simulations carried out with fully variable earthquakes (type I). Ground motions of different intensities and train traveling at several speeds were considered. Figure 13 shows maximum deck lateral accelerations in its sections on pier 9 (a) and on pier 15 (b), respectively. It was observed that the peak ground acceleration of the earthquake, i.e., earthquake intensity, was a factor that decisively influenced the deck lateral accelerations. An earthquake of PGA equal to 0.04 g caused deck accelerations that were four times greater than those caused by an earthquake of 0.01 g PGA (Fig. 13b). However, train speed did not seem to influence the bridge deck accelerations when the bridge was subjected to an earthquake.
Locomotive and passenger car maximum lateral accelerations versus train speed are shown in Fig. 14. Earthquake intensity also influenced these accelerations. This time, the acceleration curves presented some peaks. For example, in Fig. 14b, maximum lateral acceleration curves for locomotive body presented peaks at 280 km/h for all earthquake intensities tested in the simulations. At that speed a resonance phenomenon occurred due to the interaction between the bridge, the vehicle, and the earthquake. The resonance was more evident with a higher earthquake intensity. When the train speed was v = 280 km/h (77.77 m/s), the lateral acceleration of the locomotive car was higher. An explanation of this phenomenon can be that at 77.77 m/s the locomotive lateral movement, with a frequency of Fr = 0.643 Hz, had a wavelength (L = v/Fr) of L = 121 m which approximately matched the wavelength of the bridge vibration mode 43 (Fig. 14c). Figure 15 shows the higher values of two traffic safety indices, i.e., their worst values. The traffic safety indices considered are the derailment factor of a train wheelset and the maximum wheelset-track lateral force, which is known as the Prud'Homme criterion. As described in the standard EN-14363 (2005), the derailment factor DF of a train axis is the ratio between the lateral Y and vertical Q forces that one of its wheels exerted on the rail. For safety reasons, this factor should not exceed 0.8.
In the Prud'Homme index (Eq. (31)), P 0 is the static vertical force of a train wheelset in kN, and ΣY is the lateral dynamic force of a train wheelset on the track, i.e., the sum of lateral forces of the two wheels on the rails. This criterion limits the risk of the track sliding over the ballast. The Prud'Homme index should be less than one to ensure there is no risk for traffic safety.
The European standard EN-14363 (2005) explains how to calculate these two indices in terms of frequencies to be considered. This standard has been followed to calculate these two indices in the present study.
In Fig. 15, it can be observed that the most unfavorable traffic safety index is the derailment factor. That is because with an earthquake of only 0.03 g PGA there are train speeds that fail to meet the safety criterion (DF < 0.8) (Fig. 15b). Most train speeds are unsafe with an earthquake of 0.04 g PGA. We could ask ourselves what the probability of a 0.03 g PGA earthquake occurring on this bridge when a high-speed train is traveling over it at 300 km/h for example is. This refers to the occurrence probability of a situation that is incompatible with traffic safety. This bridge is located in Álora municipality, Malaga province (Spain). At that location, the value of the reference peak ground acceleration on type A ground a gR is 0.153 g, according to the Spanish Annex of the Eurocode 8 (En-1998-1 2004). The peak ground acceleration corresponding to the design seismic action is: The importance factor I of the bridge and the coefficient S of the terrain type C (firm clays) have been considered in Eq. (32). This design earthquake has a return period of T NC = 1073 years, and an occurrence probability during the bridge service life (100 years) of Pr 0 = 8.90%.
According to this, the occurrence probability of an earthquake of 0.03 g PGA in a year at the site of this bridge is Pr 1 = 24.38%, (during the bridge life is Pr' 1 = 100%), with a return period of T R1 = 3.6 years. On the other hand, if the train passage frequency (at 300 km/h) over the bridge was one train every hour (with service from 07.00 to 23.00 h), it could be estimated that the probability of one of these trains being on the bridge at a given moment would be Pr 2 = 0.37%. We can consider that the seismic event of 0.03 g PGA and the passage of the train over the bridge are independent events. So, the probability of this earthquake occurring during the bridge life while a train is passing over the bridge is Pr 3 = Pr' 1 · Pr 2 = 0.37%.
But what would be the outcome if this bridge was in a different location? For instance, in the higher seismic risk province of Granada, a gR = 0.22 g and a PGA DSA = 0.33 g.
It has been proven that this bridge is able to withstand this new design earthquake with PGA DSA = 0.33 g using the right reinforcement in the pier cross sections. In this case, the occurrence probability in a year of a frequent earthquake with PGA equal to 0.03 g would be Pr 4 = 54.43%. The occurrence probability during the whole bridge life (100 years) of a frequent earthquake with PGA equal to 0.03 g would be Pr′ 4 = 100%. If we increase the train passage frequency over the bridge (at 300 km/h) to one every 10 min, with the same service interval (from 07.00 to 23.00 h), the probability of the train being on the bridge at a given moment would now be Pr 5 = 2.21%. So the probability of the frequent earthquake of 0.03 g PGA occurring during the bridge life while a train is on the bridge located in Granada is also Pr 6 = 2.21%.
The formulation used in the calculations of the mentioned probabilities is shown in Appendix 2 (Table 6).
These numbers raise a question: what risk level must be taken in the design of a high-speed viaduct to consider the bridge safe for high-speed train traffic during frequent earthquakes?
We can calculate the probability (Pr 9 ) of a train traveling over track (outside the bridge) at 300 km/h that is at risk of derailment due to an earthquake as a reference. In the present study, numerical simulations have been carried out to obtain the intensity of the earthquake capable of causing derailment factors exceeding 0.8 in a Thalys-type train traveling at 300 km/h over track (outside the bridge). For these simulations, the track model described in Sect. 3.2 has been used. An earthquake that endangers the safety of the Thalys train traveling at 300 km/h over high-speed track has a PGA of 0.35 g. The probability of a train of this type traveling at the aforementioned speed that is at risk of derailment because an earthquake occurs is     Table 7 of Appendix 2. For the calculation of Pr 9 the following issues have been taken into account: the earthquake of PGA 0.35 g or greater occurs with an area of action of 15 km along the track in Granada province; the train frequency is one every 10 min in a service interval of 07.00 to 23.00 h.
If we consider that in the area of action (15 km long) of both earthquakes (the one of 0.03 PGA and the other of 0.35 PGA) there are two bridges like the one studied, the probability of the train traveling over a bridge at 300 km/h that is at risk of derailment because a frequent earthquake occurs is now Pr′ 6 = 4.42%. Therefore, we can state that the existence of this type of bridges on a high-speed line with the seismicity of Granada province can multiply by 4 (Pr′ 6 = 4.42% > Pr 9 = 1.11%) the risk of derailment due to seismic movements.
The authors of this paper propose the following rule for the design of a high-speed railway viaduct in seismic areas: The risk of derailment of a high-speed train over a viaduct when an     (frequent) earthquake occurs must be smaller than the risk of derailment of a train travelling over track (outside bridges) caused by an uncommon earthquake.
We wish to emphasize that the train passage frequency and the duration of the train travel over the bridge are necessary parameters when calculating the occurrence probability of this new seismic serviceability limit state proposed.

Influence of wave passage effect, loss of coherence, site effect and bridge lateral stiffness
In order to study the influence of the different effects of the non-uniform seismic action on the bridge and train responses separately, the different types of earthquakes generated (Sect. 6.2) have been tested in several dynamic simulations.

Bridge response
The values of several bridge and train response parameters (accelerations and displacements) obtained from the dynamic simulations carried out are shown in Figs. 16, 17, 18 and 19. Each simulation corresponds to a train traveling over the bridge at a certain speed when a type of seismic ground motion with PGA equal to 0.03 g occurs. The seismic action types considered are I, II, III, and IV. These earthquakes have already been described in Sect. 6.2. We can name these earthquakes as follows: type I earthquake is fully variable; type II earthquake is fully uniform; type III earthquake is uniform but considering the wave passage effect; and type IV earthquake is variable but with the same site conditions. Additionally, the type I earthquake has also been tested using a bridge model with piers half the height of those of the Arroyo de las Piedras viaduct piers. Figure 16 shows the maximum deck lateral displacements in support sections over viaduct piers numbers 9 and 15, when the train runs at several speeds and when several types of earthquakes occur.
In Fig. 16a, we observe that the maximum deck lateral displacement over pier 9 caused by the type I earthquake (line with circles) is 28 mm. This type I earthquake takes into account all the effects of the spatial variation of seismic action. The type IV earthquake differs from the type I earthquake in that the latter considers the ground under pile 9 to be loose gravel rather than stiff clays. For the type IV earthquake, Fig. 16a shows lower deck lateral displacement values than for the type I earthquake. In other words, the loose gravel considered in the type I earthquake amplified the displacements of pier 9. The type II earthquake (line with stars), which is spatially uniform, (i.e., it is the same in all the supports of the bridge), produced deck displacements values over pier 9 that were lower than those of the type I earthquake. In this Fig. 16a it is also observed that the type III earthquake (line with rhombuses), (which is like the type II earthquake but considering the seismic wave passage effect), produced greater displacements than the type II earthquake.
The simulations that produced the smallest deck lateral displacements on pier 9 were those in which the fully variable earthquake (type I) acted but on the bridge model with the shortest piers.
In Fig. 16b, we observe that the differences in deck displacements on pier 15 produced by each type of earthquake tested were small. Although in this case, the deck displacements on pier 15 that occurred due to the uniform earthquakes were slightly greater. Nevertheless, there was a great difference between the displacements obtained in the simulations carried out with the Arroyo de las Piedras viaduct model and those produced in the simulations with the shorter pier viaduct model.
From the results observed in Fig. 16b, it can be stated that the viaduct height had a great influence on the lateral displacements that occur when a frequent earthquake acted. Furthermore, earthquakes generated that considered the spatial variation of seismic action (type I) produced greater displacements in high piers than those that did not consider this variation (type II) (Fig. 16a).
It is also observed that train speed had negligible influence on the deck lateral displacements. Figure 17 shows the maximum deck lateral accelerations in support sections over several piers, when the train traveled at several speeds and when several types of earthquakes occurred. Figure 17a shows that the maximum deck lateral accelerations over pier 9 that occurred with the type I and type IV earthquakes. There were no differences in deck lateral accelerations on pier 9 produced by the type I and type IV earthquakes. The lowest accelerations in Fig. 17a are shown in the curve corresponding to the simulations with the short-pier viaduct model (line with x). Therefore, the lateral flexibility of the viaduct is a factor that has an important influence on the response in deck accelerations under frequent earthquakes. Figure 17b also shows deck lateral accelerations, but this time of its section on pier 14, which was shorter than pier 9. In this case, although the deck lateral accelerations were greater on pier 14 than on pier 9, reaching 1.78 m/s 2 , all types of earthquakes tested produced similar lateral accelerations at that point of the bridge deck.
In Fig. 17c, the bridge deck response corresponds to the section above pier 17. In this case, it was observed that the types of earthquakes that produced greater accelerations were the uniform ones. One explanation for this observed outcome is that the short pier 17 was located next to abutment 2. Vibration modes that excited this pier keep the bridge deck fixed on abutment 2. When the variable earthquake acted, the ground at abutment 2 moved differently from the ground under pier 17, and this seems to have reduced the vibration amplitude of modes that excited pier 17 to some extent.
From the results of the simulations shown in Figs. 16 and 17, we can state the following: • About bridge lateral stiffness. Of all the effects studied in this section (wave passage effect, loss of coherence, site effect, bridge lateral stiffness, and train speed), bridge lateral stiffness was the parameter that had the greatest influence on deck lateral accelerations. The laterally stiffer the bridge was, the lower the deck lateral acceleration was at all speeds when a frequent earthquake occurred. In the graphs in Fig. 17a-c, the curve corresponding to the bridge with half the height of pier and the type I earthquake, (line with x), was the one that showed the lowest values of deck lateral accelerations. • About site effect. In this work this effect was studied by means of the type I and type IV earthquakes that considered the loss of coherence in the seismic movement at the different supports. In the type I earthquake, the real terrain under the bridge supports was considered, i.e., the terrain at piers 8 to 11 was loose gravel and the terrain at the rest of the supports was firm clay. In the type IV earthquake, firm clay at all supports was assumed as the terrain. In Fig. 16a, we can observe that the loose gravel considered in the type I earthquake amplified displacements at pier 9, as mentioned above. However, this amplification of the deck response, due to a less firm terrain at pier 9, was not observed in deck lateral accelerations (Fig. 17a). • About wave passage effect. This effect can be analyzed by comparing deck lateral movement outcomes obtained with the type III earthquake (line with rhombuses) and with the type II earthquake (line with asterisks). These earthquakes are uniform, i.e., without loss of coherence. However, the type III considered the passage of the seismic wave, whilst the type II did not. In the graphs in Figs. 16a and 17a, the values of the type III earthquake (line with rhombuses) were greater than those corresponding to the type II earthquake (line with asterisks). This means that the wave passage effect produced a higher deck response (lateral displacement and accelerations) on high piers. On shorter piers, (Figs. 16b,17b,c), the wave passage effect was not appreciable in the displacements nor in the accelerations. • About loss of coherence. We can analyze this last effect by comparing deck lateral movement outcomes obtained with the type IV earthquake (line with squares) and with the type III earthquake (line of rhombuses). These two types of earthquakes differ in that type IV considered the loss of coherence in ground movements at each bridge support and type III did not. We can observe that this loss of coherence increased the bridge response at pier 9 in term of accelerations (Fig. 17a), where deck accelerations due to the type IV earthquake (line with squares) were greater than those due to the type III earthquakes (line of rhombuses). However, this did not happen with the deck movements on lower piers located near abutment 2 (Fig. 17c). Figures 18 and 19 show graphs of bridge deck lateral accelerations on pier 9 and on pier 17 in the domain of time (a) and frequency (b) when the train traveled over the bridge at 260 km/h and the uniform type II earthquake or the variable type I earthquake occurred.
In Fig. 18b, in general the variable earthquake excited more the vibration frequencies of the deck movement on a high pier such as pier 9. On the other hand, Fig. 19b shows that some frequencies of the deck movement on pier 17 were more excited by the uniform earthquake. The boundary conditions of the structure may explain this phenomenon. Additionally, pier 17 is in this area of the bridge, near abutment 2, where there were greater deck lateral accelerations with uniform earthquakes (Fig. 17c).    Figure 20 shows the maximum lateral accelerations of passenger vehicles (a) and of locomotives (b) when the train traveled at several speeds over the bridge and when seismic movements of several types occurred. All these seismic movements had an intensity of 0.03 g PGA.

Vehicle response
In the case of passenger cars (Fig. 20a), it was observed that the spatial variation of the seismic action did not influence the vehicle response. This means that the type I earthquake (fully variable) produced the same lateral accelerations in the passenger cars for all train speeds as the type II earthquake. On the other hand, the speed of the train had a great influence on the maximum accelerations of passenger vehicles. To explain this, it must be stated that in this Thalys-type train, passenger vehicles have little mass, and therefore, little inertia. Their greatest lateral accelerations occurred when these cars entered the bridge. At that moment they went from an excitation due to ground movement outside the bridge to being excited by the seismic bridge response that was greater than the mentioned ground movement. This sharp increase in the response of passenger cars did not occur in the locomotive or in transition cars, because they had greater inertia and a different suspension system.
In the case of the locomotive, uniform earthquakes produced greater lateral accelerations in the locomotive body for most train speeds (Fig. 20b). This last fact can be explained as follows: locomotives are independent vehicles in their movements from the rest of the convoy, i.e., they do not have common wheelsets with the rest of the train or longitudinal suspensions between their bodies and that of the next vehicle in the convoy. The convoy's head locomotive obtained its greatest excitement when it traveled over the end area of the bridge deck, where it had greater lateral accelerations (Fig. 17c). That area of the bridge was the one that was most excited by the type II and III (uniform) earthquakes because of edge conditions. It is also observed, in Fig. 20b, that a resonance phenomenon occurred in the behavior of the locomotive when it traveled at 280 km/h. This resonance has been explained above in the comments regarding Fig. 14 Figure 21 shows the maximum values of the Prud'Homme index for the train wheelsets (a), and the maximum derailment factor values for the train wheelsets (b), when the train traveled at several speeds over the bridge and when seismic movements of several types occurred, all with an intensity of 0.03 g PGA. It should be recalled that the most realistic earthquake type was type I, fully variable.

Traffic safety indices
We have also included simulation results obtained with the ICE3-type train when it traveled over Arroyo de las Piedras viaduct under the type I earthquake (fully variable) at several speeds, (curve with triangles) in this Fig. 21.
This train was a conventional-type train of eight vehicles. Its vehicles are characterized by being passenger cars all of which had traction. For this reason, ICE3 train cars have greater inertia than Thalys-type train passenger cars. The mechanical properties of the ICE3-type train can be found in the documents of Antolín (2013) and Olmos (2016).
It has been verified that the highest values of the traffic safety indices studied were reached when passenger vehicles entered the bridge. Thus, the moment with the greatest risk to traffic safety corresponded to the entry of passenger vehicles onto the bridge.
From the graphs in Fig. 21 we can see that the maximum values of the traffic safety indices (Prud'Homme index and derailment factor) were very similar for any type of earthquake tested.
On the other hand, the lateral bridge flexibility does not seem to influence the maximum values of the derailment factor and the lateral track-wheelset force index (Prud'Homme index). The curve with x (in Fig. 21), which corresponds to the Thalystype train traveling over a viaduct with piers that are half the height of those of the Arroyo de las Piedras viaduct, had a similar layout to the other curves.
It can also be stated that the train mechanical characterization is one of the factors that most affects traffic safety when an earthquake occurs. In both Fig. 21a, b the curve with triangles is below the rest of the curves. This curve corresponds to the simulations carried out with the ICE3-type train whose vehicles have more mass, and therefore greater inertia, than the passenger cars of the Thalys-type train. The ICE3-type train had lower derailment factor values than the Thalys-type train, and therefore, better behavior under earthquakes.
From Fig. 21b we can see that Thalys train speeds exceeding 280 km/h are not safe when an earthquake with intensity PGA = 0.03 g occurs. At those speeds the derailment factors of the passenger car wheelsets are greater than 0.8. However, the ICE3-type train can travel safely over the viaduct at any of the speeds tested.

Conclusions
A spectral-representation-based algorithm and an iterative scheme of Deodatis have been used to generate seismic ground motion (acceleration) time histories at each bridge support which are compatible with prescribed Eurocode 8 response spectra. The acceleration time histories at bridge supports are correlated according to the Eurocode 8 coherence function, and include the wave propagation effect and site effects.
A train-bridge-earthquake dynamic interaction model has been built. This model considers the non-linearities of the wheel-rail contact. Since this is a non-linear dynamic interaction model, the formulation used to obtain the train and bridge dynamic behavior corresponds to the direct method of the accurate solution methods for obtaining the response of a structure subjected to non-uniform seismic ground motions in a framework of FEM. In this formulation, the seismic action is introduced into the differential equation of bridge motion by means of ground displacement functions at each bridge support.
From the parametric study carried out, the following conclusions can be drawn: • The traffic safety of a Thalys-type articulated train when traveling along the Arroyo de las Piedras Viaduct has been studied under the action of moderate earthquakes. The numerical simulation outcomes show that for earthquakes of intensity 0.03 g PGA, and for train speeds greater than 280 km/h the train travels over the bridge are not safe. The occurrence probability of this earthquake at this bridge location when a train passes is 2.21%, considering a passage train frequency of one every 10 min in an service interval from 06.00 to 23.00 h. If we consider that there are two bridges like the Arroyo de las Piedras Viaduct in the area of action of the earthquake, the occurrence probability of the event described will be 4.42%. The risk of derailment of this train traveling over the bridge when a frequent earthquake occurs may be four times greater than the risk of derailment of the train traveling on track due to a strong earthquake. The authors have proposed a new serviceability limit state (SLS) of traffic safety for the design of railways in this paper. This SLS can be formulated as follows: the risk of derailment of a high-speed train over a viaduct when a frequent earthquake occurs must be smaller than the risk of derailment of a train travelling over track (outside bridges) caused by a strong earthquake. • The lateral flexibility of the viaduct has been shown to have a great influence on the bridge response and on the vehicle response when frequent earthquakes occur. This flexibility increases the maximum lateral accelerations of the train locomotive. A bridge with piers half as high as the viaduct studied, reduces the deck lateral acceleration when an earthquake of 0.03 g PGA occurs by 35%. • The spatial variation of the seismic action, which can be understood as a sum of several effects (wave passage effect, loss of coherence, and site effect), produces an increase in deck and pier lateral displacements and accelerations in a large part of the bridge length. This increase in bridge response is more significant in the bridge area of higher piers. Loss of coherence and site effects are the effects that influence this increase most. All this leads to an increment in the internal forces at the piers due to this spatially variable seismic action. • The factors that most affect traffic safety during a frequent earthquake are the earthquake intensity, the train speed, and the vehicle type. The vehicles that presented a worse behavior under an earthquake when traveling or entering a bridge were passenger cars, because they have less mass. The most critical time for traffic safety during an earthquake is when the passenger cars of a train enter the bridge. However, factors such as the spatial variation of the seismic action or the bridge lateral flexibility do not seem to be decisive in the maximum values of the derailment factor.

Appendix 1
This Appendix 1 shows some results of the analysis carried out to determine the most dangerous seismic event, in terms of traffic safety, of the 20 events of PGA 0.8 g generated. These 20 seismic events were generated following the Deodatis procedure described in Sect. 6. In these Figs. 22 and 23, we can observe that seismic event number 4 caused a maximum lateral wheelset-track force index or Prud'Homme index (Eq. (30)) on the Thalys-type train travelling at 240 km/h over the bridge of IPH = 1.336, and a maximum derailment factor of DF = 1.007. As a reminder, the DF derailment factor of a set of train wheels is the relationship between the lateral Y and vertical forces Q of one of its wheels exerted on the rail. These traffic safety values, obtained with the seismic event number 4, correspond to the highest value of all the 20 events in the case of the IPH,