Shedding Light on the Effect of Uncertainties in the Seismic Fragility Analysis of Existing Concrete Dams

Abstract

The seismic risk assessment of existing concrete gravity dams is of primary importance for our society because of the fundamental role of these infrastructures in the sustainability of a country. The seismic risk assessment of dams is a challenging task due to the lack of case histories, such as gravity dams’ seismic collapses, which hinders the definition of limit states, thus making the application of any conventional safety assessment approach difficult. Numerical models are then fundamental to predict the seismic behaviour of the complex dam-soil-reservoir interacting system, even though uncertainties strongly affect the results. These uncertainties, mainly related to mechanical parameters and variability of the seismic motion, are among the reasons that, so far, prevented the performance-based earthquake engineering approach from being applied to concrete dams. This paper discusses the main issues behind the application of the performance-based earthquake engineering to existing concrete dams, with particular emphasis on the fragility analysis. After a critical review of the most relevant studies on this topic, the analysis of an Italian concrete gravity dam is presented to show the effect of epistemic uncertainties on the calculation of seismic fragility curves. Finally, practical conclusions are derived to guide professionals to the reduction of epistemic uncertainties, and to the definition of reliable numerical models.

1. Introduction

Existing concrete gravity dams are fundamental assets for all countries around the world because of their multiple connections with local communities for flood control, water supply and energy production [1]. Most of the existing large concrete gravity dams have been designed in the last century without considering any seismic design criteria [2]. However, concrete dams have behaved well when subjected to earthquakes. They have never collapsed after experiencing seismic events [3] and only few of them have exhibited damage.

Despite the good seismic performance shown by existing concrete gravity dams, their seismic safety must be quantified in order to allow decision-makers to prioritize Disaster Risk Reduction (DRR) strategies, thus increasing the resilience of communities. This justifies why professionals and scientists are still investing resources in this field [4].

Two different approaches are commonly adopted to assess the seismic safety of a structure: the Load-and-Resistance-Factor Design (LRFD) [5] and the Performance-Based Earthquake Engineering (PBEE) [6,7]. In the LRFD approach, the structural safety is primarily assessed in terms of failure of individual building components. The resistance of the component under investigation is compared to the effects of the actions affecting it in a particular scenario, the so-called Limit State (LS).

The LRFD approach is the basis of numerous modern codes [8] because it allows standardising design and safety assessment of a large number of new and existing structural typologies. Code parameters, such as combination factors, characteristic values, confidence levels and safety factors, are needed to make the LRFD approach applicable to real cases. These factors allow matching a target structural reliability in a particular scenario and accounting for all the uncertainties involved in the analysis. The calibration of those parameters can be based on: expert judgment, fitting methods, code optimization or a combination of these [9].

Expert judgment was the common approach until computational tools were developed to solve complex optimization problems. In that case, the parameters of the code were initially guessed based on expert knowledge and then verified over the years.

When expert judgment cannot be derived, for instance, because of lack of historical data, and the construction features of the analysed structural typology do not allow technicians to define clusters of similar buildings characterised by common failure modes, then the application of LRFD approach is vague, as in the case of concrete gravity dams.

In the PBEE approach, structural performances are assessed at system levels in terms of risk of collapse, fatalities, direct and indirect costs; this is usually known as “dollars, deaths and downtime” [10]. The outputs are then probabilistic performance metrics which are relevant to decision-makers for the definition of decision criteria. The PBEE aims to estimate the frequency of threshold exceedance by a particular performance metric for a given design at a given location.

The PBEE framework, in the version proposed by Pacific Earthquake Engineering Research Center (PEER) [10], can be thought as being composed of four steps: hazard analysis, structural analysis, damage analysis and loss analysis. The result of each step is then combined to obtain a probabilistic description of the Decision Variable ($DV$).

The hazard analysis makes use of the Probabilistic Seismic Hazard Analysis (PSHA) for the selection of ground motions to describe the annual frequency of exceedance of a seismic excitation, for a given geographic area. This step requires the choice of an Intensity Measure ($IM$) for the seismic excitation. The $IM$ has to be a quantity that captures attributes of the ground motion hazard at a specific site. Simple scalar $IM$s which are readily available from the seismic hazard analysis, such as the Peak Ground Acceleration (PGA), may introduce a broad variability in structural analysis results, thus requiring a large number of nonlinear analyses in further steps. The choice of a proper $IM$ should involve considerations related to both site characteristics and dynamic properties of the analysed structure [11]. The selection of the best $IM$ for the seismic assessment of existing concrete gravity dams is still an important open issue [12].

During the structural analysis step, numerical models of the facilities are built to estimate the uncertain structural behaviour, which is synthetically expressed in terms of Engineering Demand Parameters ($EDP$s). $EDP$s must be selected based on the performance target and the type of system of interest in order to be relevant. In this context, there are fundamental differences in tracking damage or collapse LSs. Assuming that the computational model is refined enough to be able to properly describe the material degradation due to seismic effects, the collapse is implicitly considered within the simulations and collapse LSs become just a binary limit point. Whereas, pre-collapse damage accumulates gradually and requires further steps to relate $EDP$s to appropriate Damage Measures ($DM$s). Moreover, this aspect is an important open issue in dam engineering.

Unlike other structure typologies, dams are characterised by construction features which make the definition of dam-specific $EDP$-to-$DM$ relationships particularly complex. In fact, the great number of different dam features and the variety of the ancillary works’ arrangement prevent the generalization of the $EDP$-$DM$ relationships.

