Numerical investigation of seismic behaviour of railway embankments in cold regions

To investigate more fully seismic behaviour of the Qinghai-Tibet railway embankment, a comprehensive discussion and analysis is conducted in this paper by applying a numerical technique. Specifically, the one dimensional equivalent linear ground response analysis was conducted in permafrost regions. On this basis, the seismic response of a typical railway embankment was further studied by applying the nonlinear dynamic finite element analysis method. As a result, nonlinear behaviour of permafrost sites was determined, and the dynamic acceleration, velocity and displacement of the embankment was discussed and the quantitative assessment was approximately estimated. The results indicate that the dynamic response of the embankment has distinct nonlinear characteristics. The peak ground acceleration coefficient at the embankment shoulder is larger than the natural ground surface, marking a 73% increase compared to the coefficient on the natural ground surface. When the seismic intensity reaches a certain value, a plastic zone gradually appears in the embankment, and a continuous extension of the plastic zone can be noted with an increase in peak acceleration of the input seismic wave. The findings of this research may provide an additional insight and have significant implications for further research of cold regions.


Introduction
Permafrost, or perennially frozen ground, refers to the ground that has had a temperature lower than 0°C continuously for at least two consecutive years. The Qinghai-Tibetan Plateau (QTP) is the highest and one of the most extensive plateaus in the world with the largest expanse of elevational permafrost on the earth. The area of permafrost on the QTP is estimated to be 1.3×10 6 km 2 [1,2]. Moreover, the QTP lies in one of the world's most active seismic zones: a total of 33 Ms 6.0-6.9 earthquakes and 3 Ms 7.0-8.5 earthquakes have occurred in this region since 1980 [3,4]. The Ms 8.1 Central Kunlun earthquake of 14 November 2001 produced a nearly 400 km-long surface rupture zone, with as much as 16.3 m of left lateral strike-slip along the active Kunlun fault in northern Tibet [5]. It has been observed that the forms of devastation in areas underlain by permafrost mainly include earthquake ruptures, fractures, liquefaction, seismic subsidence, and collapses. Many researchers have studied dynamic properties of frozen soil using both laboratory and field procedures. Permafrost is a kind of special soil that is highly sensitive to temperature changes [6]. The properties of permafrost are largely dependant on temperature. Even a relatively small change in temperature may cause a drastic change in dynamic response [7]. For instance, the Young's modulus of frozen soil is expressed in magnitudes that are tens to hundreds times higher compared to those of unfrozen soil, and the frozen soil is characterized by higher shear wave velocities in comparison with unfrozen soils [8][9][10][11]. The soil temperature exerts a large influence on the shear modulus and damping ratio. Such change in dynamic soil properties may significantly alter the seismic site response, as well as the seismic response of the infrastructure. However, only limited research has been conducted so far to investigate the effects of permafrost on the seismic site response. Vinson [12] studied the seismic site response of thick permafrost and concluded that the seismic response of a thick permafrost layer would be similar to that of the rock. Sritharan et al. [13] investigated the influence of seasonally frozen ground on the seismic performance of bridges. Yang et al. [14] studied ground motion characteristics in permafrost regions of Alaska by using one-dimensional equivalent linear analysis. Wang et al. [15] studied characteristics of ground motion at a permafrost site under various temperatures and summarized the influences of earth temperature on ground motion parameters. Chen et al. [16] discussed the impact of permafrost change on the seismic site response and concluded that ground motion parameters were significantly impacted by permafrost and active layer thickness. In previous related researches in permafrost regions, most seismic responses of permafrost sites are studied based on the one-dimensional model only, and are rarely connected to the response of the infrastructure in cold regions. Therefore, It is imperative to understand how permafrost affects ground motion characteristics at various ground motion intensity levels. Hence, the quantitative analysis of the site and the infrastructures should be further studied. This paper presents a quantitative assessment of the seismic response of sites considering the active layer thermal state. An embankment section situated in Beiluhe segment of Qinghai-Tibet Railway (QTR) was selected to represent a typical soil profile. One-dimensional equivalent linear analysis was adopted to analyse the seismic site response for various seismic hazard levels, and site responses in various seasons were determined. Then the seismic response of this embankment was determined by applying the nonlinear two-dimensional finite element method, and an approximate quantitative assessment was made. The results of this research will be of great theoretical and practical significance for resisting earthquake and preventing disasters in cold regions of China.

