Assessing the Risk Probability of the Embankment Seismic Damage Using Monte Carlo Method

School of Civil and Architecture Engineering, Shandong University of Technology, Zibo, China Zibo Transportation Service Center, Zibo, China Key Laboratory of Roads and Railway Engineering Safety Control (Shijiazhuang Tiedao University), Ministry of Education, Shijiazhuang, China State Key Laboratory of Mechanical Behavior and System Safety of Traffic Engineering Structures (Shijiazhuang Tiedao University), Shijiazhuang, China Laoling Branch of Dezhou Highway Development Center, Dezhou, China


Introduction
Embankment is one of the most common structural forms of highway subgrade [1,2]. In Tangshan earthquake, Wenchuan earthquake, Kyushu earthquake, East Japan earthquake, and other previous violent earthquakes, embankments suffered from varying degrees of damages, which seriously disrupted the highway networks [3][4][5][6]. It is greatly significant for improving the highway antiseismic capacity through proposing proper embankment reinforcement method on the basis of assessing the risk probability of the seismic damage [7][8][9]. e risk probability assessment is based on the seismic fragility assessment and seismic hazard assessment and analyzes their coupling characteristics. e assessment results can be represented by the risk probability exceeding each seismic damage level for a certain embankment in a future period of time [10][11][12][13][14].
Due to the increasing in frequencies and losses of earthquakes, researches on the risk probability assessment of the embankment seismic damage have attracted attention from many scholars [15][16][17][18][19]. Enomoto and Sasaki [20] conducted a series of dynamic centrifuge model tests in order to evaluate some factors affecting the seismic performance of hillside embankments consisting of sandy or silty soils and resting on stiff base slope. Hutt et al. [21] benchmarked the antiseismic performance of a 50-story archetype office building designed following "Uniform Building Code (1973)" against "International Building Code (2015)" based on the seismic risk assessment. e mean annual frequency collapse risk of the 1973 building was 28 times greater than the equivalent 2015 building or approximately 13% versus 0.5% probability of collapse in 50 years. Liu et al. [22] performed the seismic fragility analysis of recycled aggregate concrete (RAC) bridge columns with different recycled coarse aggregate (RCA) replacement ratios subjected to freeze-thaw cycles (FTCs) by the cloud analysis method. Xu et al. [23] developed the Cumulative Absolute Velocity (CAV) with 24667 strong ground-motion records from Taiwan and the seismic hazard assessment for Taipei was then conducted using the local CAV models. Zanini et al. [24] proposed a framework for the development of seismic risk maps at nationwide scale. e seismic risk map was computed taking into account a seismogenic model of the analyzed area and properly characterizing vulnerability and exposure of an asset of interest. Scozzese et al. [25] analyzed the effectiveness of damping and isolation devices in mitigating the seismic risk by proposing a general framework for investigating the sensitivity of the seismic risk of structural systems with respect to system properties varying in a prescribed range. Bartoli et al. [26] analyzed the seismic risk of a confined historic masonry tower located in Italy and critically discussed the multilevel assessment path proposed by the "Italian Guidelines for the Assessment and Mitigation of the Seismic Risk of the Cultural Heritage." Risia et al. [27] studied how different procedures of site response analysis influence the risk quantification at urban scale and the risk was quantified in terms of direct losses incurred by the portfolio of buildings in Benevento for the 1980 Mw 6.9 Irpinia earthquake. Andreotti and Lai [28] proposed a comprehensive numerical methodology to construct fragility curves for mountain tunnels subjected to transversal seismic loading. e proposed numerical technique, which was based on fully nonlinear dynamic analyses accounting for the nonlinearity of both ground and tunnel support, allowed considering the features of the tunnel. Wang et al. [29] solved time-evolving probabilistic structural response through developed stochastic elastic-plastic finite element method (FEM); following that, formulation for seismic risk analysis was derived. Liu et al. [30] presented a novel perspective on scenario-based seismic vulnerability and hazard assessment; the findings could be a potential guide to decision-making in disaster risk reduction in rural Weinan (China). Maio et al. [31] aimed to investigate the correlation between two well-known approaches for the seismic risk assessment of stone masonry buildings located within historic centers: the vulnerability index method and nonlinear static seismic analyses. e latter were carried out by using a new three-dimensional macroelement model to numerically represent the considered sample of prototype buildings, together with the application of the N2 method. Castaldo and Alfano [32] developed incremental dynamic analyses to assess the seismic fragility considering several natural ground motions and different elastic and inelastic structural properties under the hypothesis of modeling the friction coefficients of the two surfaces of the double friction pendulum system (DFPS) as random variables.
ere are few systematic studies on the hazard, fragility, and risk probability assessment, although the embankment seismic damages are extremely serious [11,33]. Meanwhile, retaining walls and other retaining structures play positive roles in improving the embankment seismic performance, but relevant quantitative study is rarely reported [34]. In particular, there are few existing studies on the seismic damage characteristics of the embankment fill-soil foundation system and retaining structure-embankment fill-soil foundation system under different seismic actions [35]. In view of this, the Xi'an-Baoji expressway K1125 + 470 embankment was taken as the research object to assess the risk probability of the seismic damage and the positive influences of the retaining wall on reducing the risk probabilities were verified.