The loss analysis, the last stage of the PBEE process, can be seen as the probabilistic estimation of structural performances parameterized through Decision Variables ($DV$s), such as direct economic losses, down-time and life-safety. When the analysed structure is characterised by a high collapse risk, $DV$s will reflect full replacement costs, rebuild time and risk of causalities associated with the collapse. The calculation of the relationship between collapse and direct losses is fairly direct, while the quantification of causalities requires further considerations related to the number of people who will suffer due to the collapse, the rescue resources, the importance of the assets, etc.

Assuming that $p\left[X|Y\right]$ is the probability density of X, conditioned to the event Y, and $g\left[X|Y\right]$ is the frequency of occurrence of X, given Y, the PBEE framework can be represented as in Figure 1.

Whereas, if D is the facility definition, the frequency of occurrence of $DV$ for a given structure can be quantified as,

$$g[DV|D]=\int \int \int p[DV|DM,D]p[DM|EDP,D]p[EDP|IM,D]p[IM|D]dIMdEDPdDM$$

(1)

Equation (1) directly follows the total probability theorem, where uncertainties in each step of the process are described in terms of independent conditional probabilities.

$DV$s are usually controlled by less severe and more frequent earthquakes for which the collapse risk is low and the progression of damage is more gradual. In this case, the aforementioned problems related to the estimation of indirect costs still remain but also the direct loss quantification is more complicated because of the uncertainties involved in loss predictive models and in the definition of $DM$-to-$DV$ relationships for a given facility. This is exactly the case of existing concrete dams, for which the definition of reliable $DM$-to-$DV$ relationships is even more complicated by the specific use of the analysed dam. A closer collaboration among scientists, public agencies, dam owners and (re)insurance companies could help to solve these issues, thus making the first step towards the application of PBEE in dam field.

A more accurate estimation of the dam safety and a more comprehensive treatment of the involved uncertainties may be achieved by moving from the LRFD approach to the PBEE one. This is even more important considering that national and international codes relate concrete dam failure to unsatisfactory performances that could lead to uncontrolled water release [13]. The limited availability and poor quality, in terms of data aggregation, of historical damage data associated with several seismic-prone areas make the derivation of numerical fragility functions (i.e., based on computational models of the structures) an essential component of probabilistic seismic risk assessment for dams. These needs have been already highlighted by Hariri-Ardebili in their pioneering work [14], which paved the way for the application of PBEE to dam engineering, addressing some of the main issues related to the dam fragility analysis.

The main focus of this paper was to study how epistemic uncertainties affect the fragility analysis of existing concrete gravity dams. Starting from an extended literature review on the fragility curve derivation for concrete gravity dams, this study presents a classification of the main sources of uncertainties involved in the seismic analysis of such infrastructures. Finally, an illustrative case of study is presented to show how epistemic uncertainties affect the calculation of the seismic fragility, thus highlighting the need for considering them in the application of the PBEE to existing concrete gravity dams.

2. Source of Uncertainties Involved in the Seismic Analysis of Existing Concrete Gravity Dams

Uncertainties in the field of numerical modelling can be divided into two categories: epistemic and aleatory [15]. The term aleatory, which is derived from the Latin alea, means the rolling of the dice and it represents the intrinsic randomness of a physical phenomenon. The word epistemic derives from the Greek $ϵ\pi \iota \sigma \tau \eta \mu \eta$ (episteme), which means knowledge. An epistemic uncertainty is then related to the lack of knowledge. Pragmatically, the main difference between them is that epistemic uncertainties can be reduced while aleatory uncertainties cannot. In practical cases, uncertainty reduction process involves the subjectivity of the analyst who decides which uncertainties of his/her model are random and which epistemic, based on the problem typology, his/her own experience and the available information.

In the context of seismic assessment of existing structures, aleatory uncertainties are related to the variability of ground motions (the so called record-to-record variability), while epistemic uncertainties are mostly related to the variability of the mechanical characteristics of materials [16].

Numerical models simulating the behaviour of complex systems are sources of uncertainty. In this context, the choice of the best modelling approach is known as model class selection [17]. The selection of the best modelling approach is an aspect that must be considered in the seismic assessment of existing concrete gravity dams [18]. De Falco et al. [19,20] address this topic highlighting how different modelling approaches involve the variation of analysis results. The authors mainly focus on three modelling aspects for gravity dams: geometrical modelling, soil–Structure Interaction (SSI) and fluid–Structure Interaction (FSI). The geometrical modelling approach is related to the choice of 2-dimensional (2D) or 3-dimensional (3D) numerical models. De Falco et al. [21] show that only 3-dimensional models are able to properly describe the seismic behaviour of the dam-basin-soil system. From the perspective of the dam body, 3D models allow the peculiarities of each dam to be described, thus catching failure modes which cannot be reproduced in 2D models. Three-dimensional modelling also affect the FSI and SSI results. Regarding the FSI approach, the authors show the main differences between the Westergaard added mass [22] and the acoustic Finite Element (FE) modelling [19]. The comparison between simple parameterized models representing prototype dams and complex models of real dams shows that the two approaches lead to the same results only around the fundamental period of the system, while they provide completely different results for higher frequency values.

Finally, the authors investigated the effects of the SSI modelling approach by comparing several commonly adopted strategies: the massless soil approach, the spring and dashpot boundary conditions, the Perfectly Matched Layer (PML) [23] and the Infinite Elements (IEs) [24]. The results reveal that the inertial part of the SSI, the so-called radiation damping [25], is a fundamental aspect which cannot be neglected in the simulation of dam seismic behaviour. Therefore, the massless approach turns out to be inadequate for this purpose. It is worth noting that no one has ever really proven that the massless approach is fully conservative in the seismic analysis of dams.

