Metamodel-Based Seismic Fragility Analysis of Concrete Gravity Dams

Methods for the seismic analysis of dams have improved extensively in the last few decades, and the growth of computing power has expedited this improvement. Advanced numerical models have become more feasible and, thus, constitute the basis of improved procedures for design and assessment. Moreover, a probabilistic framework is required to manage the various sources of uncertainty that may impact the system performance and decisions related thereto (Ellingwood and Tekie 2001). Fragility analysis, which depicts the conditional probability that a system reaches a structural limit state, is a central tool in this probabilistic framework. Traditional vulnerability assessment methods develop fragility functions by using a single parameter to relate the level of shaking to the expected damage, which consequently produces a robustness of predictions that is highly dependent on the selected parameter. However, the estimation of the fragility of the system can be potentially improved by increasing the number of parameters; in this way, a more complete description of the properties of ground motions can be obtained (Alembagheri 2018). Furthermore, the effect of the variation of the material properties in the seismic fragility analysis of structures with complex numerical models, such as dams, is frequently overlooked due to the costly and time-consuming revaluation of the numerical model. The seismic response and vulnerability assessment of key infrastructure elements often require a large number of nonlinear dynamic analyses of complex finite-element models (FEMs). The substantial computational time may be reduced by using machine learning techniques to develop a surrogate or metamodel, which is an engineering method used when an outcome of interest cannot be easily directly measured; thus, a model of the outcome is used instead (Forrester et al. 2008). In addition, if the outcome of interest comes from nonlinear FEM analysis that reflects the dynamic behavior of the structure under seismic loading, the metamodel will emulate this behavior. To this end, these algorithms include several features in their mathematical formulation to help capture this highly nonlinear behavior. Among them, the use of higher order sparse polynomials, partitioning the sample space and fitting a series of models and then combines them into an ensemble with an overall better performance and by mapping into high dimensional feature space, between others.

Such a challenge is particularly relevant to the case of large-scale infrastructures such as dams subjected to seismic loads. Thus, the main goal of this study is to explore the applicability of metamodels for the seismic assessment of gravity dams and present a methodology to develop parameterized multivariate fragility functions through the use of the latter. The secondary goal is to explicitly account for the effect of the model parameter variation in the seismic fragility analysis. The proposed methodology has the added asset of properly depicting the seismic scenario likely to occur at a specific site, enhancing the accuracy of the seismic fragility analysis. The proposed methodology is applied to a case study gravity dam located in northeastern Canada.

In the past two decades, univariate fragility functions, which depict the potential for limit state exceedance conditional on a ground motion intensity, have been readily adopted and developed for the seismic assessment of structures, as it can been found in several state-of-the-art reviews (Ghosh et al. 2017; Hariri Ardebili and Saouma 2016b; Muntasir Billah and Shahria 2015). However, given the recognized limitations of such univariate fragility curves, the use of multivariate fragility functions to assess the vulnerability of a given structure or portfolio of structures is beginning to be used progressively (Brandenberg et al. 2011; Koutsourelakis 2010; Pan et al. 2010; Gehl et al. 2009; Grigoriu and Mostafa 2002). Noted advantages of these parameterized or multivariate fragility functions include the potential for efficient posterior uncertainty propagation, exploring sensitivities or the influence of design parameter variation, and enabling application across a portfolio of structures. Nevertheless, given the large number of simulations required, the development of fragility surfaces or multivariate fragility functions that leverage numerical models impose high computational burdens.

The combination of numerical models, probabilistic approaches, and machine learning has gained considerable interest in the literature in recent years for engineering design and structural reliability (Seo and Rogers 2017; Yu et al. 2014; Sudret 2012; Wang and Shan 2007; Simpson et al. 2001). This combination is justified by the significant randomness that characterizes not only the earthquake excitation but also the structural system itself (e.g., stochastic variations in the material properties, degradation due to aging, and temperature fluctuation, etc.). Surrogate modeling techniques within a seismic fragility framework have found recent applications for the safety assessment of buildings and bridges, among other structures (Mangalathu and Jeon 2018; Sichani et al. 2017; Kameshwar and Padgett 2014; Ghosh et al. 2013; Seo and Linzell 2013; Seo et al. 2012). Even though many of these studies considered several seismic intensity measures (IMs) and model parameters (MPs) for building the metamodels to predict the response of the structure, most of them do not clearly depict the influence of all the considered parameters in the form of multivariate fragility functions from the respective metamodels. For the specific case of dam-type structures, an extensive comparison between machine-learning data-based predictive models for monitoring the dam behavior can be found in Salazar et al. (2015a, b). More recently, Hariri-Ardebili and Pourkamali-Anaraki (2017a, b), and Hariri-Ardebili (2018) have used machine learning techniques to perform reliability analysis applied to gravity dams against flooding, earthquakes, and aging, considering, in some cases, explicit limit state functions and simplified FEM in others.