Methodology.
According to the differences in data acquisition and calculation method, the seismic fragility assessment is divided into empirical method and theoretical method [36,37]. e theoretical method reflects the probabilistic relation between the seismic damage level and corresponding seismic ground-motion intensity through IDA and PSDA analysis [38,39]. is method was adopted to assess the seismic fragility of the research object; the reasons are as follows [40][41][42].

Easy to Control Loading Levels.
is method applies computer loading and can select different types and intensities of the seismic ground-motion records, thus completely reflecting the uncertainties of the ground-motions and reducing the influences of the subjective factors on the assessment results.

Realizing Quantitative Control of the Seismic Damage.
Seismic damage levels can be output through software, avoiding the errors in judging the seismic damage levels by manual surveys, thus realizing quantitative control of the seismic damage.

Realizing Accurate Study on a Monomer Structure.
Finite difference model of a monomer structure can be established by computer automatically, thus completely reflecting the antiseismic capacity of different structures and realizing accurate study on a monomer structure.

Data Preparation.
By referencing on Castaldo et al. [43], the finite difference model of the embankment fill-soil foundation system of the research object was established by adopting Flac software, wherein the vehicle loads were converted into the thickness of the embankment fill according to the elastic lamination theory [44,45], as shown in Figure 1.
In modeling the embankment fill-soil foundation system, an elastoplastic constitutive relationship was employed, and the Mohr-Coulomb criterion was used as the yield criterion [43]. e mechanical properties of the embankment fill and soil foundation were determined as follows: shear modulus was 17.91 MPa and 15.67 MPa, respectively, density was 1970.00 kg/m 3 and 1630.00 kg/m 3 , respectively, elastic modulus was 48.00 MPa and 42.00 MPa, respectively, Poisson's ratio was 0.34 and 0.36, respectively, bulk modulus was 50.00 MPa and 43.75 MPa, respectively, internal friction angle was 33.00°and 28.00°, respectively, and cohesive force was 34.00 KPa and 31.00 KPa, respectively [46,47]. e drained conditions were not assumed.
Under the action of the ground acceleration € u g , the fundamental motion equation of the embankment fill-soil foundation system is shown in the following [48]: where M refers to the total mass matrix containing the added vehicle mass, C refers to the total damping matrix, K refers to the total stiffness matrix, J refers to the indicator matrix of each seismic component, and € u, _ u, and u refer to the acceleration array, speed array, and displacement array of the nodes, respectively [49,50].
e Rayleigh damping was selected as the model damping, PGA was selected as the seismic intensity parameter, and ε max was selected as the seismic damage parameter, as shown in the following: where d max is the maximum lateral permanent displacement on the embankment surface and D is the embankment width. e reasons for selecting ε max are as follows: (1) ε max cannot only characterize the damage to the whole embankment, but also characterize the most vulnerable parts; (2) ε max can reflect the seismic effect on traffic capacity: the larger the ε max the more serious seismic damage and, therefore, the lower the traffic capacity. e corresponding relationship between the seismic damage level and ε max was established on the basis of the displacement damage criterion, as shown in Table 1. 10 seismic ground-motion records provided by PEER for IDA analysis were selected [51]. e epicentral distance (E d ) of these records ranges from 11.7 km to 50.9 km, the magnitude (M w ) ranges from 6.0 to 7.6, and the original PGA ranges from 0.094 g to 0.714 g, as shown in Table 2.
Due to the limited space, the time-history curves of acceleration of the No. 1-No. 3 seismic ground-motion records are listed in Figure 2.