Once a particular deterministic model, i.e., a model class $\mathcal{C}$, has been defined, epistemic uncertainties related to numerical model parameters arise. They produce a variability of the model output leading to a mismatching between results and real observations [26]. The discussion above suggests a natural division of epistemic uncertainties into two categories based on their origin: model class and model parameter uncertainties, (Figure 2).

Epistemic uncertainties should be reduced as much as possible in order to provide a reliable prediction of the system behaviour [27]. The process of reducing uncertainty is generally known as model parameter calibration, and it is performed by observing the actual behaviour of the system. Being an ill-posed inverse problem in the Hadamard’ sense [28], the model parameter calibration requires a regularisation to be solved.

Both deterministic and probabilistic procedures can be adopted for this purpose. Deterministic procedures [29] aim to minimize objective functions expressing the relationship between observations and model output. Probabilistic procedures, such as the Bayesian inference [30], exploit prior information about the uncertain parameters (i.e., prior distributions) to regularise the inverse problem. Deterministic procedures require low computational performances and they enable deterministic error estimators related to the quality of the fitting to be calculated. On the other hand, probabilistic approaches are more computationally demanding than deterministic ones, but they allow deriving more comprehensive estimations of model parameters, i.e., Probability Density Functions (PDF) of model parameters and errors. The comparison between posterior distributions of errors related to different modelling approaches is the basic tools for the model class selection.

In this context, Sevieri et al. proposed a Bayesian approach for the reduction of the uncertainty related to elastic parameters of the materials based on static [31] and dynamic [26] measurements recorded by dam monitoring system. The authors proposed probabilistic hybrid predictive models for static and dynamic Quantity of Interest (QIs), such as displacements, frequencies, mode shapes, in which the dam behaviour is represented by the results of FE models. The computational burden is strongly reduced by using the general Polynomial Chaos Expansion technique (gPCE) [32] to surrogate the numerical results. These important examples show the possibility to apply probabilistic approaches to real dams.

3. Seismic Fragility Assessment of Existing Concrete Gravity Dams

Fragility functions are an important tool in earthquake engineering to compute the probabilities of different damage states as a function of seismic response of every type of structures and infrastructures. Originally created to characterize the seismic risk of nuclear power plants [33], fragility curves were developed for large systems such as buildings and bridges (e.g., [34]) and also for individual components (e.g., [35]).

Several recent studies focus on the seismic fragility assessment of existing concrete gravity dams. Most of them address the problem by analysing simplified 2D FE models. The reason for that is the massive computational burden related to the use of 3D models within probabilistic procedures, such as uncertainty propagation techniques. The analysis of refined numerical models with a large number of dynamic degrees-of-freedom is prohibitive without high performance computing.

Some studies assume nonlinear material models, otherwise the nonlinearities are concentrated within the interfaces, be they lift or contraction joints. The aleatory uncertainty is always related only to the variability of ground motions, while epistemic uncertainties are related to constitutive model parameters, basin level, effective uplift area and drain/grout curtain effectiveness. In the absence of case histories and precise indication about LSs in the regulation codes, the definition of dam failure is different from one author to another. LSs are usually identified as scenarios that are most likely based on common sense, their definition is then arbitrary and variable. They generally represent damage states that are easy to manage from a computational point of view, as the cracking of a certain section (e.g., dam base, neck), the attainment of a limit deformation or strength, or the limit extent of damaged areas. Thresholds are set in the form of limit stress ranges, limit crest-to-hell displacements or maximum damaged areas, mostly related to the overall size of the dam. In conclusion, the available studies are more focused on the choice of analysis method, material models, selection of uncertain parameters and methods for fragility curve derivation rather than on the definition of proper LSs.

One of the first probabilistic seismic assessment of existing concrete gravity dams is performed by de Araújo, J.M. and Awruch, A.M. [36], who combine aleatory uncertainties related to the record-to-record variability of ground motions and epistemic uncertainties related to the materials characteristics, by using a Monte Carlo simulation. The computational burden is reduced by modelling seismic excitation as a non-stationary stochastic process having two basic random variables, the acceleration amplitude and the phase angle. The dam concrete compressive strength is modelled as random variable, while other material properties, such as tensile strength, Young’s modulus and adhesion of the dam-foundation interface are deduced through deterministic equations. Instead of calculating proper fragility curves, the authors derive the Cumulative Distribution Function (CDF) of safety factors against slide at the base, concrete crushing at the toe and tensile cracking at the heel for 50 simulations.

Tekie and Ellingwood [37] perform a probabilistic analysis of the Bluestone gravity dam in West Virginia (US). In this study, the aleatory uncertainties are only related to the ground motion. The deterministic numerical model is composed of rigid blocks with frictional interfaces, while dam-foundation interface is characterized by a perfectly plastic Mohr-Coulomb law. The epistemic uncertainties are related to the mechanical properties of the concrete and foundation, the drain and the grout curtain effectiveness, the water elevation, the effective uplift area. The effects of epistemic and aleatory uncertainties are assessed through a Monte Carlo simulation using the Latin Hypercube Sampling (LHS) method. Different LSs are considered in this work, (1) cracking of the neck; (2) compressive failure of the foundation material at the toe; (3) sliding at the dam-foundation interface and (4) deflection of the crest to the heel, that was limited to a percentage of the dam height. Fragility curves are log-normally distributed for different LSs.