Most of the prior studies on the seismic assessment of concrete gravity dams via machine learning techniques are limited to the consideration of a single metamodel, simplified FEMs, and univariate fragility functions. Therefore, they do not explore the most suitable metamodel for fragility analysis of this type of structure, nor they explore the influence of the variation of the model parameters on the seismic fragility analysis. Moreover, none these studies discuss the proper definition of the seismic scenario likely to occur at a specific site in a probabilistic manner.

To address the aforementioned gaps, this paper aims to identify the most viable metamodel, from the subset of machine learning methods considered, for the seismic fragility assessment of gravity dams, provide an overview on the importance of the parameters influencing the dam performance, and formulate design recommendations from the analysis. The major contributions of this paper, in order of presentation, can be listed as follows: (1) perform a comparative analysis to determine the best performing regression metamodel to predict the base sliding of gravity dams from six metamodels with different basis function configurations, yielding a total of 14 different regression techniques, for the first time; (2) present a methodology to fit parameterized fragility surfaces, as a function of IMs and MPs, from the metamodels; (3) consider the correlation between the seismic IMs from a probabilistic seismic hazard analysis (PSHA) to generate the samples to characterize the seismic scenario where the metamodel will be evaluated; (4) gain insight into the influence of the model parameters affecting the dam performance and explicitly quantify their effect with the generation of multivariate fragility functions; and (5) formulate model parameter design recommendations from the analysis, e.g., appropriate range of parameters to achieve target risk.

Dam seismic assessment is a complex task due to the uniqueness of each of such structures and to the interaction between the different components of the system. Similarly, the seismic response of dam-type structures involves nonlinear dynamic analysis of complex FEMs, often requiring prohibitively high computation times. Machine learning describes a series of methods that allow learning from data, i.e., what relationships exist between the quantities of interest. When performing probabilistic studies related to structural reliability, it may be interesting to replace the FEM by a regression model built on a set of simulated responses, for the purpose of computational efficiency (Goulet 2018). Accordingly, the motivation of applying statistical algorithms to develop a seismic probabilistic demand model or metamodel is to expedite this safety assessment process. Such a model of a model is based on machine learning techniques to allow the algorithm to learn from the data and because observations are the output of a simulation, the observation model does not include any observation error. Within this context, a probabilistic seismic demand model expresses an approximate relationship between an uncertain seismic response, e.g., a dam’s maximum relative base sliding, and a set of parameters that influence the response. The basic idea is for the surrogate or metamodel to act as a curve fit to the available data so that the results may be predicted without requiring costly simulation. In the general case, the metamodel can be described as follows:where the surrogate model $g(·)$ statistically predicts the response of the structure, $y_i$, for a given set of intensity measures and model parameters, ${\mathbf{x}}_i$; and $v$ = error due to the lack of fit of the surrogate model.

The metamodels considered herein are within an adaptive scheme, i.e., the functions in the metamodels can change according to the input data to reduce the burden of manual selection of several parameters in the metamodel. Therefore, three adaptive metamodels and an interpolation scheme will be considered. In addition, the performance of two other statistical learning algorithms based on kernels and decision trees will also be addressed. Among the different regression techniques, this paper focuses on the following: (1) polynomial response surface models of order 2 to 4 with stepwise regression (PRS ${\mathcal{O}}^{2–4}$); (2) multivariate adaptive regressive splines (MARS) with linear and cubic splines; (3) adaptive basis function construction (ABFC); (4) radial basis function (RBF) interpolation with multiquadratic, thin plate splines and Gaussian basis functions; (5) support vector machines for regression (SVMR) with linear, quadratic, cubic, and radial basis function kernels; and (6) random forest for regression (RFR).