Embankment Seismic Fragility Assessment Results.
e 120 seismic ground-motion records were input to the finite difference model shown in Figure 1, and 120 dynamic response analysis was undertaken [53,54]. ε max of each analysis as well as their mean value and standard deviation was computed, as shown in Table 3. e total energy dissipation of each analysis was also computed, as shown in Figure 3. e logarithms of the PGA and ε max were taken according to the dynamic response analysis results of the embankment fill-soil foundation system so as to obtain the probabilistic seismic demand model [32,55,56], as shown in the following: ln ε max � ln 1.1396 + 1.1821 ln PGA. (3) Formula (3) was substituted into the classical calculation formula of the seismic fragility to obtain the following: where P j refers to the exceeding probability of the embankment seismic damage level j and j � 1 represents basically intact, j � 2 represents slight damage, j � 3 represents moderate damage, j � 4 represents severe damage, j � 5 represents destruction, and S j refers to the structural performance level shown in Table 1, namely, S 2 � 0.20 when j � 2, S 3 � 0.40 when j � 3, S 4 � 0.80 when j � 4, and S 5 � 1.20 when j � 5. e embankment seismic fragility curves were plotted according to formula (4), as shown in Figure 4. e exceeding probability of each seismic damage level of the embankment fill-soil foundation system was calculated according to Figure 4, as shown in Table 4. e exceeding probabilities of the severe damage and destruction are 0.760416 and 0.458867 when PGA � 1.0 g, respectively, which are basically consistent with the survey results of the embankment seismic damage near the epicenter of the Wenchuan earthquake (PGA � 0.977 g) and further show that the embankment seismic fragility assessment method adopted in this paper is scientific and rational.

Risk Probability Assessment of the Embankment Seismic Damage
3.1. Methodology. Risk probability of the embankment seismic damage numerically equals the convolution between the seismic hazard of the site and the embankment seismic fragility, that is, the integral of exceeding probabilities of all damage levels under all seismic events in the next t years [57], as shown in the following: where Pj,t refers to the risk probability exceeding the seismic damage level j in the next t years and ft(S) refers to the probability density function of the comprehensive effect on the embankment from the risk events, that is, the probability Advances in Civil Engineering density function of the values of ε max in the next t years [58]. e calculation method of ft(S) is shown in the following: where E refers to the earthquake event and f(S E) refers to the conditional probability density of ε max under any given E, that is, the results of the embankment seismic fragility assessment. Ft(E) refers to the probability density function of the seismic event, that is, PGA distribution model of the embankment site, being results of the seismic hazard assessment [59]. According to (5) and (6), the calculation method of Pj,t is shown in the following:

Seismic Ground-Motion Distribution Model.
Seismic ground-motion distribution model includes the probability distribution model of the site seismic intensity and transformational relationship between the seismic intensity and PGA.

Probability Distribution Model of the Site Seismic
Intensity. Liu and Zhang [60] analyzed the seismic hazard    assessment results of dozens of sites in the " ree North Areas" in China. According to the results, the seismic intensities of a certain site in the next 50 years comply with the extreme value type (III) distribution; the distribution function is shown in the following: where ω refers to the maximum seismic intensity value. According to China's current seismic intensity   Advances in Civil Engineering 5 classification method, ω � 12. ε refers to the normal seismic intensity, that is, the seismic intensity with exceeding probability of 63% in the next 50 years. K refers to the shape parameter, which can be determined by 3 methods, that is, fractile method, least square method, and maximum likelihood method [61]. From the view of accuracy of fitting, the least square method is optimal; however, the K value corresponding to the basic seismic intensity (exceeding probability of 10% in the next 50 years) can satisfy the requirements from the view of engineering practice [62]. According to (8), the probability distribution function of the site seismic intensity in the next t years is shown in the following: Taking the derivative of (9), the probability density function of the site seismic intensity in the next t years was obtained, as shown in the following:

