Fragility curves accounting for uncertainties in material parameters and ground motion characteristics using a data driven surrogate model

Material parameters

The material model adopted here for the RC members and diagonal masonry struts is defined by seven parameters for concrete and masonry: the peak compressive strengths (f_c′, f_m′); the strains at the peak strength (ε_1c, ε_1m); the residual strengths (f_rc, f_rm); the strains at the onset of the residual strength (ε_2c, ε_2m); the parameters lambda (λ_c, λ_m) defining the ratio of the unloading stiffness at the onset of the residual strength to the initial stiffness; the tensile strengths (f_tc, f_tm); and the tension softening stiffnesses (E_tsc, E_tsm). The damping ratio, ξ, of the structure is numerically modeled using the mathematically convenient Rayleigh formulation, considering 3% of the critical damping for the first two modes. To incorporate the effect of the uncertainty associated with damping, it is also considered as one of the variables.

Sensitivity analyses are conducted to quantify the effect of each of these 15 parameters on the structural response. In the sensitivity study, one parameter is varied at a time; once to a lower and once to a higher value. The variation range of each variable is selected to reflect the level of confidence in the selected value. The strength, deformation, and stiffness parameters are perturbed by 20%, 30%, and 40%, respectively, from the calibrated values of the numerical model, as those values are selected based on information available from material tests in Nepal^[31,32]. However, no information is available on the calibration of the parameter lambda (λ_c, λ_m) and the damping ratio; therefore, these are varied by 80% of the calibrated values to investigate their effects on the structural response. The calibrated properties, along with their percentage variations in the sensitivity analysis, are presented in Table 1. A total of 30 nonlinear models are created and, in each case, a nonlinear model is subjected to the two components of the ground motion recorded at the Univ Grants Comm., Sanothimi, Bhaktapur (THM) station, the closest station to the school building^[33,34]. The modeling uncertainties associated with the parameters used to represent reinforcing steel are neglected in this study, as they are not expected to vary considerably.

The results of the sensitivity analysis are summarized in Figure 2. The figure demonstrates the variability in the peak first-story drift values in the two orthogonal directions for the response history analyses compared to the values obtained from the baseline model. The first story drift is selected for this comparison, as it is the most critical engineering demand parameter for these structures^[10]. The results indicate that the response is sensitive to eleven parameters, which are shaded in Figure 2. These parameters are considered as input variables in the surrogate model developed in this study.

Figure 2. Variations in the peak ISD model in two orthogonal directions as compared to the values obtained from the calibrated model during the sensitivity analysis.

Ground motion selection and parameters

To incorporate the ground motion uncertainties, the bi-directional ground motion time-histories from the mainshock and the three major aftershocks of the 2015 Gorkha earthquake^[28,29] recorded at five stations in the Kathmandu valley [KATNP (Kanti Path, Kathmandu, Nepal), PTN (Pulchowk Campus, Tribhuvan Univ., Patan), THM, KTP (Municipality Office, Kirtipur), and TVU (Dept. Geology, Tribhuvan Univ., Kirtipur)] are considered to obtain a set of 20 ground motions. Due to the lack of other recordings from Nepal, FEMA P695^[35] is used to obtain additional motions in this study. FEMA P695 provides 22 far-field pairs of ground motions representing magnitudes in the range of 6.5-7.6 recorded on firm soil. Five of these 22 pairs of recorded motions and five from the set of 20 ground motions recorded in Nepal, summarized in Table 2, are employed here to run dynamic analyses to obtain data for the calibration of the surrogate model. These ground motions are selected to represent different soil types, frequency contents, amplitudes, etc.

The individual ground motions are first normalized by their peak ground velocities to remove unwarranted variability between records due to inherent differences in event magnitude, distance to the source, source type, and site conditions, while still maintaining the inherent record-to-record variability. To achieve statistical stability, the records are scaled to minimize the difference with the response spectrum used in Nepal (NBC 1994). The selected scale factors also ensure that the median of the geometric mean response spectra of all ground motions matches the response spectrum used in Nepal. Figure 3 shows the comparison between the NBC 1994 design spectrum and the median spectrum of the ten selected ground motions, scaled to the design hazard level considered in Nepal.

Figure 3. Acceleration response spectra of the selected ground motions.