Within the extent of this study, the first comparative analysis of metamodels in the context of seismic assessment of dams is performed. Using the results from the finite-element simulations, this study develops metamodels for approximating the seismic response of gravity dam-type structures using the three steps outlined in Fig. 1. The subsequent sections will detail such steps.

To minimize the associated cost of running dynamic nonlinear FEM of dams under seismic excitation while analyzing an adequate number of loading conditions and structural system configurations, an appropriate experimental design method should be used. Structural and material properties likely to affect the seismic response of the structure should be considered, and their associated ranges should be based on experimental data or values found in the literature. The Latin hypercube sampling (LHS) experimental design method is adopted to generate $n_f$ sample points representing the different configurations of the dam under study. This sampling technique was selected because of its ability to divide the desired range of values for each parameter into $n$-equiprobable intervals and then select a sample once from each interval.

The following subsections provide a brief overview of the different metamodels tested in this study for the ability to offer viable metamodels of the seismic response of dams. Only the most relevant features of each technique will be presented given that exhaustive mathematical formulation can be found in the literature. Alternatively, relevant details concerning the model fitting and the algorithm settings are provided.

The goodness-of-fit estimates depict the discrepancy between the observed values from the FEM simulation and the estimated value with the meta-model in question. The root mean square error (RMSE) provides a measure of global error, quantifying the difference between the responses predicted by the metamodel and actual data and is computed as follows:where $y_i$ = response values in the dataset; ${\widehat{y}}_i$ = predicted values; and $n_f$ = total number of points in the dataset (FE simulations). Additionally, the coefficient of determination $R^2$ can be calculated as follows:where ${\sigma }^2$ = variance of the response in the dam response dataset. Similarly, the relative maximum absolute error (RMAE) measures the extent of the local fitting error and is the ratio of the maximum absolute difference between the metamodel and test data responses to the standard deviation of the actual response

In the present study, the predictive capability of the metamodels will be assessed through 5-fold crossvalidation (5-CV). The dataset is randomly divided into 5 sets, and the metamodel is trained using $k-1$ sets, with the remaining set used as test data. This procedure is repeated 5 times; thus, 5-CV provides an estimate of the predictive accuracy of the model for unknown data. The average $R^2$ value resulting from 5-CV will be used along with RMSE and RMAE to compare the different metamodels.

As stated by Ghosh et al. (2013), single-parameter demand models suffer from two potential drawbacks: (1) the inability to assess the impact of structural model parameter variation on structure performance during earthquakes without costly reanalysis for each new set of parameter combinations, and (2) the lack of flexibility to incorporate field instrumentation data from monitoring of existing structures to enable the updating of seismic fragility estimates. Consequently, the use of multivariate fragility functions enables the efficient uncertainty propagation of the random variables and allows for the exploration of the effects of design parameter variation on the vulnerability of the structure. Thus, the goal is to identify the role of the most influential ground motion IMs and MPs on the induced damage in the structure to build multivariate fragility functions from the metamodel output to provide a more complete and accurate view of the vulnerability of the structure.

Similar to fragility curves, multivariate fragility functions offer the conditional probability of exceeding different limit states given the occurrence of an earthquake of a certain intensity. The only difference is that the specific limit state is characterized with $n$ parameters $p_1,p_2\dots p_n$ instead of one parameter, as is the case with fragility curves. Hence the probability of limit state exceedance is conditioned on the resulting set of critical parameters. The fragility function corresponding to the limit state $l$ is defined as follows:where LS = limit state damage index; and ${\mathrm{LS}}_l$ = value corresponding to the $l$th limit state.

The proposed methodology in this study is applied to a case study gravity dam in Quebec, Canada. It possesses 19 unkeyed monoliths, a maximum crest height of 78 m, and a crest length of 300 m. The dam was chosen for its simple and almost symmetric geometry, its well documented dynamic behavior, and the availability of forced vibration test results used to calibrate the dynamic properties of the numerical model (Proulx and Paultre 1997).