Site description
Based on in-situ investigations and experimental analyses, combined with previous engineering experience in permafrost regions, the mean annual ground temperature (MAGT) is wildly used in permafrost regions to evaluate thermal stability of permafrost. The MAGT is the ground temperature at the depth of a profile where the annual temperature range approximates zero. Moreover, the volume ice content (iv) serves as another key factor influencing mechanical properties of frozen soil [17]. Presently, the MAGT of −1°C is defined as the temperature boundary between warm and cold permafrost regions [18]. The permafrost with more than 20 % of volumetric ice content is defined as ice-rich permafrost [19,20]. In the light of this division, the QTR crosses about 550km of permafrost regions, including 275km of warm permafrost regions, 221km of ice-rich permafrost regions, and 124km of warm and ice-rich permafrost regions [21]. The Beiluhe segment of the QTR (34°51.26 'N, 92°56.35 'E) is located in the arid climate region of QTP, the freezing period of which spreads from September to April. The maximum, minimum and mean annual temperatures in this region are 23.2 °C, -37.7 °C and -3.8 °C, respectively [22]. The MAGT of the Beiluhe segment varies from -1.41 to -1.68 °C [23]. In environments containing permafrost, the active layer is the top layer of soil that thaws during summer and freezes again during autumn. Beneath the active layer lies permafrost, and the natural permafrost table of the segment is between -2 m and -3 m. The lithology in this segment is characterized by gravelly sand, silty clay, and mudstone. Thick ice underground layer near the permafrost table is abundant along the entire test section, and the temporal and spatial variations of thermal regimes of permafrost exist in the site. From the point of view of engineering geology, the Beiluhe segment can be classified as an unfavourable and poor engineering geological section [24].  The ground temperature and shear wave velocity (VS) curve of a typical borehole in Beiluhe section K1137+700 of the QTR, measured in July, 2012, is shown in Figure 2. It can be observed that the shear wave velocity of soil increases significantly near the upper limit of the permafrost. Moreover, the ground temperature at a depth of 12m is -1.5°C and the VS of the frozen soil can reach up to 620 m/s.

Seismic response model
Proper understanding of seismic site effects is crucial for improving assessment of seismic hazards. The one-dimensional equivalent linear method was adopted in this study [25,26]. It has been proven that calculation results are technically acceptable when the shear strain amplitude of soil in site is less than about 1-2 %, and the peak ground accelerations are less than about 0.3-0.4g [27]. This method is implemented on the basis of the elastic wave propagation theory, assuming vertical propagation of horizontal shear waves in soil deposits. Then the seismic wave propagation can be described mathematically by solving the wave equation. The soil profile is discretized into a multidegree-of-freedom shear beam system, as shown in Figure 3. Each layer m is assumed to be homogenous and isotropic and is characterized by the shear modulus G m , soil density ρ m , damping ratio λ m and thickness h m , and the normalized shear modulus G m / G max (G max is the maximum shear modulus) is used in presenting the modulus reduction curve. Then the iterative calculation is carried out when the definite calculating parameters are used, and thus the nonlinear problem is changed into a linear problem.

Numerical computation scheme
A numerical scheme is presented on the basis of the geology of Beiluhe section in order to compute seismic response in permafrost regions ( Figure 4). As shown in Figure 4, this one-dimensional site response analysis model contains three layers: the upper active layer, the ice-rich permafrost layer, and the ice-poor permafrost layer. Measured temperature curves of the Beiluhe section at different time in 2012 are shown in Figure 5. The active layer is frozen in cold season, and then it thaws in warm season, while permafrost is always frozen and the temperature is nearly -1.5 °C. Therefore, the thermal state of the active layer is taken into account in the following calculation. The active layer temperature in cold seasons is -1 and -5 °C, and it is around 5 °C in warm season.   Relevant parameters of the one-dimensional model are illustrated in Table 1. It is assumed that soils with the same lithology have the same ranges of shear moduli and damping ratios, but have different densities and different shear wave velocities at various depths [28]. Shear wave velocities of soil layers at different depths, shown in Table 1, are determined according to the field and laboratory shear wave velocity tests. Then the seismic site response is calculated by using the equivalent linear method.

Dynamic parameters
The shear modulus and damping ratios are two most important parameters in dynamic analysis, and variation of these parameters with shear strain amplitude can then be tested in the laboratory. A series of dynamic triaxial tests of frozen soil was performed at the State Key Laboratory of Frozen Soil Engineering, Chinese Academy of Sciences. Then the normalized shear modulus G m / G max , as well as the damping ratios, were determined. The variation of dynamic parameters of the typical frozen soil are shown in Figure 6., the dynamic shear modulus ratio G m /G max decreases with an increase in shear strain γ and damping ratios λ.