The selected records are used here to conduct incremental dynamic analysis (IDA)^[36] to correlate a performance indicator (PI), such as the peak inter-story drift (ISD), with an intensity measure (IM) of the seismic excitation. The peak first story drift ratios in the two orthogonal dimensions are selected as PIs, as discussed in the sensitivity study. A number of IMs have been used in the literature^[37-39] besides the commonly used peak ground acceleration (PGA). Each of these IMs considers different characteristics of the ground motion with varying levels of influence on the structural response. Anastasopoulos et al. in their study on the response of bridge piers concluded that it is not possible to estimate the structural response satisfactorily using a single IM^[40]. Hence, a series of statistically significant IMs are considered in this study, including the following. In all cases, the IMs are considered for the ground motions in both the orthogonal directions.

• the PGA, peak ground velocity (PGV), and peak ground displacement (PGD);

• the Arias Intensity, I_A, proportional to the integral of the squared acceleration A(t) time history defined according to Equation 1;

• the spectral IMs in terms of acceleration (S_a), velocity (S_v), and displacement (S_d);

• the pseudo-spectral IMs in terms of acceleration (PS_a), velocity (PS_v), and displacement (PS_d);

• the Housner Intensity, I_H, integral of the pseudo-velocity response spectrum over the period range of 0.1 to 2.5 s^[41] according to Equation 2;

• the root-mean-square measures in terms of acceleration (A_RMS), velocity (V_RMS), and displacement (D_RMS) as the square roots of the mean values over the duration of the ground motion record;

• the characteristic Intensity, I_C, defined according to Equation 3;

• the sustained maximum measures in terms of acceleration (SMA) and velocity (SMV) defined as the third highest absolute peaks in the time histories^[42];

• the spectrum IMs in terms of acceleration (ASI) and velocity (VSI), defined according to Equation 4 as the integral of the spectral acceleration (S_a) and velocity (S_v), respectively, over the period range of 0.1 to 0.5 s^[43];

• the acceleration parameter, A95, defined as the level of acceleration which contains up to 95% of the Arias intensity^[44];

• the strong motion duration or significant duration, D_sig, defined as the interval of time between the accumulations of 5% and 95% of Arias intensity^[45];

• the Predominant Period, T_P, estimated using the 5% damped acceleration response spectrum at its maximum, as long as T_P > 0.2 s^[46];

• the mean Period, T_mean, obtained from the Fourier amplitude spectrum, C_i for each frequency f_i within the range of 0.25 to 20 Hz, according to Equation 5^[47];

• the effective peak acceleration (EPA), defined as the mean of the spectral acceleration over the period range of 0.1 to 2.5 s according to Equation 6^[48];

• the spectral acceleration measure, S*, defined in terms of pseudo-spectral accelerations PS_a(T₁) and PS_a(T₂), at the first and second mode periods, respectively, according to Equation 7^[49];

• the cumulative Absolute Velocity, CAV, defined in terms of the minimum and maximum accelerations according to Equation 8^[50].

where N is the number of 1-second windows in the time history, PGA_i and t_i are the maximum recorded acceleration (g) and start time of the window i, respectively, A_min is an acceleration threshold (user-defined, taken as 0.025 g typically) to exclude low-amplitude motions, and H(x) is the Heaviside step function;

• the standardized CAV (SCAV), defined in terms of the CAV according to Equation 9^[50];

• the normalized energy density (E_D), defined in terms of the velocity time history according to Equation 10^[43];

• the IM, I_zdefined in terms of the acceleration history and PGA, PGV (Cosenza and Manfredi, 1998) according to Equation 11^[46];

• the uniform duration, D_U, the summation of the time windows where the ground acceleration is above a certain threshold, which is considered 0.05 g in this study^[44];

• the number of peaks, defined as the number of times the acceleration time history crosses threshold acceleration values -0.025 g (N_p025), 0.05 g (N_p050), and 0.1 g (N_p100).

Statistical modeling

Nonlinear regression equations are developed herein to predict the seismic response of the four-story school building with functions of the material properties and ground motion characteristics. To this end, the 3SLS statistical modeling approach used by Anastasopoulos et al. to examine the effectiveness of various IMs in predicting the structural response of a bridge pier is adopted^[40]. This method correlates the selected damage indices with a number of statistically important IMs and material properties. It also takes into account the errors associated with bias due to omissions of variables and the regression properties, such as heteroscedasticity, autocorrelation, and exogeneity of the regressors, among others^[51].

Assessment of the surrogate model