The tallest monolith of the dam is selected as representative of the structure and was modeled with the explicit finite-element software LS-Dyna version R. 9.1.0, as shown in Fig. 4, following the recommendations of the United States Bureau of Reclamation (USBR) (Mills-Bria et al. 2013). Only one load case combination was considered, which included self-weight, hydrostatic thrust, uplift, hydrodynamic effects, and seismic load. The proposed model takes into account the different interactions among the structure, reservoir, and foundation. The reservoir is modeled with compressible fluid elements, whereas the concrete dam and the rock foundation are modeled with linear elastic materials to which a viscous damping is associated. Given that the model should remain stationary after the static loads are applied, two loading phases were considered: (1) a dynamic relaxation phase for static loads, and (2) a dynamic phase for the seismic loads, each with different boundary conditions. In the first loading phase, a symmetric boundary condition was applied where the normal displacements are zero. For the dynamic phase, nonreflective boundaries were included to prevent artificial amplification of the seismic waves because of the finite length of the foundation and reservoir. Concerning the contact surfaces between the different components of the system, sliding contact with zero friction was used to model the dam-reservoir interface. For the reservoir-foundation interface, tied contact was applied, except near the upstream face of the dam, where sliding contact with zero friction was used to maintain reservoir load during the sliding of the dam. Preliminary linear analyses identified the concrete-rock interface at the base of the dam and the concrete-concrete interface at the crest of the dam as high areas of tensile stresses and, therefore, where cracking and sliding could occur. Consequently, the model nonlinearity was constrained to these two areas only, using tiebreak contact elements with a tension-shear failure criterion. Further details of the modeling assumptions and the validation of the numerical model can be found in Bernier et al. (2016).

Table 1 presents the parameters that were considered random variables in the numerical analysis of the dam response and for which the uncertainty was formally included through their probability distribution functions (PDFs). All of the remaining input parameters were kept constant and represented by their best estimate values. For the case study dam, due to limited availability of material investigations, the probability distributions were defined using empirical data of similar dams. The uniform distribution was used for most parameters except for damping, for which a log-normal distribution was adopted as proposed by Ghanaat et al. (2012). By posing the resulting metamodels and fragilities as functions of MPs that have a significant impact on the behavior of the dam, future studies may integrate these models with emerging estimates of these parameters or updated PDFs.

In recent years, typical damage modes that could lead to the potential collapse in dams after a seismic event have been identified, and seismic damage levels can be established. Preliminary analyses have confirmed sliding as the critical failure mode for the case study dam (Bernier et al. 2016), and other failure modes would only occur after sliding had already been observed. As a result, base sliding at the concrete-to-rock interface damage limit states proposed in Tekie and Ellingwood (2003) is considered in this study and presented in Table 2. The incipient sliding limit state should only cause minor damage since well-dimensioned dams should be able to undergo slight deformations or displacements while remaining stable. The moderate damage limit state can be considered as the onset of nonlinear behavior where material cracking occurs, deformations may become permanent and the drainage system begins to be affected. Displacement greater than 50 mm could cause differential movements between the blocks and potentially lead to the loss of control of the reservoir and very extensive damage, while displacements greater than 150 mm represent a complete damage state and high probability of collapse of the structure.

A probabilistic seismic hazard analysis (PSHA) was performed at the dam site with the computer software OpenQuake version Engine 2.8 to characterize the possible earthquake scenarios at different intensity levels. The hazard levels were defined in terms of horizontal spectral acceleration at the fundamental period of the structure [$S{a}_H(T_1)=0.1\mathrm{g}:0.1:1.0\mathrm{g}$] to conveniently cover the range of spectral accelerations corresponding to return periods from 700 to 30,000 years.

To proceed with the selection of a representative set of ground-motion time series (GMTS), the generalized conditional intensity measure (GCIM) approach (Bradley 2010) was adopted. The purpose of using the GCIM approach is to include the most influencing seismic IMs with respect to the structural response. For the case of gravity dam-type structures, peak ground velocity (PGV) was found to be one of the best performing ground-motion structure-independent scalar IM to correlate damage (Hariri Ardebili and Saouma 2016a). Similarly, the vertical spectral acceleration $S{a}_V$ is also expected to be relevant in heavy structures of this sort. As a result, the set of considered IMs in the GCIM are $S{a}_H(T)$, $S{a}_V(T)$, and PGV, where $S{a}_H(T)$ and $S{a}_V(T)$ are computed at 20 vibration periods in the range of $T=\{0.2T_1-2T_1\}$, as proposed by Baker (2011), leading to a total of 41 IMs to be considered in addition to the conditioning IM, $S{a}_H(T_1)$. The GCIM distribution computed with the abovementioned IMs was then used to simulate and select 250 ground motions. The records were selected from the PEER NGA-West2 database (Ancheta et al. 2013) due to the limited availability of strong ground motion records in the PEER NGA-East database (Goulet et al. 2014). Further details on the PSHA and the record selection procedure can be found in Segura et al. (2018).