Transformational Relationship between the Seismic
Intensity and PGA. e transformational relationship between the seismic intensity and PGA proposed by Liu and Zhang [60] was adopted in this paper, as shown in the following: where I refers to the seismic intensity. Although the seismic intensity values are integer discrete variables between 1 and 12, in consideration of the calculation accuracy and model applicability, the seismic intensity values in (11) are considered as continuous variables.

Risk Probability Assessment Based on the Monte Carlo
Method. According to (3), (10), and (11), complicated integral operations will be inevitably conducted when the numerical method is directly applied to calculate the risk probability, and convergence in the field of integral cannot be guaranteed. erefore, the Monte Carlo method was applied to assess the risk probability of the embankment seismic damage.

Basic Idea.
MATLAB software was used to conduct n times of sampling for the seismic intensity based on (10). e seismic intensity and PGA were transformed based on (11); then, the occurrence probability of an individual PGA in overall samples was 1/n. Pj(PGAi) was defined as the risk probability exceeding the seismic damage level j under the action of PGAi, which can be calculated based on the fragility assessment results. When n is large enough, (7) can be changed to (12) according to the law of large numbers [63,64]. where

Case Study.
According to the seismic hazard assessment results along the Xi'an-Baoji expressway, PGAs with certain exceeding probabilities in the next 50 years were calculated [11]. According to (11), the corresponding seismic intensities were calculated, as shown in Table 5. According to Table 5 and Figure 5, the shape parameter K was determined based on the seismic intensities with exceeding probabilities of 10% and 63% in the next 50 years, as shown in the following: e solution is K � 8.397. e probability distribution of the seismic intensity of the embankment site in the next 50 years is shown in the following formula and Figure 6.
12 − x 5.615 8.397 . (14) Taking the derivative of Formula (14), the probability density of the seismic intensity of the embankment site in the next 50 years is shown in the following formula and Figure 7. . (15) In consideration of the calculation accuracy and model applicability, the best results can be obtained when the number of simulations is 100000 when the Monte Carlo method is conducted [65]. According to the probability density of the seismic intensity of the embankment site shown in Formula (15), the MATLAB software was used in sampling and 100000 random numbers I 1 , I 2 , I 3 , . . .. . ., I 100000 were selected, which indicated the 100000 possibilities of the seismic intensities; besides, the above random numbers complied with the seismic hazard assessment results. After inputting I 1 , I 2 , I 3 , . . .. . ., I 100000 to Formula (11), 100000 possible PGAs of the embankment site in the next 50 years PGA 1 , PGA 2 , PGA 3 , . . .. . ., PGA 100000 were calculated.

Risk Probability Assessment Results of the Retaining Wall-Embankment Fill-Soil Foundation System.
In order to verify the positive influences of the retaining wall on the risk probability of the embankment seismic damage, a retaining wall was assumed to be built at the left side of the research object. An isotropic elastic constitutive relationship was used to model the retaining wall, and the mechanical properties were: elastic modulus 3000.00 MPa, Poisson's ratio 0.17, bulk modulus 1515.15 MPa, shear modulus 1282.05 MPa, and density 2300.00 kg/m 3 (Yin et al. [11]), and the structural form of the retaining wall-embankment fill-soil foundation system is shown in Figure 5. 120 dynamic response analysis were conducted on the retaining wall-embankment fill-soil foundation system. e total energy dissipation of each analysis is shown in Figure 8, the fragility curves are shown in Figure 9, the assessment results are shown in Tables 7 and 8.
It can be seen that the risk probability of the embankment seismic damage can be significantly reduced through the retaining wall; for example, the risk probabilities of the severe damage and destruction of the retaining wall-embankment fill-soil foundation system are 29.07% and 7.62%, while risk probabilities of the embankment fill-soil foundation system are 19.30% and 3.46% and there is a reduction of 33.61% and 54.60%, respectively. Figure 10 lists the seismic damages of the embankment fill-soil foundation system and the retaining wall-embankment fill-soil foundation system under the actions of No. 1 seismic ground-motions when PGA � 0.2 g, 0.4 g, 0.6 g, 0.8 g, 1.0 g, and 1.2 g, respectively.