Seismic input motions
Ground motion characteristics should be taken into account when site response analyses are conducted. Considering the seismic intensity distribution and ground motion parameters in the QTP, artificial seismic waves, including three different seismic intensity levels, are used in this paper as inputs to the bedrock [29].
Three intensity levels are defined as exceedance probabilities, i.e. 63 % in 50 years, 10 % in 50 years, and 5 % in 50 years. Figure 7 shows the input acceleration time histories and their Fourier spectra characteristics. Peak accelerations of 26 cm/s 2 , 167 cm·s −2 and 326 cm/s 2 refer to the small earthquake (a), medium earthquake (b) and large earthquake (c), respectively.

Seismic site response characteristics
Ground motion acceleration time histories at Beiluhe section of QTR, under seismic motion with exceedance probabilities of 10.0 % in 50 years, are illustrated in Figure 8(a), while spectrum characteristics of these waveforms are illustrated in Figure  8(b). It can be intuitively seen that the phase difference of the time history is relatively small, and spectral characteristics are almost the same under different thermal conditions. The peak ground acceleration (PGA) is 198.35 cm/s 2 when the active layer is 5°C, while this value decreases to 187.50 cm/s 2 at -1°C. The PGA increases when the temperature of the active layer increases, which can be expected due to energy trapping and dissipation in a thermally softened layer. The normalized acceleration response spectra are shown in Figure  9. It can be seen that the DAF (dynamic amplification factor) of the response spectra in cold regions is about 3.0. This index denotes the ratio of the maximum acceleration of the mass The data shown in Figure 9 reveal that the DAF is almost the same when the active layer has different temperatures. It can be seen that the variation in thermal conditions of the active layer does not change inherent characteristics of the site. The stiffness variation of active layer may have little effect on the response spectra. is presented in this paper in order to quantify and estimate this effect. This coefficient is defined as the ratio of peak ground acceleration (PGA) to input seismic acceleration (ISA). The distribution of the PGA Camp at ground surface and the permafrost table with the active layer temperature are shown in Figure 10, considering the artificial seismic waves with different exceedance probabilities in 50 years. It has been established that the peak ground acceleration is the greatest when the active layer thaws in warm seasons. The peak acceleration at the permafrost table has nearly the same values, which indicates that the thermal state of the active layer mainly influences the dynamic response of the shallow soil above the permafrost table in permafrost regions. Figure 10 probably represents nonlinearity of the PGA camp. The intensity of seismic motion has a significant influence on the seismic amplification effect. The PGA Camp decreases from 1.3 to 1.1 with an increase in the input seismic acceleration from 26 cm/s 2 to 326 cm/s 2 . While at higher input seismic acceleration, it is necessary to stress that the PGA Camp may further decrease, and even have no amplification effect.

Seismic response of embankment
According to the above analysis and discussion, the seismic site response of permafrost regions is strongest in warm seasons. The seismic site amplification effects are also studied. Thus, on this basis, by applying the nonlinear dynamic finite element analysis method, the embankment in warm season is selected as a typical case, and the seismic response is studied with the ABAQUS software.

Numerical model and parameters
The numerical finite element model of the typical traditional railway embankment at Beiluhe section of the QTR is built on the basis of the field geology survey and embankment structure. The numerical analysis model is shown in Figure 11. This finite element model comprises five layers (the ballast is one layer) which are consistent with the stratigraphy shown in Figure 1. In the numerical case, the roadbed filling and the active layer is thawed and the underlying soil layer is permafrost. Various densities and mechanical parameters, corresponding to the one-dimensional model, are also taken into account.
In this simulation, numerical calculation is deduced based on the plane strain assumption, and a typical dimension of the elements is approximately 1×1 m.  In the dynamic finite element analysis, the soil is regarded as an elastic-plastic material and its shear modulus (Gmax) can be determined according to the shear wave velocity (V s ) and the mass density (ρ) of soil by the following formula (1): (1) In the above mentioned formula, the shear wave velocity (VS) and the mass density (ρ) are defined as related to soil parameters of the one-dimensional model in Table 1. Moreover, an elastic perfectly plastic constitutive model with a Mohr-Coulomb failure criterion is employed [30]. It has been proven that the Mohr-Coulomb criterion is suitable for frozen soil, as witnessed by a large number of mechanical tests of frozen soil [31][32][33]. Physical and mechanical soil properties in numerical calculation, as shown in Table 2, are determined by dynamic triaxial tests of frozen soil [34][35][36][37][38]. Figure 12 illustrates contour plots of horizontal acceleration at the time of peak acceleration, under excitation of a large earthquake with 2 % probability of exceedance in 50 years. It is indicated that the acceleration response of the embankment is greater than the ground surface. The energy accumulates in the embankment due to the reflection and refraction of seismic wave during the wave propagation in the embankment. In this paper, the results of the seismic site response analysis are considered acceptable. Then the 2D simulation results are compared with the 1D model results. In this study, the acceleration time-history curve of the observation point G on the natural ground surface is extracted, and the maximum value is compared with the PGA obtained from the seismic site response. Figure 13 shows the ground motion acceleration time histories of observation points G and S under the seismic motion with different probabilities of exceedance in 50 years. Through the comparison, it can be seen that the 2D simulation results are quite similar to 1D model results. Moreover, when comparing the observation points, the maximum acceleration values at embankment shoulder are larger than those registered at the natural ground surface. The distribution of PGA amplification coefficients at different observation points with the exceedance probability is presented in Figure 14. The peak ground acceleration amplification coefficient on the natural ground surface is 1.5 during a small  In the case of medium and large earthquakes, the peak ground acceleration amplification coefficient on the natural ground surface almost maintains its value, marking a 58 % increase as well.