Lin and Adams [38] present a set of seismic fragility curves based on expert judgement for concrete gravity dams located in eastern and western Canada. Two-dimensional dam models are used to derive fragility curves in terms of damage state according to [39]. Log-normal CDFs are then used to fit discrete points representing damage states.

Mirzahosseinkashani and Ghaemian [40] present a study on the Pine Flat dam, where the only record-to-record variability is considered, while epistemic uncertainties are neglected. The soil in the FE model is elastic and massless and the concrete of the dam body is simulated via a smeared crack model. Two LSs are considered, the one is attained for a given crack length at the base of the dam and the other for a given amount of the cracked area on the dam faces.

Lupoi and Callari [41] present a new probabilistic procedure for fragility assessment of existing concrete dams by considering epistemic and aleatory uncertainties separately. The LS function is linearized with respect to the model parameters (i.e., random variables). The reliability problem is thus characterized by a general cut-set formulation, and the failure probability is calculated by using a Monte Carlo simulation. The authors apply the proposed procedure to the case of the Kasho gravity dam, using an elastic FE model simulating the full fluid structure interaction. Only operational LSs are considered: (1) crest-to-heel displacements; (2) cracking at the neck; (3) base sliding governed by tangential stress and (4) cracking at the upstream face. Epistemic uncertainties are related to only two parameters: the concrete strength and the soil strength. Two aleatory uncertainties are considered: the seismic action and the water level.

Hebbouche et al. [42] adopt the procedure proposed by Tekie and Ellingwood to consider both epistemic and aleatory uncertainties in the fragility calculation. The authors assume a linear behaviour for the dam body, while the soil is modelled as perfectly plastic material with Mohr-Coulomb yield criterion. The dam–soil interface is governed by Coulomb’s friction law. The LHS method is used to sample from six uncorrelated random variables: friction angle, cohesion, foundation dilation angle, Young’s moduli of concrete and soil and concrete compressive strength. Four LSs are considered, (1) cracking at the neck; (2) sliding at the dam-foundation interface; (3) crest to heel displacement and (4) crushing at the toe. Six near-fault ground motions are used and scaled with regard to the pseudo-spectral acceleration (PSA) from 0.2 to 2 g.

Ghanaat et al. [43] provide fragility curves for different concrete gravity dams. Both material parameters and seismic input are assumed as random variables and the sampling approach is based on the LHS method. The authors analyse the seismic performance of the Mühleberg’s gravity dam using a detailed 3D FE model having only a nonlinear contact surface between dam and foundation. The following uncertainties are considered: concrete elastic modulus, concrete damping, rock elastic modulus, cohesion and friction angle. They perform Incremental Dynamic Analysis (IDA) [44] with 30 three-components ground motion records. Given the high computational cost, each analysis is stopped when the first convergence failure occurs. In this way, any possible “resurrection” [44] at higher intensity is avoided. Finally, the fragility curve is obtained by fitting a log-normal CDF through the data points, via the least-square approach.

In Ghanaat et al. [45] a simplified version of the previous method is presented. In this case the number of trial analyses for the combination of epistemic and aleatory uncertainties is smaller with respect to the previous study. A different scaling approach and different directional factors are also used. Fragility curves are then calculated by considering a Weibull distribution rather than a log-normal one. The proposed procedure is applied to the highest non-overflow monolith of a gravity dam with shape similar to that of the Folsom dam. The 3D model has a nonlinear dam-foundation interface and an additional upper lift joint at the neck. Local and global failure modes are then considered. Elastic modulus and damping coefficients of concrete and rock are treated as random variables together with tensile strength, friction angle and cohesion of the nonlinear joints.

Finally, in Ghanaat et al. [46], a 3D FE model with linear elastic materials, massless soil, contraction and peripheral joints with nonlinear behaviour, is adopted for the seismic fragility analysis of a concrete dam. Thirty ground motions are selected and scaled in order to achieve five LSs. Each analysis is performed twice, treating the parameters both as random variables and as deterministic. The model parameters are: elastic modulus, compressive and tensile strength of concrete, maximum aggregate size, elastic modulus and tensile strength of rock, cohesion and friction angle of the interfaces.

Kadkhodayan et al. [47] analyse an arch gravity dam through nonlinear IDA, considering the nonlinearity concentrated within the contraction joints. The selected DM is the percentage of overstressed area on the dam faces, initially proposed by Ghanaat [48]. A linear relationship between $DM$ and $IM$ is assumed. Nine three-components ground motion records are selected and three potential $IM$s are considered: PGA, Peak Ground Velocity (PGV) and Spectral Acceleration at the structure’s first-mode period (${\mathrm{S}}_{\mathrm{a}}\left({\mathrm{T}}_1\right)$). Three LSs are then defined based on different damage levels.

Bernier et al. [49,50] derive seismic fragility curves for concrete gravity dams through the analysis of a 3D FE model with nonlinearity concentrated within the dam-foundation interface and the lift-joint at the neck of the dam. The parameters of the nonlinear interfaces, cohesion, tensile strength, friction angle and damping ratio are treated as epistemic uncertainties, whereas, the seismic variability is the only aleatory uncertainty. The vertical component of the ground motion is also considered, by scaling the horizontal component with a random factor ranging between 0.5 and 0.8. Two LSs are assumed: sliding at the dam-foundation interface and sliding at the lift–joint interface. After the seismic analyses, the authors compare three distributional models, e.g., normal, log-normal and Weibull, applying two fitting methods: the Sum of Squares due to Error (SSE) and the Maximum Likelihood Estimator (MLE). Finally, the effect of the spatial variation of the friction angle in the dam-foundation interface is investigated.

