Embankment Seismic Fragility Assessment under the Near-Fault Pulse-like Ground Motions by Applying the Response Surface Method

School of Civil and Architecture Engineering, Shandong University of Technology, Zibo 255049, China Zibo Transportation Service Center, Zibo 255000, China State Key Laboratory of Mechanical Behavior and System Safety of Traffic Engineering Structures, Shijiazhuang Tiedao University, Shijiazhuang 050043, China Key Laboratory of Roads and Railway Engineering Safety Control (Shijiazhuang Tiedao University), Ministry of Education, Shijiazhuang 050043, China School of Transportation and Vehicle Engineering, Shandong University of Technology, Zibo 255049, China Shandong Dongtai Engineering Consulting Co., Ltd., Zibo 255000, China Zibo International Academy at High-Tech Zone, Zibo 255000, China Laoling Branch of Dezhou Highway Development Center, Dezhou 253000, China


Introduction
Strong earthquakes have occurred frequently in the world; for example, the 2004 Sumatra earthquake in Indonesia, the 2008 Wenchuan earthquake in China, the 2010 earthquakes in Haiti and Chile, the 2011 east Japan earthquake, and the 2017 Jiuzhaigou Valley's earthquake in China resulted in severe structural damages, causing major human casualties and economic losses [1][2][3]. Near-fault ground motions generally occur at sites located near active faults and are characterized by the following main features: (1) large amplitude and long-period velocity pulses; (2) high Peak Ground Velocity/Peak Ground Acceleration (PGV/PGA) and Peak Ground Displacement/Peak Ground Acceleration (PGD/PGA) ratios; and (3) concentration of energy in one or few pulses. e pulse-like ground motions are caused primarily by the forward directivity effect, which are observed at a site when the fault rupture propagates towards the site with a velocity close to shear wave velocity [4][5][6]. e near-fault ground motions with long-period pulses have much more energies input to long-span linear structures than the far-fault ground motions. erefore, researches on the fragility assessment of structures such as embankments under the near-fault pulse-like ground motions are of great significance for performance-based seismic design [7][8][9].
Uncertainties of the ground motions and structural parameters are the main factors that limit the accuracy of the seismic fragility assessment [10,11]. e uncertainties of the ground motions are influenced by aleatory uncertainties (uncertainty due to the random nature of the processes under consideration such as ground motion record-to-record variability) and epistemic uncertainties (uncertainty due to incomplete knowledge and data such as model parameter uncertainty, omissions, and errors) [12][13][14][15]. ere are two main methods to solve the problem of the uncertainties of the ground motions. e first selects the actual recorded ground motions provided by the Pacific Earthquake Engineering Research Center (PEER). For example, Yin et al. [16] selected 15 ground motions in medium-hard sites to carry out the embankment seismic fragility assessment and verified the impact of the retaining wall on improving the embankment seismic performance. Bao et al. [17] selected 10 near-fault ground motions and 10 far-fault ground motions to carry out the fragility assessment of the nuclear reactor shell structure under the main shock and aftershocks, respectively. Li et al. [7] investigated the dynamic response characteristics of the 1/30 scale model of the Wuhan Shimao Center subjected to the El Centro wave and Taft wave by analyzing the shaking table test data. He et al. [18] adopted 22 actual recorded ground motions according to the recommendations of the Applied Technology Council (ATC) in the FEMA P695 report in the seismic fragility analysis of multiage buried steel pipes. Since it is difficult to find the ground motions that fully comply with the lithology properties of where the structure is located and with the law of seismic wave propagation, the accuracy of the seismic fragility assessment based on actual recorded ground motions is relatively low [19,20]. e second is synthesizing ground motions. For example, Sheng et al. [21] calculated the fragility indexes of the frequent earthquake, fortification earthquake, and rare earthquake of group structures combined with the synthesized ground motions and empirical earthquake damage indexes. Salami et al. [22] carried out the seismic fragility assessment of low-rise reinforced concrete structures under the main shock and aftershocks based on the synthesized ground motions. Abyani et al. [23] used a genetic algorithm to analyze the impact of the number of synthesized ground motions on the fragility assessment results of jacket offshore platforms. For near-fault pulse-like ground motions, the maximum values of the parameters in different directions vary significantly. For example, Yang and Zhou [24] studied the impacts of the pulse-like ground motions on the response spectrum in different directions and pointed out that the ground motions on the strongest pulse direction comprehensively considered the impacts of the pulse period, peak velocity, and duration, and the energies input to the structure were maximized. However, existing artificial near-fault ground motion models do not take into account the directivity of pulse-like ground motions in history analysis, and the seismic performance may be overestimated in the seismic fragility assessment [15,[25][26][27].
Some scholars ignore the impacts of the uncertainties of the structural parameters on the seismic performance and build models using deterministic values directly, while some introduce the response surface method to characterize the mapping relationships among the structural parameters and seismic responses [28].
ere are 3 steps of the seismic fragility assessment based on the response surface method: (1) experimental design, that is to determine the sample points and structural response values [29]; (2) select the response surface function that meets the structural functional requirements [30]; and (3) fit the response surface function and verify the model accuracy [31,32]. For example, Saha et al. [33] used the Monte Carlo method to study the discreteness of the parameters of the tank structure and carried out the seismic fragility assessment based on the response surface method. Li and Li [29] proposed a uniform design response surface method to study the impacts of the multiple correlations of the design variables on the bridge seismic fragility. Tran et al. [34] carried out the seismic fragility assessment of nuclear power plants based on the response surface method and plotted the fragility curves by the maximum likelihood estimation method and linear regression method, respectively. Since the response surface functions often take the form of summation of several complete polynomials, when there are too many structural parameters, the applicability of the model greatly reduces and the solution of the response surface function tends to fall into a "curse of dimensionality" [30,[35][36][37].
In view of this, the Xi'an-Baoji expressway K1125 + 470 embankment was taken as the research object, the stochastic pulse models with different rupture fault distances were established by considering the directional impacts of the pulse-like ground motions, and the linear correlation of the structural parameters was eliminated based on the principal component analysis. In addition, embankment seismic fragility assessment using 15 artificial pulse-like ground motions with the rupture fault distances ranging from 1 to 15 km was carried out by applying the response surface method, aiming to provide a theoretical basis for the performance-based design of embankment.
wavelet transformation results of two horizontal components (fault-normal velocity history and fault-parallel velocity history) of the pulse-like ground motions was used to build the motion model on the strongest pulse direction, and the following methods were applied for stochastic modeling and artificial synthesizing of the near-fault pulse-like ground motions [24,26,[39][40][41]: (1) Apply wavelet transformation to synthesize the faultnormal velocity history and fault-parallel velocity history of the pulse-like ground motions into the velocity history on the strongest pulse direction and divide the velocity history into low-frequency components and high-frequency components.
(2) Fit the low-frequency components using the Gabor wavelet analytic function and propose the correlations among the pulse parameters, i.e., seismic moment, rupture fault distance and site condition, and seismological parameters. e regression model between the PGV on the strongest pulse direction and rupture fault distance is calculated, and the lowfrequency stochastic pulse function is established after obtaining the probability distribution of the low-frequency pulse parameters.
(3) Differentiate the residual velocity history containing the high-frequency components to obtain the residual acceleration history and employ the improved Kanai-Tajimi power spectrum model to fit the statistical power spectrum containing the high-frequency components to obtain the power spectrum parameters.
(4) Use the envelope function to fit the residual acceleration history after Hilbert transformation to reflect the nonstationarity of the high-frequency components and obtain the probability distribution of the corresponding envelope parameters via statistics.
(5) Generate the velocity history of the near-fault pulselike ground motions by superimposing the longperiod velocity pulse and velocity history with the high-frequency components and obtain the acceleration history through differentiation, as shown in Figure 1.
e pulse-like ground motions of the station CHY080 of the Chi-Chi earthquake on September 21, 1999, in Taiwan, China, were adopted to synthesize 15 ground motions with the rupture fault distances ranging from 1 to 15 km, respectively. e rupture fault distance of the CHY080 station was 10.9 km, and the PGA was 1.193 g [42][43][44].
e regression model between the PGV on the strongest pulse direction and rupture fault distance was calculated according to the recorded ground motions, as shown in Figure 2.
For the synthesized ground motions were generated from the same earthquake event, i.e., Chi-Chi earthquake, the magnitude M w � 7.6 was not included in the regression model. e standard deviation σ 1 reflects epistemic modeling uncertainty regarding the factors controlling variety of PGVs with the same rupture fault distance R and follows a log-normal distribution. A representative sample under the scenario of rupture fault distance R � 15 km was used to generate the synthesized ground motions according to Figure 1, as shown in Figure 3. It is indicated that the synthetic velocity history represents the pulse characteristic of the near-fault ground motions properly, and the corresponding acceleration history can be used as the stochastic excitation input for the embankment seismic fragility assessment. It is worth noting that the seismic fragility assessment in this paper is only based on the Chi-Chi earthquake, and the results are valid for earthquakes with the same characteristics as the Chi-Chi earthquake.

Response Surface Method
3.1. Embankment Model. FLAC can simulate structures made of soil, rock, and other materials that flow plastically when reaching their yield limits [45]. FLAC adopts the finite difference scheme to solve the governing differential equation, which can accurately simulate the yield, plastic flow, softening, and large deformation of materials, especially has unique advantages in the fields of elastic-plastic analysis, large deformation analysis, and simulation of construction process [46]. A finite difference model of the embankment fill-soil foundation system of the research object was established by using FLAC. e width of the embankment fill was 24.5 m, the left side was 6.1 m high, the right side was 2.6 m high, and the slope ratio was 1:1.5. e slope of the soil foundation was 24°, the thickness was 30 m, and the width was 120 m. Among them, the vehicle loads have been converted to the thickness of the embankment fill according to the elastic layer theory, as shown in Figure 4.
A total of 12 structural parameters were selected as the design variables, namely, the elastic modulus, bulk modulus, shear modulus, density, cohesion force, internal friction angle of the embankment fill, and soil foundation. Considering the uncertainties of the material sources and construction qualities, the values of the above parameters were highly random [29]. According to the geological survey data of the Xi'an-Baoji expressway, the laboratory test results of the embankment fill and soil foundation, and other existing researches [18,47], the elastic modulus, bulk modulus, shear modulus, cohesion force of the embankment fill and soil foundation followed the log-normal distribution; the density, internal friction angle of the embankment fill, and soil foundation followed the normal distribution; the probability distribution parameters of these structural parameters are shown in Table 1, and the probability distribution curves of these structural parameters are shown in Figure 5.

Selection of the Embankment Seismic Damage Parameter.
According to Song et al. [48] and Li et al. [49], PGA is selected as the intensity parameter of the ground motions, and the embankment seismic damage levels are divided into 5 grades: basically intact, minor damage, moderate damage, severe damage, and destruction. Since the embankment is a Shock and Vibration 3 statically indeterminate structure, damages caused by earthquakes would be reflected by the displacement. erefore, ε max is selected as the embankment seismic damage parameter based on the displacement failure criterion, and the calculation method of ε max is shown in the following equation: where d max is the lateral maximum permanent displacement on the top surface of the embankment and D is the embankment width; D � 24.5 m. According to the investigation results of the Wenchuan earthquake, embankments at the epicenter (Yingxiu town) suffered from the most severely damage, i.e., destruction, and ε max reached 2.059% [50][51][52]; therefore, considering ε max � 2.0% as the boundary of severe damage and destruction was reasonable. Besides, the boundaries of ε max among other embankment seismic damage grades were further established based on the equidistant classifying method [17,53], as summarized in Table 2.

Principal Component Analysis Results.
Since the solution of the response surface function easily falls into a "curse of dimensionality" when there are too many structural parameters and the parameters may not be linearly independent, the principal component analysis was used to reduce the number of variables [36,54]. According to the probability distribution parameters listed in Table 1    calculations. e calculation method is shown in the following equation: where X i is the ith structural parameter, i is a natural number from 1 to 12, X imax and X imin are the maximal and minimal values of X i , and u X i is the average value of X i , as shown in Table 1.
and v(X j ) is calculated according to the following equation: where , and σ v(X j ) are the mean values and standard deviations of v(X i ) and v(X j ), respectively, and   Cohesion force of the embankment fill Log-normal distribution 34.00 kPa 0.20 X 10 Cohesion force of the soil foundation Log-normal distribution 31.00 kPa 0.20 X 11 Internal friction angle of the embankment fill Normal distribution 33.00°0.20 X 12 Internal friction angle of the soil foundation Normal distribution 28.00°0.20  Shock and Vibration other symbols have the same meaning. e calculation results are shown in Table 3. e weight vector, eigenvalue, original information contribution rate, and cumulative original information contribution rate of each structural parameter corresponding to the principal component were further calculated based on SPSS22.0, and 3 principal components F 1 , F 2 , and F 3 were extracted according to the eigenvalue not less than 1, as shown in Table 4. e equations of F 1 , F 2 , and F 3 are shown in equations (4)-(6) according to Table 4.

Shock and Vibration
e following conclusions can be obtained from equations (4) to (6): (1) F1 mainly reflects the elastic modulus, density, and cohesion force of the soil foundation, and the cohesion force of the embankment fill; F2 mainly reflects the elastic modulus and bulk modulus of the embankment fill, and the bulk modulus of the soil foundation; F3 mainly reflects the shear modulus and internal friction angle of the embankment fill, and the shear modulus of the soil foundation. (2) When X 1 to X 12 all take the mean values, F 1 , F 2 , and F 3 can also take the mean value 0. When X 4 , X 6, X 7 , and X 12 take the minimal values and other structural parameters take the maximal values, F 1 takes the maximal value, or otherwise, the minimal value. When X 1 , X 8 , and X 10 take the minimal values and other structural parameters take the maximal values, F 2 takes the maximal value, or otherwise, the minimal value. When X 2 , X 3 , X 7 , X 8 , and X 11 take the minimal values and other structural parameters take the maximal values, F 3 takes the maximal value, or otherwise, the minimal value.

Calculation Results of the Response Surface Method.
Mapping relationships between ε max and F 1 , F 2 , F 3 , PGA, and standard deviation σ 2 of ε max and F 1 , F 2 , F 3 , PGA were  Shock and Vibration 7 calculated based on the uniform design response surface method [17], as shown in the following equations: where α 0 -α 14 , β 0 -β 14 are undetermined coefficients and ε 1 and ε 2 are regression coefficients of ε max and σ 2 , respectively. e main object of this paper was to reveal the dynamic response regulars of the embankment under the near-fault pulse-like ground motions with different rupture fault distances.
In dynamic response analysis, more accurate results are often obtained when more than one earthquake are used, but according to the stochastic modeling and artificial synthesizing method of the near-fault pulse-like ground motions, it is necessary to calculate the attenuation relationship of PGV with the rupture fault distance, which requires a large number of seismic records [55]. In this paper, we calculated the attenuation relationship of Chi-Chi earthquake with M w � 7.6, and the 15 synthesized ground motions generated from the Chi-Chi earthquake were employed to the dynamic response analysis and the results were valid for earthquakes with the same characteristics [56]. Dynamic response analysis based on more earthquakes would be conducted in future researches. Since the PGA of the Chi-Chi earthquake was 1.193 g, three values of 0.2 g, 0.7 g, and 1.2 g were assumed to cover the recorded PGA of the Chi-Chi Table 3: Correlation coefficient matrix ρ ij of v(X i ) and v(X j ).  earthquake. PGAs of the 15 synthesized ground motions were adjusted to 0.2 g, 0.7 g, and 1.2 g, respectively, and 45 ground motions were obtained. 26 combinations were acquired through the experimental design of 3 levels of F 1 , F 2 , F 3 and PGA, the above 45 ground motions were used to conduct 390 dynamic response analysis for each combination. In addition, a 27th combination was added, and the 15 original ground motions were used for dynamic response analysis, and a total of 405 dynamic response analysis were conducted. e analysis results are shown in Table 5.
According to the dynamic response analysis results of the first 26 combinations, equations (7) and (8) were solved, and the regression coefficients are shown in Table 6.

Embankment Seismic Damage Analysis
4.1.1. Impact of PGA. Based on the dynamic response analysis results, the embankment seismic damage grades are basically intact and minor damage when PGA � 0.2 g, but due to the impacts of the structural parameters and rupture fault distances, ε max has a high dispersion, ranging from 0.1219 to 0.4532. With the increase of PGA, the embankment seismic damage is intensified, most of the damage grades are minor damage and moderate damage when PGA � 0.7 g, and ε max ranges from 0.5791 to 1.0816; the damage grades are all severe damage and destruction when PGA � 1.2 g, and ε max ranges from 1.5007 to 2.3301. Figure 6 shows the embankment seismic damages of R � 10 km for the 11th, 14th, 18th, and 27th combinations. It can be seen that when the structural parameters and rupture fault distances are constant, the damage grade and PGA have a positive correlation. Embankment seismic damages are mainly manifested by the lateral displacement on the top surface of the embankment and cracks (red solid lines in the figures). Among them, all the cracks originate from the deformation of the slope and extend to the interior of the embankment. e lateral displacements and length of the cracks are also positively correlated [57][58][59]. Figure 7 shows the embankment seismic damages of the 11th combination with R � 1 km, 4 km, 7 km, 10 km, 13 km, and 15 km. Combined with the analysis results of other combinations, it is noticeable that the correlation between the rupture fault distance and ε max is not obvious when the structural parameters and PGA are constant. Figure 8 shows the embankment seismic damages of the 27th combination with R � 1 km, 4 km, 7 km, 10 km, 13 km, and 15 km. It can be seen that different rupture fault distances correspond to different seismic damages to the same structure in the same earthquake, i.e., the smaller the rupture fault distances, the more severe the damages; and the larger the rupture fault distances, the milder the damages.    Figure 9 shows the embankment energy dissipations of the 27th combination with R � 1 km to 15 km, respectively. It can be seen that the energy dissipations are between 3.048 × 10 6 and 4.626 × 10 6 N·m; the smaller the rupture fault distances, the greater the energy dissipations; and the larger the rupture fault distances, the smaller the energy dissipations [59].

Embankment Seismic Fragility
Curves. By inputting the 10,000 sets of structural parameter combinations obtained in Section 3.3 to equation (7), obtain the probabilistic seismic demand model of the embankment [60][61][62], as shown in the following equation: ln ε max � ln 1.5897 + 1.3457 ln PGA.
Substituting equation (7) into the classical calculation formula of the seismic fragility, the following equation is obtained: where R f is the exceeding probability of the embankment seismic damage grade f, S f is the structural performance level shown in Table 2, that is, S f � 2.00 when f is destruction, S f � 1.20 when f is severe damage, S f � 0.60 when f is moderate damage, and S f � 0.30 when f is minor damage. e embankment seismic fragility curves plotted from equation (10) are shown in Figure 10.
Exceeding probabilities of each embankment seismic damage grade under different PGAs were calculated according to Figure 10 [63,64], as shown in Table 7.
It can be seen from Table 7 that the exceeding probability of the severe damage is 0.713094 and of the destruction is 0.323044 when PGA � 1.0 g, while according to Yin et al. [47], the two exceeding probabilities were 0.634172 and 0.265767, respectively. According to Yu et al. [65], Zhang et al. [66], and Luo et al. [67], the exceeding probability of each seismic damage grade of embankments at the epicenter of the Chi-Chi earthquake (PGA � 0.977 g) is shown in Table 8.
By comparing Tables 7 and 8, it can be seen that the research results of this paper are consistent with the actual embankment seismic damage conditions of the Chi-Chi earthquake, which indicates that the method proposed in this paper is scientific and reasonable. It also shows that it would obviously overestimate the seismic performance in  the embankment seismic fragility assessment without considering the uncertainties of the ground motions and structural parameters.

Conclusion
In this paper, stochastic modeling and artificial synthesizing method of the near-fault pulse-like ground motions was proposed, and a total of 12 structural parameters were selected as the design variables by taking the Xi'an-Baoji expressway K1125 + 470 embankment as the research object. In order to eliminate the linear correlations of these parameters, 3 principal components with large impacts on the embankment seismic fragility were extracted based on the principal component analysis. Mapping relationships among the principal components and embankment seismic damages were analyzed using the uniform design response surface method, and the seismic fragility assessment was carried out and the fragility curves were plotted. e following conclusions can be drawn: (1) Since the orientation of the apparatus in a seismic station is arbitrary, the recorded ground motions at the time of the earthquake are not necessarily the maximum, especially for the near-fault pulselike ground motions, the seismic intensities in different directions are obviously different. e stochastic pulse models of different rupture fault distances were established by considering the directional impacts of the pulse-like ground motions, and 15 ground motions with the rupture fault distance ranging from 1 to 15 km were synthesized by taking the Chi-Chi earthquake in Taiwan, China, as an example.
(2) Embankment seismic damages are mainly manifested by the lateral displacement on the top surface of the embankment and cracks. All the cracks originate from the deformation of the slope and extend to the interior of the embankment, and the lateral displacements and length of the cracks are positively correlated. It can also be seen that different rupture fault distances correspond to different seismic damages and energy dissipations to the same structure in the same earthquake, i.e., the smaller the rupture fault distances, the more severe the damages and the greater the energy dissipations; and the larger the rupture fault distances, the milder the damages and the smaller the energy dissipations. (3) e research results of the embankment seismic fragility assessment are consistent with the actual embankment seismic damage conditions of the Chi-Chi earthquake, indicating that the method proposed in this paper is scientific and reasonable. It also shows that it would obviously overestimate the seismic performance in the embankment seismic fragility assessment without considering the uncertainties of the ground motions and structural parameters. (4) Supporting structures such as retaining wall and antislide pile have positive effects on improving the embankment seismic performance. Defining the characteristics of the damages of the supporting structure-embankment fill-soil foundation system under different ground motions is a prerequisite to uncover the dynamic coupling mechanism between the supporting structure and the embankment, but relevant researches have not been carried out yet.

Data Availability
All the data generated or analyzed during this study are included within this article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the work submitted.