Numerical results and analyses
The amplification coefficient will decrease with an increase in intensity. When the intensity reaches a certain level, this coefficient will almost maintain its value.  Numerical investigation of seismic behaviour of railway embankments in cold regions The analysis of the above mentioned dynamic properties such as acceleration, velocity and displacement shows that dynamic response of the embankment has distinct nonlinear characteristics. Figure 17 shows the plastic strain magnitude of the embankment at the time of peak acceleration under seismic excitation, with 10 % and 2 % probabilities of exceedance in 50 years. When the seismic intensity reaches a certain value, the plastic zone gradually appears in the embankment, and the plastic zone extends continuously with an increase in earthquake intensity. It can be seen that the soil at embankment toe enters the plastic state, the plastic strain magnitude reaches 5·10 -3 at the time of peak acceleration under the seismic motion, with 10 % probabilities of exceedance in 50 years, while this value reaches 1.9·10 -2 under the seismic motion with 2 % probabilities of exceedance in 50 years as well.

Discussion of results
The dynamic response of the shallow soil above the permafrost table was mainly influenced by the thermal state of the active layer. The dynamic response of permafrost is inhibited with the freezing of the active layer, mainly because of the change of dynamic properties of soil during the soil freezing-thawing process. Meanwhile, the dynamic changes of the soil mechanics parameters with seasons may not be sufficient to alter the site property in permafrost regions, expressing that the site response spectra are almost the same with different active layer thermal state. In addition, the seismic site effects are strongly related to the seismic wave itself, such as the near field and farfield earthquake accelerograms which have different frequency and amplitude characteristics. This problem still needs to be discussed in greater detail. The above mentioned seismic response analysis of embankment shows that the embankment has higher acceleration than the natural ground surface, which is due to the effect of energy accumulation in the embankment. Moreover, the velocity and displacement of the embankment have distinct nonlinear characteristics with an increase in amplitude of seismic motion. Because of the damage of soil structure, the soil in the embankment enters the plastic state, and then the subgrade slope reinforcement should be considered in practical engineering. The dynamic stability of the embankment under large earthquakes is another issue that needs to be thoroughly studied.
It should be noted that one specific site situated in Beiluhe segment of the Qinghai-Tibet Railway (QTR) was selected and that the corresponding studies were conducted in this paper, which is why the results can not be generalized to all sites and all soils in the region.

Conclusion
The seismic site response in permafrost sites is analysed in this work, considering the thermal state of the active layer. On this basis, the seismic response of embankment in the warm season is further studied by applying the nonlinear dynamic finite element analysis method. The following conclusions can be drawn from the results of this study: The active layer may have little influence on the seismic site effect in permafrost regions. The thermal state of the active layer mainly influences the dynamic response of the shallow soil above the permafrost table, and thermal state of the active layer has little influence on the response spectra. The response of the embankment is much greater than that of the ground. The peak ground acceleration at embankment shoulder is larger than that at the natural ground surface, marking a 73 % increase compared to the coefficient at the natural ground surface. The dynamic response of the embankment has distinct nonlinear characteristics. When the seismic intensity reaches a certain value, the plastic zone gradually appears in the embankment, and there is a continuous extension of the plastic zone with an increase in earthquake intensity. Because of the damage to soil structure, the soil in the embankment enters the plastic state, and then the reinforcement of the embankment toe should be considered in practical engineering.