Hariri-Ardebili et al. [51] assess the seismic fragility of concrete gravity dams via the Multiple Stripe Analysis (MSA). Nine ground motions are selected and three different intensity levels are considered for each of them, thus resulting in 27 transient analyses. The authors perform a set of linear elastic analyses and a set of nonlinear analyses, where the nonlinearity is modelled by rotating smeared crack and Mohr-Coulomb base joints. Only the record-to-record variability is considered. The fragility curves are derived considering the following $DM$s, in the case of linear analyses: Demand Capacity Ratio (DCR) [52], Cumulative Inelastic Duration (CID) [53], Cumulative Inelastic Area (CIA) and Damage Spatial Distribution Ratio (DSDR). In the case of nonlinear analyses the following DMs are considered: Joint opening damage index, Joint sliding damage index and Crack-based damage index.

In Hariri-Ardebili and Saouma [54] the authors use the Cloud Analysis (CLA) [55] for the seismic assessment of concrete gravity dams. The authors analyse the tallest non-overflow monolith of a 122 m high dam via a 2D model, in which the interface between dam and foundation is the only source of nonlinearity. Uplift pressures are then automatically adjusted with regard to the crack lengths. Epistemic uncertainties related to the model parameters are neglected, but a large set of ground motions of 100 records is considered. Moreover, fragility curves and surfaces are developed for 70 different $IM$s in order to compare them and to identify the one with minimum dispersion.

In Hariri-Ardebili and Saouma [12] a numerical model similar to the one developed in Hariri-Ardebili and Saouma [54] is considered, this time adopting the concrete smeared crack model. Two different ground motions combinations are used, horizontal component only or horizontal and vertical components, in two different loading scenarios, empty and full reservoir. Twenty-one earthquake records are applied for each combination of loading scenarios and ground motion assumptions. For each case a log-normal CDF is fitted through the model outputs.

In Hariri-Ardebili and Saouma [56] the Endurance Time Analysis (ETA) method is adopted for the derivation of the fragility curve. Since only one ETA is performed, the aleatory uncertainty related to the record-to-record variation is neglected. A 2D model with nonlinearity concentrated within the interface between dam and foundation soil is used in this work. The elastic parameters of concrete and rock, and those of the lift joints, including the specific fracture energy, are treated as random variables. Two probabilistic Monte Carlo simulations, using LHS sampling, are applied to propagate the uncertainties through the numerical model. In one case the random variables are considered correlated, while in the second case they are considered uncorrelated. Log-normal CDFs are finally fitted through the deterministic results by using the MLE approach to derive the fragility curves.

In Hariri-Ardebili et al. [57] the Koyna dam is analysed by means of a 2D elastic FE model with massless soil. In a first analysis the dam concrete is assumed homogeneous, so the elastic parameters are modelled as random variables. Whereas, in a second analysis the dam concrete is assumed heterogeneous, and the elastic parameters are modelled as random fields. The structural safety is evaluated during a post-process step, in which several performance indices are considered. The authors vary the correlation length of the random field showing the effect on the model output in terms of crest and neck displacement, principal stress at the base and performance indices.

4. Illustrative Application: The Effect of Epistemic Uncertainties on the Seismic Fragility Analysis of Existing Concrete Gravity Dams

4.1. The Dam

An existing Italian large concrete gravity dam is selected as case study to show how epistemic uncertainties related to material mechanical parameters affect the seismic fragility analysis. The dam is 55.5 m tall and 199.20 m in length. As no vertical contraction joints are present and the dam is curved in plan, a 3D FE model is built to properly simulate the structural dam behaviour.

4.2. The Numerical Model

The original technical drawings and the orographic map of the region are elaborated to build a 3D Computer-Aided Drawing (CAD) of the dam-reservoir-soil system (Figure 3). ABAQUS 6.14 [58] is then used to perform FE analyses.

The mesh of the FE model is made up of 13,280 quadratic tetrahedral C3D10 elements for the dam body, while 7671 linear tetrahedral AC3D4 acoustic elements simulate the reservoir at its maximum level (52.9 m). As discussed in Section 2, both SSI and FSI must be modelled to achieve a proper refinement and subsequent accuracy of the solution. Therefore, 237 linear hexahedral CIN3D8 infinite elements reproduce the unboundedness of the soil domain, thus considering the radiation damping, while low reflecting boundary conditions are applied at the boundary of the acoustic domain, thus avoiding wave reflections. It is worth noting that the minimum mesh size is smaller than one tenth of the smallest wavelength of the seismic input, considering a significant frequency range of the seismic motion.

The model reproduces the foundation soil which is composed of two geological formations, an arenaceous mass and a marl mass, which have completely different mechanical characteristics.

Only the collapse LS of the dam body is considered. Special attention is then paid to the modelling of the dam concrete and its parametrization for the solution of the forward problem.

The concrete damage plasticity constitutive model [59] is suitable to simulate the behaviour of the dam concrete under cyclic loadings. In addition to the elastic parameters, this model requires the definition of plastic flow parameters and the concrete post-elastic behaviour both in compression and in tension. Plastic flow parameters commonly available in the scientific literature [60] were used in this study (

Table 1. Plastic flow parameters.

$\mathit{\psi }$ [°]	$\mathit{\epsilon }$	${\mathit{\sigma }}_{\mathit{b}0}/{\mathit{\sigma }}_{\mathit{c}0}$	$\mathsf{\Lambda }$
36.31	0.1	1.16	0.66

). In particular, $\psi$ is the dilation angle, $\epsilon$ is the eccentricity that defines the rate at which the flow potential approaches the asymptote, ${\sigma }_{b0}/{\sigma }_{c0}$ is the ratio of initial equibiaxial compressive yield stress to initial uniaxial compressive yield stress and $\mathsf{\Lambda }$ is the ratio of the second stress invariant on the tensile meridian to that on the compressive meridian.