To assess its accuracy, the statistical model is used to estimate the PIs for the models and ground motions used for its training. Figures 4 and 5 compare the interstory drift values obtained from the detailed nonlinear time history analyses and those estimated using the surrogate models. The comparisons indicate reasonably good matches between the statistical and the FE models. The fits are better for the higher drift levels, which are of interest when structural damage and/or collapse are estimated. This is particularly encouraging as the primary focus in this study is the estimation of the collapse potential of the structure, which occurs at larger drifts.

Figure 4. Observed and predicted values in terms of the peak ISD for the case-study building. (A) Peak ISD along the X direction; (B) Peak ISD along the Y direction.

Figure 5. Observed and predicted values in terms of the residual ISD for the case-study building. (A) Residual ISD along the X direction; (B) Residual ISD along the Y direction.

To further assess the confidence in the surrogate model, it is used to predict the dynamic responses of numerical models with all 11 MP values different than those presented in Table 3 used in the training process. These models are also subjected to ground motions not considered in the development of the surrogate model. These “unseen” ground motions are taken randomly from the set of 42 ground motions, which comprises the far-field suite of FEMA P695 and motions from the 2015 Gorkha earthquake, with the exclusion of those used to train the surrogate models. The performance of the surrogate model is evaluated for ten such cases using five “unseen” ground motions and two sets of “unseen” MPs. The two sets of MPs obtained by estimating the average of Values 1 and 2 and the average of Values 3 and 4 in Table 3, respectively for all MPs. For example, the values used for f_c′ are 1.19 ksi and 1.61 ksi. The five “unseen” ground motions used here are presented in Table 6. The results of these analyses, marked as blue stars in Figures 4 and 5, are rather encouraging, considering the complex nonlinear behavior of the structure.

Table 6

“Unseen” ground motions selected to evaluate the surrogate model

Earthquake			Recording station	PGA (g)
Mw	Year	Name
7.3	2015	Gorkha	Municipality Office, Kirtipur (KTP)	0.069
6.7			Univ Grants Comm., Sanothimi, Bhaktapur (THM)	0.097
6.7	1994	Northridge	Canyon County - WLC	0.482
6.9	1995	Kobe	Shin-Osaka	0.243
6.9	1989	Loma Prieta	Gilroy Array #3	0.555

Uncertainties in model parameters

The statistical models summarized in Table 4 indicate that only six of the eleven sensitive parameters identified in Figure 2 are statistically significant: the peak compressive strengths and corresponding strains of concrete and masonry, the tensile strength of masonry, and the damping coefficient. Therefore, these are considered in the MC simulations. The probability density functions of these parameters are estimated based on data from material tests on concrete and masonry available in the literature^[54-58]. In previous studies, the distributions are typically assumed to be lognormal^[20,59]; however, to remove any possible bias due to the selection of the probability density functions and investigate which distribution is more suitable for these MPs, both Weibull and lognormal distributions are fitted here, as shown in Figure 6.

Figure 6. Fitting cumulative and probability distributions to available test data. (A) f_c′ from the lab tests in Taiwan^[54]; (B) f_c′ from the actual buildings in the USA^[56,57]; (C) f_c′ from the lab tests in UB, USA^[58]; (D) ε_1cfrom the actual buildings in the USA^[56,57]; (E) ε_1cfrom the lab tests in UB, USA^[58]; (F) f_m′ from an actual building in the USA^[56]; (G) f_m′ from the lab tests in UB, USA^[58]; (H) ε_1mfrom an actual building in the USA^[56]; (I) ε_1mfrom the lab tests in UB, USA^[58]; (J) f_m′ from the lab tests in IITK, India^[55]; (K) ε_1mfrom the lab tests in IITK, India^[55]; (L) f_tm from the lab tests in IITK, India^[55].

It can be observed in Figure 6 that both distributions fit the test data equally well. Hence, a lognormal distribution is assumed in line with previous studies. Table 7 presents the median, e^μ, and standard deviation, σ for a lognormal distribution adopted for all the model parameters. In the case of the damping ratio, ξ, for which test data is very limited, the parameters are taken from the literature^{[20,54,60-62]}. The standard deviation values of the datasets in Table 7 that have mean values closest to the parameters of the calibrated model summarized in Table 1 are employed here to generate samples for the MC simulations. These are indicated with bold text in the table.

Results