By selecting different configurations of the model parameters in Table 1, $n_f=250$ samples of the FEM were generated with LHS and paired with the selected ground motions. The maximum relative sliding at the base (${\delta }_{\max }$) was computed from nonlinear simulations, and 14 regression metamodels were fitted to the structural response. An initial prescreening of the covariates or predictors was made before selecting the algorithm, based on visual inspection of the scatter plots of the possible predictors with respect to the response of the structure and taking into account the parameters affecting the dynamic behavior of the structure. This initial set of predictors (input) was used to train all the considered metamodels to perform a comparative analysis. In addition to the model parameters listed in Table 1, several seismic intensity measures were considered in the starting set of variables, such as spectral acceleration at the fundamental period [$S{a}_H(T_1)$], spectral velocity at the fundamental period [$S{v}_H(T_1)$], peak displacement, velocity and acceleration (PGD, PGV, PGA), spectrum intensity (SI), earthquake angular frequency (${\omega }_{eqk}$), significant duration (${\mathrm{D}}_{595}$), Arias intensity ($I_a$), and peak ground acceleration and spectral acceleration at the fundamental period in the vertical direction [${\mathrm{PGA}}_V$, $S{a}_V(T_1)$]. Given that some algorithms already perform an internal selection of the predictors (with forward and backward iterations), the selected final set of predictors (output) in each metamodel is a subset of the initial set of predictors.

Regarding the generation of the samples where the metamodel will be evaluated to predict the dam’s response and derive point estimates of fragility, independence among all the MPs was considered to generate $5\times {10}^5$ samples with LHS. Nevertheless, it should be noted that for a specific site, the seismic IMs are correlated. For the case study dam and considering the 250 GMTS used to train the metamodels, a linear correlation is to be expected, as seen from Fig. 8. Fig. 9 shows the histogram of the logarithm of the seismic IMs that follow an approximate normal distribution. Based on this, the samples of these parameters are taken from a jointly log-normal distribution with their respective correlation coefficients. In addition, to consider a range of values corresponding to return periods from 700 to 30,000 years at the dam site, the possible values of PGV, $I_a$, and ${\mathrm{PGA}}_V$ were bounded, respectively, as $0.8\mathrm{cm}/\mathrm{s}\le \mathrm{PGV}\le 25\mathrm{cm}/\mathrm{s}$, $0.0\mathrm{m}/\mathrm{s}\le I_a\le 2.5\mathrm{m}/\mathrm{s}$, and $0.01\mathrm{g}\le {\mathrm{PGA}}_V\le 0.25\mathrm{g}$.

The fragility analysis was performed following the MSA methodology depicted in Fig. 2. The range of each of the parameters of the fragility surface was divided into 100 intervals. As a result, ${10}^4$ fragility point estimates were generated for each limit state. Moreover, to display how the variation of the model parameters alters the seismic fragility, fragility surfaces as a function of PGV and each of the model parameters considered in the metamodel (CRF, DR, and CRC) were generated. It is noteworthy that PGV was selected as the IM to be displayed in the fragility surface not only because of its practicality of database availability and value accessibility but also because, as previously mentioned, it was found to be one of the best performing ground-motion IM to correlate with the proposed damage state (Hariri Ardebili and Saouma 2016a).

Parametric fragility surfaces for each limit state and for each MP were generated, implementing the methodology explained in Fig. 3. The resulting fragility surfaces are depicted in Figs. 10–12, and the goodness-of-fit between the proposed parameterized fragility surfaces and the fragility point estimates is presented in Table 4.