The post-elastic compressive behaviour (Figure 4a) is represented by three different branches: the hardening branch, the first softening branch and the second softening branch. Let us consider the concrete compressive strength $f_{C,c}$, the undamaged concrete elastic moduli $E_C$, the compressive stress value beyond which the material shows a non-linear behaviour $f_{C,c0}=0.85f_{C,c}$ and the compressive stress value beyond which the material strength exponentially decreases $f_{C,cu}=0.3f_{C,c}$ (these two latter correspond to ${ϵ}_{C,c0}=f_{C,c0}/E_C$ and ${ϵ}_{C,cu}=0.005$ respectively). The post-elastic compressive behaviour can be then described in terms of stress ${\sigma }_c$ and strain $ϵ$ by the following equations:

• Hardening branch: if ${ϵ}_{C,c0}<{ϵ}_c\le {ϵ}_{C,c}$,

  $${\sigma }_c({ϵ}_c)=f_{C,c0}+(f_{C,c}-f_{C,c0})\sqrt{\frac{2({ϵ}_c-{ϵ}_{C,c0})}{{ϵ}_{C,c}-{ϵ}_{C,c0}}-\frac{{({ϵ}_c-{ϵ}_{C,c0})}^2}{{({ϵ}_{C,c}-{ϵ}_{C,c0})}^2}}.$$

  (2)
• First softening branch: if ${ϵ}_{C,c}<{ϵ}_c\le {ϵ}_{C,cu}$,

  $${\sigma }_c({ϵ}_c)=f_{C,c}-(f_{C,c}-f_{C,cu}){\left(\frac{{ϵ}_c-{ϵ}_{C,c}}{{ϵ}_{C,cu}-{ϵ}_{C,c}}\right)}^2.$$

  (3)
• Second softening branch: if ${ϵ}_c>{ϵ}_{C,cu}$,

  $${\sigma }_c({ϵ}_c)=f_{C,cr}+(f_{C,cu}-f_{C,cr})exp\left(2\frac{f_{C,cu}-f_{C,c}}{{ϵ}_{C,cu}-{ϵ}_{C,c}}\frac{{ϵ}_c-{ϵ}_{C,cu}}{f_{C,cu}-f_{C,cr}}\right).$$

  (4)

where $f_{C,cr}$ is the residual value of the compressive strength, herein assumed equal to $0.1f_{C,c}$ and ${ϵ}_{C,c}$ is the strain value related to $f_{C,c}$. Once $f_{C,c}$ and $E_C$ are sampled from their PDF, the compressive post-elastic behaviour is completely defined.

The parametrization of the post-elastic tensile behaviour (Figure 4b) of the concrete requires more attention. According to Hillerborg et al. [61] the fracture energy $G_f$ is the energy required to open a unit area of crack and it can be thought related to the area under the post-elastic stress-fracture curve. Therefore, defining the tensile constitutive law by the stress ${\sigma }_t$ and displacement $u_t$ (or cracking displacements $u_t^{ck}$ if one considers only post-elastic displacements) is a particularly convenient choice to reduce the mesh dependence of the problem solution. In this study, the tensile post-elastic behaviour of the concrete is described by the law

$${\sigma }_t\left(u_t^{ck}\right)={\alpha }_{t}exp(-{\beta }_t{u}_t^{ck}).$$

(5)

where ${\alpha }_t$ and ${\beta }_t$ are the shape parameters of the exponential function. Once a value of the tensile strength $f_{C,t}$ is sampled from their distributions, ${\alpha }_t$ and ${\beta }_t$ are calculated in order to have $G_f$ equal to 150 N/m. This value for $G_f$ is commonly adopted in dam engineering [62]. The calibrated post-elastic tensile behaviour is finally used to determine the tensile damage law (Figure 4c) which describes the variation of the tensile damage variable $d_t$ versus $u_t^{ck}$.

The reduction of the concrete Young moduli due to damage related to the development of micro-cracks is expressed by a scalar value $d_t$. The tensile damage law is then simulated by a second order polynomial curve (Equation (6)) and its parameters are calibrated by assuming that the related plastic displacements $u_t^{pl}$ (Equation (7)) must be positive and descending,

$$d_t\left(u_t^{ck}\right)=a_{d_t}+b_{d_t}{u}_t^{ck}+c_{d_t}{u_t^{ck}}^2,$$

(6)

$$u_t^{pl}=u_t^{ck}-\frac{d_t}{(1-d_t)}\frac{{\sigma }_t{l}_0}{E_0}.$$

(7)

A proper Matlab [63] code was developed to automatize the calculation of the tensile behaviour parameters for every sampled tensile strength value.

The values of mechanical characteristics of the dam concrete and the soil are deduced from the experimental campaigns conducted over the years. The results of in situ tests on dam concrete and foundation soil performed in 2012 are reported in