MC simulations are conducted here using the statistical models to investigate the effects of the uncertainties associated with the random variables (IMs and MPs) on the fragility curves. To generate the models for these analyses, the MPs are sampled N_samp times, generating a set of realizations consistent with the assumed lognormal distribution. Then, two approaches are considered to generate the cases analyzed with the statistical models.

Approach-1: The realizations are set up such that the first random variable is varied from 1 to N_samp while keeping the other random variables constant. Then, the second variable is varied from 1 to N_samp without changing the other variables, and so on. Thus, this approach can consider the interactions between the random variables. However, the number of models (N_mod) for each simulation increases exponentially as $$N_{mod}=\left(N_{samp}\right)^{N_{MP}}$$, where N_MP is the number of input MPs varied.

Approach-2: With this approach, N_samp sets of random variables are generated, with each set representing one realization. In this case all MP values are different between the different models. Thus, solving the problem N_samp times gives information on the randomness in the output response of the system to each set of input variables. Thus, the number of models (N_mod) in this approach is equal to N_samp.

All the material properties considered here for developing the set of structural models are assumed to be uncorrelated. For each set of realizations of MPs, the statistical surrogate model is used to calculate the peak ISDs for the ten selected ground motions. The IMs are increased incrementally to perform the additional IDAs incorporating both modeling and record-to-record uncertainties. The geometric mean of the spectral acceleration at the first modal period of the building, S_a[T₁] of the two components of each ground motion is selected as the IM, while the maximum drift ratio is selected as the PI to develop the IDA curves. The periods of the structure can differ from one model to another as some MPs can affect the stiffness; however, for consistency, the period of the calibrated model is used in all cases.

To examine the effects of modeling uncertainties at different levels of the structural behavior, distinct limit states in the form of damage grades (DGs) are considered for the case study building. For the damage classification, a critical aspect of post-earthquake damage assessment^[63-65], the damage classification used by Brzev et al. to assess low-rise RC buildings in Nepal affected by the 2015 Gorkha earthquake is adopted here^[65]. This damage classification scheme includes five DGs ranging from DG1 to DG5 based on the EMS-98 scale, to characterize the severity of damage in structural components. This damage classification is correlated to maximum drifts reached during an earthquake by using data from laboratory specimens^[59,66] and by consulting with the field engineers from the National Society of Earthquake Technology (NSET), Nepal^[67] who implemented it after the 2015 earthquake. As proposed by Hilly^[67], the limit states, DG1 to DG5, correspond interstory drift ratios of 0.3%, 0.6%, 1.2%, 2.0%, and 3.0%, with DG5 representing severe damage and partial or total collapse of the building. The medians and standard deviations of S_a[T₁] for the five DGs obtained from both approaches of MC simulations are summarized in Figures 7 and 8. It can be observed that the medians and standard deviations converge to the same values of S_a[T₁], representing the five DGs for both approaches. However, the first approach needs 46,656 models for all median and standard deviation values to converge compared to the 100 analyses needed for the second approach. Therefore, the second approach is adopted here for further analyses.

Figure 7. Medians and standard deviations of the S_a[T₁] at various damage grades obtained from MCS - Approach 1. (A) Medians; (B) Standard deviations.

Figure 8. Medians and standard deviations of the S_a[T₁] at various damage grades obtained from MCS - Approach 2. (A) Medians; (B) Standard deviations.

To capture the uncertainties associated with record-to-record randomness but not the MPs, IDA is also performed for the baseline model and the ten ground motions considered here. The IDA curves of the deterministic model and the extended IDA plot for N_mod = 100 are presented in Figure 9.

Figure 9. IDA curves for the calibrated model and from Monte Carlo Simulation using the surrogate model. (A) Baseline model; (B) Monte Carlo simulations Approach 2, N_mod = 100.

Table 8 presents the medians and standard deviations obtained from the MC simulations, the baseline FE model, and the baseline surrogate models. For all DGs, the methods considered here result in similar trends for the median and dispersion values, thus further validating the use of the statistical models for MC simulations. However, one can observe that for all DGs, the median values of the spectral acceleration are lower in the case of the MC simulations, which account for the variability in the MPs. Moreover, one can note that the standard deviation increases by approximately 10% in the case of the surrogate model compared to the FE model. The increase is larger in the case of the MC simulations, and especially for the larger DGs. The substantially higher scatter can be expected, considering that in the case of the MC simulations, the uncertainties in the MPs are also included. This is evident in Figure 9, which demonstrates a higher scatter in the data points.