Influences of the Retaining Wall on the Embankment Seismic Damages.
It can be seen from Figure 10 that the seismic damages mainly include cracks and lateral displacement on the top surface of the embankment fill, wherein the damages of the retaining wall-embankment fill-soil foundation system are obviously lower than those of the embankment fill-soil foundation system; namely, the ε max is reduced by 18.11% and the crack length is shortened by 21.83% on average. However, the seismic damages of the retaining wall itself are also very serious; especially after a violent earthquake, the  engineering disturbance of the retaining wall is reduced greatly or even lost completely. Figure 11 lists the seismic damage evolution processes of the embankment fill-soil foundation system and the retaining wall-embankment fill-soil foundation system under the action of No. 1 seismic ground-motion when PGA � 1.2 g, it can be seen that the differences of the seismic damages between the two systems are not obvious at the initial stage of a violent earthquake (t � 2 s and 4 s). e seismic damage of the retaining wall-embankment fill-soil foundation system is obviously lower than that of the embankment fill-soil foundation system since the retaining wall bears larger passive Earth pressure at the middle stage of a violent earthquake (t � 6 s and 8 s). At the final stage of a violent earthquake (t � 10 s and 12 s), the seismic damage of the retaining wall is increased and the passive Earth pressure is reduced with continuous development of the cracks and lateral displacement on the top surface of the embankment fill, and the seismic damage increasing range of the retaining wall-embankment fill-soil foundation system is basically consistent with that of the embankment fill-soil foundation system.      Figure 11: Seismic damage evolution processes. (a) Embankment fill-soil foundation system, t � 2 s, ε max � 0.0927, and crack length is 1.24 m. (b) Retaining wall-embankment fill-soil foundation system, t � 2 s, ε max � 0.0854, and crack length is 1.19 m. (c) Embankment fill-soil foundation system, t � 4 s, ε max � 0.3264, and crack length is 5.37 m. (d) Retaining wall-embankment fill-soil foundation system, t � 4 s, ε max � 0.3004, and crack length is 5.12 m. (e) Embankment fill-soil foundation system, t � 6 s, ε max � 0.5846, and crack length is 9.67 m. (f ) Retaining wall-embankment fill-soil foundation system, t � 6 s, ε max � 0.4237, and crack length is 7.51 m. (g) Embankment fill-soil foundation system, t � 8 s, ε max � 0.8064, and crack length is 15.85 m. (h) Retaining wall-embankment fill-soil foundation system, t � 8 s, ε max � 0.5932, and crack length is 11.24 m. (i) Embankment fill-soil foundation system, t � 10 s, ε max � 1.2982, and crack length is 20.37 m. (j) Retaining wall-embankment fill-soil foundation system, t � 10 s, ε max � 0.9567, and crack length is 16.58 m. (k) Embankment fill-soil foundation system, t � 12 s, ε max � 1.5843, and crack length is 24.84 m. (l) Retaining wall-embankment fill-soil foundation system, t � 12 s, ε max � 1.3548, and crack length is 20.33 m.

Conclusions
(1) e Xi'an-Baoji expressway K1125 + 470 embankment was taken as the research object, the risk probability assessment of the embankment seismic damage was conducted on the basis of the seismic fragility and seismic hazard assessment, the positive role of the retaining wall in improving the embankment seismic performance was verified, and the seismic damage evolution processes of the retaining wall-embankment fill-soil foundation system and the embankment fill-soil foundation system were analyzed. (2) e risk probabilities of the severe damage and destruction of the embankment fill-soil foundation system are 29.07% and 7.62% in the next 50 years, while those of the retaining wall-embankment fillsoil foundation system are 19.30% and 3.46% and are reduced by 33.61% and 54.60%, respectively. e influence of the retaining wall on the embankment seismic performance mainly occurs at the middle stage of a violent earthquake, but the seismic damages of the retaining wall itself are also very serious, and the engineering disturbance of the retaining wall is reduced greatly, or even lost completely at the final stage of the violent earthquake.
(3) e research results can be used as the theoretical basis of the embankment antiseismic design, and the proposed risk probability assessment method can also be applied to similar structures. However, some aspects of this paper can be further improved; for example, the study on the dynamic coupling mechanism of the embankment and the seismic ground-motion can continue to be developed in the future.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest with respect to the research, authorship, and publication of this article.