Table 2. Material test results.

	${\mathit{E}}_{\mathit{C}}$	${\mathit{\upsilon }}_{\mathit{C}}$	${\mathit{\rho }}_{\mathit{C}}$	${\mathit{E}}_{\mathit{A}}$	${\mathit{\upsilon }}_{\mathit{A}}$	${\mathit{\rho }}_{\mathit{A}}$	${\mathit{E}}_{\mathit{M}}$	${\mathit{\upsilon }}_{\mathit{M}}$	${\mathit{\rho }}_{\mathit{M}}$	${\mathit{f}}_{\mathit{C},\mathit{c}}$	${\mathit{f}}_{\mathit{C},\mathit{t}}$
	[Mpa]		[kg/${\mathbf{m}}^3$]	[Mpa]		[kg/${\mathbf{m}}^3$]	[Mpa]		[kg/${\mathbf{m}}^3$]	[Mpa]	[Mpa]
Mean	34,620	0.16	2270	21,000	0.16	1800	7000	0.22	1800	16.7	2
s.d.	12,498	0.06	189.17	9545	0.07	-	3181.8	0.1	-	5.16	1.1

in the form of mean values and the standard deviations (s.d.). The Young modulus E, the Poisson’s ratio $\upsilon$ the density $\rho$, the compressive and the tensile strengths, $f_c$ and $f_t$ respectively of concrete, arenaceous mass (orange in Figure 3) and marl mass (green in Figure 3), indicated with subscripts C, A and M, respectively, are reported below.

4.3. Uncertainty Quantification (UQ) in the Modal Analysis of Concrete Dams

In this section, the epistemic uncertainties related to the elastic parameters of the materials are propagated through the numerical model to evaluate how they affect the dam dynamic behaviour in terms of modal analysis (sensitivity analysis).

Mean values and standard deviations of the material Young moduli E and Poisson’s ratios $\upsilon$, derived during the 2012 experimental campaign (Table 2), are used to define the PDF of the bulk K and shear G moduli of concrete (i.e., $K_C$ and $G_C$, respectively) as well as those of arenaceous mass (i.e., $K_A$ and $G_A$) and marl mass (i.e., $K_M$ and $G_M$). These random variables are assumed log-normally distributed in the present study.

The gPCE belongs to the family of the spectral methods for the propagation of uncertainties through deterministic models [32]. Let us assume that the uncertain model parameters are collected in $\mathit{\theta }$, in this approach a polynomial expansion ${\widehat{u}}_P^{\mathrm{gPCE}}$, also called response surface, of the uncertain model output $u(\mathit{\theta })$ is provided by using orthogonal basis functions ${\mathsf{\Phi }}_{\mathbf{i}}(\mathit{\theta })$, that is,

$$u(\mathit{\theta })\approx {\widehat{u}}_P^{\mathrm{gPCE}}(\mathit{\theta })=\sum _{\left|\mathbf{i}\right|\le P}{\mathit{u}}_{\mathbf{i}}{\mathsf{\Phi }}_{\mathbf{i}}(\mathit{\theta })$$

(8)

where P is the degree of the polynomial expansion, $\mathbf{i}$ is a finite multi-index set and ${\mathit{u}}_{\mathbf{i}}$ is the matrix of the combination coefficients. The selection of the orthogonal basis functions is based on the PDF of the random parameters $\mathit{\theta }$. Whereas, the combination coefficients are calculated with regard to some reference solutions of the deterministic model by applying interpolation, regression or the Bayes rule. The Bayesian approach for the gPCE combination coefficient derivation proposed by Rosić and Matthies [64] is adopted in this study. Once a reliable gPCE is provided, the forward problem (i.e., the calculation of the model output statistics) is easily solved by processing the combination coefficients collected in ${\mathit{u}}_{\mathbf{i}}$. In addition, the variance-based sensitivity analysis (Sobol’ method, [65]) can be performed without any addition computational cost [66]. As already mentioned, the gPCE has been already applied in dam engineering both to propagate uncertainties and to build physic-based predictive models for the static and dynamic structural control [26,31].

The presence of the SSI in the FE model leads to a large number of numerical modes related to the soil motion only, while the propagation of the uncertainties changes the relative positions of modes once a set of new parameters is considered. These two issues are part of the so-called mode matching problem [26]. Being interested in the quantification of the uncertainty for specific dam modes, let us say those which mobilise the largest amount of dam mass, a selection criterion is needed to discard modes related to the soil only. The procedure proposed by Sevieri [60] is used to this purpose. Figure 5 shows the distributions of the frequencies related to the first three modes of the system which mobilise the highest amount of dam mass (Figure 6). The results show that the variation of the elastic parameters leads to a variability in the first three frequencies, which slightly increases toward higher modes. The mean values are 6.32 Hz, 7.8 Hz and 9.2 Hz, while the Coefficient of Variations (CoVs) are 0.1, 0.12 and 0.15, respectively.

The influence of each random parameter on the first seventeen frequencies of the system is evaluated through the Sobol’ Indices [65]. Figure 7 shows that each random parameter influences the result of the modal analysis. Therefore, none of them can be neglected in the seismic fragility analysis.

4.4. Fragility Analysis

A preliminary study allows the collapse LS to be calibrated. The LS is expressed as the attainment of a limit value of the difference between crest displacements recorded at the middle of the dam length and the corresponding point at the heel., $\mathsf{\Delta }{\delta }^{\mathrm{top}-\mathrm{heel}}$. As it was demonstrated that in this case displacements greater than $0.003\%$ of the dam height lead to a diffused damage of the dam body, this value can be related to the uncontrolled water release, and can also be associated with the achievement of the collapse LS (Equation (9)).

$$\mathsf{\Delta }{\delta }^{\mathrm{top}-\mathrm{heel}}\ge 16.65\mathrm{mm}.$$

(9)