The resulting fragility functions can also be used to assess the seismic performance of the dam by formulating recommendations with respect to the model parameters. To achieve a desired seismic performance, boundaries of model parameters for an adequate performance under extreme limit states can be formulated. The usable range of values is determined by ensuring that the probability of exceedance given that an extreme event is in line with the current guidelines for the minimum provisions for life safety. To establish the admissible values of the model parameters, a return period of 10,000 years, which corresponds to the maximum considered earthquake (MCE) (ASCE 2016; FEMA 2015; CDA 2007), was used. The peak ground velocity associated with the MCE (${\mathrm{PGV}}_{MCE}$) was determined from the PSHA and from the fragility surfaces corresponding to the complete damage limit state (LS3), and the fragility curves as a function of the model parameters were extracted for the ${\mathrm{PGV}}_{MCE}$ value.

According to what it is proposed by ASCE 7 (ASCE 2016), for a risk category III (dam-type structures) and total or partial structural collapse, the maximum probability of failure should be less than 6%. Consequently, and as shown in Fig. 16, the boundaries of model parameters were established to provide a probability of exceedance less than 6% for an MCE seismic scenario. As a result, a range of values of $50°\le \mathrm{CRF}\le 55°$, $0.35\le \mathrm{DR}\le 0.66$, and $0.87\mathrm{MPa}\le \mathrm{CRC}\le 2.0\mathrm{MPa}$ should be considered to ensure that the performance of the dam is in line with the minimum values for life safety. Nonetheless, it should be acknowledged that these parameter ranges are derived from the single parameter at a time evaluation, and the joint interaction of the model parameters should be further examined.

The main objective of this study was to identify the most viable metamodel among those considered for the seismic fragility assessment of gravity dams, present a methodology for the development of multivariate fragility functions displaying the effect of the model parameter variation on the dam seismic performance and formulate design recommendations from the analysis.

The methodology was applied to a case study dam located in northeastern Canada. PSHA was performed to characterize the seismic hazard at the dam site and to select 250 ground motions, consistent with the latter, using the GCIM approach. The sets of selected ground motion records were paired with the 250 samples of numerical models of the dam generated with LHS, representing different material and loading configurations of the system. The dataset used to train the metamodels was generated by performing nonlinear dynamic analysis of these samples, with an FEM that considered the fluid-structure-foundation interaction and by extracting the maximum relative sliding at the base of the dam. Six different types of metamodels with different configurations (basis and interpolation functions) were fitted to the seismic response of the case study dam to predict the sliding limit state at the base of the dam. The 4th-order PRS emerged as the best performing metamodel based on local and global goodness-of-fit estimates from 5-CV, and it was used to generate fragility surfaces as a function of PGV and each of the model parameters. It is observed that the variability of the concrete-rock cohesion model parameter has the most influence on the fragility analysis estimates, followed by the drain efficiency and to a lesser extent, the concrete-rock angle of friction. Finally, a seismic assessment was performed to determine the model parameter boundaries for adequate performance under extreme events, such as the MCE, to respect the minimum safety margins proposed by the current guidelines.

It should be noted that machine learning techniques are indispensable when assessing the vulnerability of structures with computationally expensive FEM such as gravity dams subjected to seismic loading. Similarly, the use of surrogate models allows for the exploration of the impact of different parameters in the fragility without the costly reevaluation of the FEM simulations. Regarding the fragility functions, as evident for the case presented herein and the goodness-of-fit values, the fragility estimates are well depicted by the methodology suggested in this study to fit parametric fragility surfaces.

It is expected that the results of this study can lead to more accurate planning and retrofitting policies to expedite the safety assessment of dams under seismic loads while supporting the decision-making process and to guide the preliminary design of future gravity dams. In future studies, additional insight into the correlation between the parameters defining the model configurations should be made, including other relevant limit states for gravity dams. Additionally, the model parameter variations in the fragility analysis should be further explored to provide parametric fragility functions, including the joint interaction of these parameters using classification meta-modeling techniques.

The authors acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC), the Fonds de recherche du Québec Nature et technologies (FRQNT), the Centre d’études interuniversitaire des structures sous charges extrêmes (CEISCE), and the Centre de recherche en génie parasismique et en dynamique des structures (CRGP). Computational resources for this project were provided by Compute Canada and Calcul Quebec. The authors also thank Carl Bernier and all the students in Padgett Research Group at Rice University for their collaboration and valuable remarks on various aspects of this study. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the sponsors.