In addition to the elastic parameters of materials, whose importance is highlighted in Section 4.3, concrete strength parameters are also modelled as random variables in the fragility analysis. They are considered normally distributed with statistics derived from Table 2.

Considering both epistemic and aleatory uncertainties in the fragility assessment of civil structures and infrastructures is a challenging task from the computational perspective. This is even more demanding when refined FE models are involved, thus leading to an unsustainable computational burden. Therefore, resolution strategies are needed.

In this study, the Multi-Record-IDA (MR-IDA) approach is adopted, as it enables the record-to-record variability to be considered by performing series of IDA with different selected ground motions. The Approximated Second Order Second Moment (ASOSM) resolution strategy proposed by Liel et al. [67] is used in conjunction to the MR-IDA [44] to further reduce the computational burden.

The IDA requires performing multiple nonlinear dynamic analyses of the dam FE model under a suite of ground motion records, each scaled to several levels of seismic intensity. The scaling levels are appropriately selected to study the structural behaviour from elastic to inelastic and finally to global dynamic instability, where the structure essentially experiences collapse. Tracing algorithms can be used to reduce the computational burden and to automatize the procedure. In this study the Hunt & Fill algorithm [44] is adopted.

Modelling SSI requires the seismic motion to be deconvolved before applying it to the boundary of the model. In this regard, Sooch and Bagchi [68] propose a suitably defined numerical procedure for ground motion deconvolution.

The first step of the ASOSM method requires the calculation of the fragility curve related only to the record-to-record variability, i.e., assuming the mean values of the random parameters. The effects of the epistemic uncertainties are then considered by calculating approximated gradients of the LS function for each random variable. The approximated gradients are evaluated by perturbing each random variable with plus and minus 1.7 times their standard deviation. This method also enables correlations between random parameters to be considered, so it is particularly suitable for studying the effect of a possible correlation between variables. Therefore, the correlation between concrete elastic/strength parameters is considered through the correlation matrix reported in

Table 3. Correlation Matrix of the concrete mechanical parameters.

	${\mathit{K}}_{\mathit{C}}$	${\mathit{G}}_{\mathit{C}}$	${\mathit{f}}_{\mathit{c},\mathit{C}}$	${\mathit{f}}_{\mathit{t},\mathit{C}}$
$K_C$	1
$G_C$	0	1
$f_{c,C}$	0.83	0.83	1
$f_{t,C}$	0.83	0.83	0.9	1

.

In the present study, 15 ground motions (

Table 4. Ground Motion Records.

#	Event	Year	Station ID	Mw	Epicentral Distance [km]	EC8 Site Class
1	Campano Lucano	1980	ST96	6.9	32	A
2	Friuli	1976	ST20	6.5	23	A
3	Campano Lucano	1980	ST98	6.9	25	A
4	Bingol	2003	ST539	6.3	14	A
5	South Iceland	2000	ST2558	6.4	5	A
6	Duzce 1	1999	ST1252	7.2	34	A
7	Tabas	1978	ST54	7.3	12	A
8	Umbria Marche	1997	ST238	6	21	A
9	Montenegro	1979	ST64	6.9	21	A
10	Basso Tirreno	1978	ST49	6	34	A
11	Golbasi	1986	ST161	6	29	A
12	Duzce 1	1999	ST3136	7.2	23	A
13	South Iceland	2000	ST2556	6.5	35	A
14	Izmit	1999	ST575	7.6	9	A
15	Friuli	1976	ST36	6	28	A

) are selected. For each of them the magnitude Mw [69] is higher than 6 and the site class is A.

Figure 8 shows the results of the fragility analysis first considering only the record-to-record variability (dashed line) and then also the uncorrelated (dotted line) and correlated (continuous line) epistemic uncertainties.

The mean values of the three fragility curves are very close: 0.661 g if only the record-to-record variability is considered, and 0.652 g if also epistemic uncertainties are modelled. The fragility curve CoV is equal to 0.015 when only the record-to-record variability is modelled. Considering epistemic and aleatory uncertainties together leads to a higher variability of the fragility curve, the CoV is 0.0327 for the uncorrelated case and 0.0394 for the correlated case. Therefore, the fragility curve variability is particularly sensitive to the epistemic uncertainty effect, while mean values are fairly constant. In addition, in this case study, modelling epistemic uncertainties as correlated does not really affect the fragility curve derivation.

These results clearly show that epistemic uncertainties related to the mechanical parameters of the materials cannot be neglected in the seismic fragility analysis of concrete dams.

5. Conclusions

This paper addresses fundamental issues related to the seismic risk assessment of existing concrete gravity dams.

The application of the Probabilistic Earthquake Engineering (PBEE) to existing dams is an interesting methodology which would allow overcoming the problems related to the use of the Load-and-Resistance-Factor Design approach in this field. The fragility analysis is a fundamental step for the quantification of the seismic risk associated to existing dams.

The focus of the paper is the quantification of the effect of the epistemic uncertainties related to the variability of the material mechanical parameters on the seismic fragility analysis of dams.

This paper discusses the main sources of epistemic uncertainty in dam engineering and presents a state-of-the-art review on the seismic fragility analyses of existing concrete dams. The effects of epistemic uncertainties are quantified through the analysis of a case study.

The results show that epistemic uncertainties strongly affect the simulation of dam dynamic behaviour and ultimately the quantification of its seismic performance. The comparison between fragility curves calculated under different assumptions show that epistemic uncertainties have a massive impact on the variation of the results.

Nevertheless, the application of the PBEE to real concrete dams still requires further studies aimed to solve practical issues, such as the development of damage-to-loss curves and the selection of significant damage scenarios.