Received: 28 August 2023 Revised: 31 October 2023 Accepted: 5 December 2023

DOI: 10.1002/eqe.4064

RESEARCH ARTICLE

Quantifying modeling uncertainties in seismic analysis of
dams: Insights from an international benchmark study

Mohammad Amin Hariri-Ardebili

College of Computer, Mathematical and
Natural Sciences, University of Maryland,
College Park, Maryland, USA

Correspondence
Mohammad Amin Hariri-Ardebili,
College of Computer, Mathematical and
Natural Sciences, University of Maryland,
College Park, MD 20742-3281, USA. PREP
Researcher, Earthquake Engineering
Group, National Institute of Standards
and Technology (NIST), Gaithersburg,
MD, USA.
Email: amin.hariri@nist.gov

Abstract
Advances in nonlinear dynamic analysis of dams have not completely resolved
concerns over modeling confidence and analysis accuracy. Verification and val-
idation offer accuracy assessment, but uncertainties persist during performance
evaluation due to both epistemic (modeling) and aleatory (parametric) sources.
Epistemic uncertainties arise from simplifications and modeling techniques.
This paper addresses epistemic uncertainties in a gravity dam seismic analysis
using data from the International Comnission on Large Dams (ICOLD) bench-
mark study.While the benchmark formulation included the finite elementmodel
of the dam, mechanical material properties, and dynamic loads, participants
retained the flexibility to opt for best-practice modeling assumptions, simplifica-
tions, and other specifics. Notable response variability emerged, particularly in
crack profiles and damage predictions. This study examines sources of variability,
quantifyingmodeling uncertainty for the benchmark problem.More specifically,
the modeling variability is quantified using the logarithmic standard deviation,
also known as dispersion. This metric enables its incorporation into other seis-
mic risk assessment and fragility studies. Under relatively low-intensity motion
(peak ground acceleration [PGA] of 0.18 g in this case), modeling dispersion
of 0.45, 0.30, 0.32, and 0.30 were calculated for the maximum dynamic crest
displacement, maximumhydrodynamic pressure at the heel, heel and crest max-
imum acceleration, respectively. Additionally, the dispersion of the failure PGA
was determined to be 0.7. Findings underscore the need for systematic seismic
response modeling in dam engineering to enhance prediction accuracy. A better
understanding of the sources andmagnitudes ofmodeling uncertainties can help
improve the reliability of dam seismic analysis and contribute to the development
of more effective risk assessment and mitigation strategies.

KEYWORDS
benchmark study, concrete dams, modeling variability, seismic analysis, uncertainty quantifi-
cation

This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any
medium, provided the original work is properly cited and is not used for commercial purposes.
© 2023 The Authors. Earthquake Engineering & Structural Dynamics published by John Wiley & Sons Ltd.

1168 wileyonlinelibrary.com/journal/eqe Earthquake Engng Struct Dyn. 2024;53:1168–1194.

mailto:amin.hariri@nist.gov
http://creativecommons.org/licenses/by-nc/4.0/
https://wileyonlinelibrary.com/journal/eqe
http://crossmark.crossref.org/dialog/?doi=10.1002%2Feqe.4064&domain=pdf&date_stamp=2023-12-19


HARIRI-ARDEBILI 1169

1 INTRODUCTION

The numerical simulation of engineered structures has perpetually posed challenges for analysts due to the numerous
unknown parameters inherent in themodeling process. These uncertainties encompass various aspects, includingmodel-
ing techniques of physical assets, the level of model fidelity, solution algorithms, loading definitions, material parameters,
and more. Among engineering structures, dams hold a distinctive position owing to the intricacies entailed in their
response. For instance, a concrete dam represents a massive concrete structure characterized by considerable heterogene-
ity across its various sections. Unlike building systems, the dam’s response is significantly influenced by interactions with
the foundation rock and the reservoir, rendering it a coupled system problem.1 This coupling introduces an additional
layer of complexity, as coupled system mechanics diverge from solid mechanics principles.
Within modern analysis procedures (and often in design contexts), verification, validation, and uncertainty quan-

tification (UQ) stand as foundational components. Verification pertains to the mathematical accuracy of numerical
solutions, while validation concerns the efficacy of mathematical models in emulating real-world systems.2 UQ combines
these facets with the inherent variability present in real-life systems to yield an overarching quantification of predictive
uncertainty.3,4
In the realm of probabilistic analysis for structural models, UQ involves the “systematic study and reduction of uncer-

tainties, both in computational simulations and real-world scenarios.” Its primary aim is to assess the likelihood of specific
outcomes, accounting for incomplete knowledge about certain facets of the system.
In the field of Civil Engineering, the uncertainty classification framework introduced by [5] (Figure 1) is widely

employed. This taxonomy categorizes uncertainties into abstracted uncertainties stemming from elements of the actual
or modeled system, and cognitive uncertainties arising from subjective abstractions of reality. Engineers and analysts
can conveniently classify uncertainties into abstracted and non-abstracted forms based on their expertise, background,
and general knowledge of the system. Non-abstracted sources encompass physical randomness, vagueness, human and
organizational errors, as well as information conflicts and confusion. The framework by [5] represents a classic and
comprehensive classification widely adopted in Civil Engineering applications.
In dam engineering, a prominent historical challenge over the past three decades revolved around validating and ver-

ifying numerical models utilized for simulations, with limited emphasis on UQ until recently. Finite element1 model
validation and verification have been previously explored in the context of structural health monitoring,6 and system
identification frameworks,7 encompassing both static and transient scenarios. While these methodologies offer insights
into dam seismic response, fully validating and verifying a finite element model against an actual dam’s behavior under
earthquake ground motion remains a formidable task.

F IGURE 1 Illustration of a comprehensive uncertainty classification in Civil Engineering, adapted from [5]. This taxonomy provides a
structured framework for categorizing uncertainties, encompassing various aspects within the field.

1 From here on, the term “finite element” will be used as a main substitute for “numerical modeling,” as the majority of dam structural analyses are
conducted using the finite element method, with only a few accessible studies exploring finite difference, extended finite element, mesh-free, and
other techniques.

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1170 HARIRI-ARDEBILI

Previous endeavors included shake table tests8 as well as centrifuge tests9 on small-scale dam models. In shake table
tests, reduced-scale dammodels are subjected to dynamic excitations based on similitude principles. However, these tests
don’t account for gravity effects on material behavior. Centrifuge testing, on the other hand, simulates the influence of
high-speed rotation on dams, accounting for gravity’s impact on material behavior. Dynamic centrifuge tests accurately
reproduce the prototype stress field and mimic prototype material properties.
Hence, ensuring the fidelity of a developed finite element model to actual behavior entails cross-validation against pub-

lished research, benchmark workshops, and blind prediction contests. Several instances of such cross-validation exist,
including those organized by ICOLD and national organizations like USSD. A significant example is the 15th ICOLD
International Benchmark Workshop10 and the 2018 USSD benchmark workshop.11 The author of this paper contributed
to the formulation of both workshops and summarized their outcomes in reports. This paper’s focus lies in highlighting
the significance of addressing modeling uncertainties in seismic analysis of gravity dams, utilizing data collected from the
2019 ICOLD benchmark study in Italy. It is worth emphasizing that the field of finite element model UQ in concrete dams
remains quite constrained, with only a handful of studies addressing this area.12–14 These limitations primarily arise from
the intricate nature of modeling the dam-foundation-water coupled system, a paucity of research efforts in this special-
ized field, limited benchmark studies, the computational intensity of resources required, scarcity of field and test data,
insufficient standardization, the inherent variability and complexity of modeling theories (for example, distinct concrete
damage models, diverse fluid-structure interaction (FSI) models, and various seismic wave propagation mechanisms in
foundation rock), and a dearth of comprehensive knowledge among analysts.
Themajority of publishedwork in probabilistic assessment of dams has predominantly focused onmaterial randomness

and ground motion record-to-record (RTR) variability. In contrast, other disciplines within structural engineering have
addressed the critical issue of modeling uncertainty. In a study by [15], uncertainties in construction details were incor-
porated on top of material properties randomness using a Bayesian framework. The findings revealed an increase in both
the mean and standard deviation of the demand-to-capacity ratio for the structure. Similarly,16 introduced uncertainty in
themodeling of structures to assess the collapse risk of buildings. The inclusion of modeling uncertainty demonstrated an
approximate 1.8-fold increase in the mean annual frequency of collapse compared to analyses neglecting modeling uncer-
tainty, particularly for a high-seismic site. Formodels of reinforced concrete frameswith infills,17 incorporated uncertainty
in the modeling parameters describing the nonlinear force-deformation response of equivalent infill struts. The study
showed a shift in the median story drift ratio in the range of about ±20%, while the dispersion increased by approxi-
mately 10-40% with the incorporation of modeling uncertainty. In the context of aging highway bridges,18 incorporated
a wide range of modeling-related uncertainties in fragility assessments. Their findings emphasized that the inclusion of
only groundmotion uncertainty is insufficient for this type of structure. Finally,19 quantified the propagation of modeling
parameter type uncertainty to the overall dispersion of the demand parameters. The mean ratio of the total dispersion to
only RTR dispersion was found to be 1.27.
In the context of nonlinear seismic analysis applied to concrete dams, several factors contribute to modeling-related

uncertainties. These include: (1) Variability in the chosen smeared crack model for concrete. (2) Variability in the discrete
damage model applied to concrete and rock joints. (3) Challenges in translating measured material properties into con-
stitutive models. (4) Variability in the FSI model. (5) Consideration or neglect of hydrodynamic pressure inside opened
joints/cracks. (6) Uncertainty in the formulation of wave-absorbing boundary conditions for water and foundation’s far-
end. (7) Uncertainty in the formulation of wave propagation and soil-structure interaction. (8) Uncertainty in the aging
model for concrete and its spatial distribution. (9) Consideration or neglect of pre-seismic loadings, such as stage construc-
tion, thermal loads, and long-term environmental loads. (10) Variability in different aspects of finite element modeling,
including mesh size, time integration size and type (explicit or implicit solution), element type, element aspect ratio, con-
vergence criteria and its threshold. (11) Consideration of geometric (large deformation) nonlinearity in the analysis. (12)
The influence of human factors and the absence of proper quality control measures.
Benchmark problems and blind prediction exercises are commonly employed to assess the precision of numerical sim-

ulations and to quantify the extent of uncertainties in predictions. For example, the impact of modeling uncertainty on
performance-based risk assessments has been systematically evaluated by researchers such as.20 They utilized outcomes
from the 7-story reinforced concrete building blind prediction contest to account for variations arising from software plat-
forms, solution algorithms, element types, andmodel parameter computation or selection. In the context of structural and
earthquake engineering, several benchmark exercises have been undertaken. Notable examples include blind prediction
tests for the seismic evaluation of unreinforced masonry buildings,21 shake table assessments of stone masonry building
aggregates,22 shake table testing of rocking podiums,23 blind predictions of the seismic response of woodframe houses,24
and a blind benchmark for the seismic assessment of slender masonry towers.25

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1171

This paper presents a pioneering effort in quantifying uncertainties related to modeling in nonlinear damage analysis
of gravity dams. The primary objective is to establish a benchmark problem, subsequently identifying modeling uncer-
tainties through collected information and quantifying their impact on structural response and dam failure. The ensuing
sections will briefly discuss various sources of uncertainties within the context of concrete dam engineering, followed by
an exploration of the case study and the information provided to participants. Subsequently, the assumptions made by
different teams will be reviewed, along with a summarized presentation of results from around 20 papers. The quantity
and nature of uncertainty in response will be quantified, and recommendations for future studies on nonlinear seismic
dam analysis will be presented.

2 CLASSIFICATION OF UNCERTAINTY SOURCES

The role of UQ holds paramount importance within probabilistic risk assessment frameworks,26 a significance that mag-
nifies when addressing critical structures27 and complex systems like dams. Various classifications of uncertainty sources
exist, each offering diverse perspectives. Here, several of these classifications will be examined.
One notable classification by [28] delineates uncertainty into three categories: (1) parameter uncertainty pertains to

uncertainty in quantifying a model with a predetermined functional form, (2) modeling uncertainty encompasses the
uncertainty surrounding the suitability of the structural or mathematical model’s form, and (3) completeness uncertainty
signifies the extent to which the model encompasses all phenomena relevant to the simulated system.
Another classification segregates variables into two fundamental categories: epistemic uncertainties and aleatory

uncertainties.29 Aleatory uncertainty, also referred to as objective uncertainty,30 characterizes the inherent randomness
of phenomena. In the context of seismic analysis for concrete dams, uncertainties linked to earthquake events (e.g., inten-
sity, timing, and recurrence) assume paramount importance as significant uncertain parameters.31 Aleatory uncertainty
is commonly quantified through random variables (RVs) using the mathematical framework of probability theory.32 Epis-
temic uncertainty, or subjective uncertainty, emerges from knowledge gaps. An instance is the uncertainty related to
material properties, such as concrete fracture energy. This category encompasses various sources, includingmodel and sta-
tistical uncertainty.33 For continuous RVs, epistemic uncertainty can be quantified through probability density functions
(PDFs).

2.1 Ground motion RTR variability

Ground motion RTR variability is important in seismic evaluation of dams, arising from the unpredictability of earth-
quake events. Methods like Incremental Dynamic Analysis (IDA)34–36 and Multiple Stripe Analysis (MSA)37–39 account
for RTR variability by utilizing scaled groundmotion records, while the cloud analysis typically uses the un-scaled ground
motion records.40–44 Seismic fragility functions are valuable tools to convey the outcomes of extensive probabilistic seismic
simulations,45 and enhance the understanding of dam responses under varying seismic conditions.

2.2 Material randomness

Material uncertainties significantly impact seismic analyses of dams. These uncertainties in mechanical properties of
material are commonly modeled with statistical distributions. When multiple RVs are involved, accounting for correla-
tions becomes essential. Multivariate correlations,46 like 𝜌(𝐸𝑐,𝑓𝑐)

, and spatial correlations47 are two important ingredients
to precisely quantify material uncertainty in dam engineering.

2.3 Modeling uncertainty

In the context of seismic analysis, the uncertainty in response modeling primarily stems from the inability of an idealized
numericalmodel to precisely predict the actual behavior of an engineered systemunder the influence of groundmotions.20
elucidate thatmodeling uncertainty embodies the uncertainty linked to themodel’s efficacy in accurately representing the

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1172 HARIRI-ARDEBILI

F IGURE 2 Illustrative examples representing the four distinct types of modeling uncertainties in the context of seismic response
analysis for concrete dams; Constitutive Model; Adapted from [48].

authentic structural response. However, Bradley48,49 further extends this classification into the subsequent sub-categories,
illustrated in Figure 2:
Level 1 Uncertainty: This category encompasses uncertainties in the measurement of physical quantities, such as the

tensile strength of concrete or the estimation of soil shear stiffness. Level I uncertainty hinges solely on the quality of
physical experiments and doesn’t directly impact numerical model preparation. Nonetheless, since these measurements
serve as input parameters in numerical simulations, they exert influence on the final response quantities.
Level 2 Uncertainty: In this tier, uncertainty arises from the correlation between measurable physical quantities and

the constitutivemodel parameters. This form of uncertainty has been observed by researchers like.50–53 For instance, while
concrete compressive strength is commonly measured from core samplings, several constitutive models may necessitate
knowledge of the tensile strength directly. As evident the level 2 modeling uncertainty is similar to material randomness
discussed in Section 2.2.
Level 3 Uncertainty: This level encapsulates the uncertainty inherent in selecting a constitutive model, involving

project-specific assumptions and simplifications.54,55 Constitutive models are often crafted either empirically based on
observations or theoretically derived via specific assumptions. In the context of concrete dams, a variety of options exist,
such as fixed/rotating smeared crack models,56 plasticity models,57 damage mechanics models,58 and fracture mechanics
models,59 to address concrete cracking and crushing during seismic loads.
Level 4 Uncertainty: This category involves uncertainty in the overall methodology of the idealized model. For

instance, it encompasses decisions such as employing 3D models for dams with narrow valleys versus using 2D plane
strainmodels.60 Other examples include considerations like accounting for dam-foundation or dam-reservoir interactions,
with potential fidelity differences. For dam-reservoir interaction, choices might involve usingWestergaard hydrodynamic
pressure or directly employing pressure-based fluid elements. Assumptions regarding boundary conditions, damping
formulations, and ground motion input, including wave propagation assumptions, also contribute to this category.

3 MULTIPLE UNCERTAINTY SOURCES

Uncertainty sources can be studied individually or collectively. The former approach examines independent epistemic
and aleatory uncertainties (or other classifications like modeling and RTR variability), necessitating their adept inte-
gration to portray the complete uncertainty inherent in the system. The latter approach directly encompasses the
correlations among various uncertainty sources, albeit often entailing higher computational expenses, particularly in
simulation-based methods.

3.1 Decoupled aleatory and epistemic uncertainties

It is important to recognize that within the context of structural and earthquake engineering, a prevailing understanding
suggests that the uncertainty arising from aleatory variables (particularly ground motion RTR variability) tends to
outweigh the epistemic randomness (stemming frommaterial and modeling variability). The stochastic nature of ground
motion, coupled with its inherent unpredictability, and the possibility of calibrating incomplete information about

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1173

0 1 2 3 4
0

0.2

0.4

0.6

0.8

1

(A) Combining decoupled uncertainties with confidence

interval approach

0 1 2 3 4 5 6 7
0

0.2

0.4

0.6

0.8

(B) Comparison of different scenarios in uncertainty com-

bination

F IGURE 3 Accounting for multiple sources of uncertainties in seismic assessment of structural systems.

material and modeling randomness are the fundamental reasons for such a belief. As a result, the spread of results
due to epistemic uncertainty is coupled with a model that explicitly addresses the effects of aleatory uncertainty. For
instance, consider the case of a seismic fragility curve, which primarily hinges on aleatory uncertainty but incorporates
adjustments for epistemic uncertainties.61
One approach to implementing this treatment is the confidence interval method.62,63 Initially, the cumulative distri-

bution function (CDF) of aleatory uncertainties is established using conventional techniques such as IDA or MSA. This
results in a fragility curve with a median capacity denoted as 𝜂 and an RTR dispersion designated as 𝛽𝑅. Subsequently, the
PDF corresponding to the median of epistemic uncertainties is determined. This can be derived through diverse meth-
ods like expert judgment, First-Order Second-Moment, direct Monte Carlo simulations, and others. For simplicity, let’s
assume that the RV 𝜂 follows a lognormal distribution2, with a median estimate of capacity denoted as 𝜂 and a dispersion
labeled as 𝛽𝑈 . The influence of epistemic uncertainty is incorporated by integrating the PDF with the median of the CDF.
While the 50% confidence level (CL) aligns with a scenario featuring a sole aleatory uncertainty, the median capacity with
CL confidence is defined as 𝜂𝐶𝐿 = 𝜂 exp(−𝛽𝑈.𝐾𝐶𝐿), where 𝐾𝐶𝐿 is the standardized Gaussian variate associated with the
probability 𝐶𝐿 of not being exceeded. Some of frequently used {𝐶𝐿, 𝐾𝐶𝐿} combinations are: {50%, 0.00}, {84%, 1.00}, and
{90%, 1.28}. To predict outcomes at a designated confidence level, 𝐶𝐿, the aleatory distribution must be adjusted to the
(1 − 𝐶𝐿∕100)th percentile of the epistemic distribution. Figure 3A illustrates hypothetical curves depicting 50%, 84%, and
90% confidence intervals.

3.2 Integrated aleatory and epistemic uncertainties

The integration of uncertainty sources provides a more comprehensive depiction of the system, avoiding unwarranted
assumptions in their amalgamation. Presently, a significant portion of combined uncertainty studies in dam engineering
encompass material randomness and RTR variability, or address level 1 and 2 modeling uncertainties combined with RTR
variability. However, relying solely on these simplified approaches (instead ofmore complex ones like levels 3 and 4model-
ing uncertainties) can be problematic. Such simplified approaches consider uncertainties confined to material and model
parameters within the constitutive model, disregarding the higher-fidelity uncertainties such as different constitutive
models and modeling strategies.
Illustrated in Figure 3 are hypothetical scenarios wherein various uncertainty sources are combined at varying degrees

of fidelity.49 Three cases are examined, each showcasing the distribution of engineering demand parameters (EDP),
for example, drift ratio: (1) a deterministic model with only aleatory uncertainty, (2) aleatory uncertainty integrated
with low-fidelity epistemic uncertainties (akin to levels 1 and 2 in modeling uncertainty), and (3) aleatory uncertainty
combined with both low- and high-fidelity epistemic uncertainties. While low-fidelity epistemic uncertainties may
marginally augment system uncertainty (owing to model parameter uncertainty in simplified numerical models and
predefined failure modes), the inclusion of high-fidelity sources of uncertainties (like constitutive model parameter

2 The lognormal distribution is chosen due to its suitability for modeling positive RVs, a characteristic advantageous for parameters like material prop-
erties. Its skewness accommodates asymmetry in probability distributions, aligning well with uncertainties in structural and earthquake engineering,
where factors such as material strength or damping may exhibit skewed behavior. Additionally, the logarithmic transformation simplifies multiplicative
processes commonly encountered in material properties and structural responses.

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1174 HARIRI-ARDEBILI

F IGURE 4 Downstream view and cross-section geometry of monolith No. 16 of Pine Flat Dam; adapted from [72]. Both the 2D and 3D
sliced finite element models of this cross-section are used for modeling UQ.

uncertainties) significantly amplifies system uncertainty. Low-fidelity models typically involve fewer unknown variables
compared to high-fidelity models. For instance, consider the dam-foundation joint interface: a low-fidelity model might
adopt the Mohr-Coulomb constitutive model, which comprises only two parameters (cohesion and friction angle). In
contrast, a high-fidelity model might utilize a zero thickness, fracture mechanics inspired cohesive crack model with up
to 14 parameters,64 some of which are either practically challenging or very expensive to obtain through experimental
tests. The increase in parameters in high-fidelity models, especially those associated with constitutive model parameters,
contributes significantly to amplifying system uncertainty.
It is worth noting that these approaches might exhibit considerable bias when compared to the actual seismic

performance of the system (that necessitates model validation and verification).

4 PROBLEM STATEMENT AND BENCHMARK DESCRIPTION

Considering the preceding discourse, which emphasizes that the most significant modeling uncertainty lies in the uncer-
tainty associated with the model methodology itself, the omission of such modeling uncertainty raises a critical question:
“Does the utilization of modeling technique, which overlook the majority of the fundamental physics underlying the
problem, align with the principles of rigorous uncertainty analysis?”
To shed light on this query, this paper undertakes, for the first time, an examination of the outcomes from a meticu-

lously structured benchmarkworkshop. Theworkshop design deliberately isolates the levels 3 and 4modeling uncertainty
sources by providing participantswith comprehensive data. The benchmark problemencompasses an intricate damgeom-
etry alongside its corresponding finite element mesh. Essential material mechanical properties for the mass concrete and
foundation rock, including concrete strength and fracture parameters, as well as shear and compression wave velocities
within the rock, are supplied. Notably, variations in reservoir water level corresponding to winter, summer, and normal
pool conditions are detailed. The document explicitly outlines the requisite loading scenarios for the finite element simula-
tion, encompassing seismic loads. Both seismic excitation at the free surface and its deconvolved version at the foundation
base are furnished. Additionally, the benchmark includes recommendations for damping values and formulation.
By providing this extensive dataset to participants, the benchmark effectively eliminates ground motion RTR variabil-

ity, material randomness, and levels 1 and 2 modeling uncertainties. Participants are granted the freedom to select any
finite element software for analysis, explore diverse approaches for modeling the foundation and reservoir sections and
their interactions with the dam, consider wave propagation complexities, implement various boundary conditions, opt
for different numerical solvers (such as implicit or explicit methods), and decide on concrete constitutive models and
damage behaviors, among other factors. These considerations collectively fall under the umbrella of levels 3 and 4
modeling uncertainties.
The case study is the tallest non-overflow monolith (#16) of Pine Flat Dam. Constructed by the US Army Corps of

Engineers in 1954, this dam comprises 36 monoliths, of which 35 are 15.25 m-wide and 1 measures 12.2 m-wide. The dam
extends linearly for a distance of 561 m and boasts a maximum height of 122 m in the case of the tallest non-overflow
monolith, as depicted in Figure 4. Pine Flat Dam has been the subject of extensive study during the 1970s and 1980s at the
University of California at Berkeley.65–71 The reference coordinate system is situated at the dam heel and the mid-width
of the monolith. The y-axis aligns with the downstream direction, while the z-axis points upward. It is assumed that the

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1175

TABLE 1 Materials properties in the finite element model of Pine Flat Dam72.

Parameter Symbol Unit Quantity
Concrete Modulus of elasticity 𝐸𝑐 MPa 22,410

Poisson’s ratio 𝜈𝑐 − 0.2
Mass density 𝜌𝑐 Kg/m3 2,483
Compressive strength 𝑓𝑐 MPa 28.0
Tensile strength 𝑓𝑡 MPa 2.0
Fracture energy 𝐺𝑓 N/m 250
Compressive strain at peak load 𝜖𝑐 − 0.0025
Tensile strain at peak load 𝜖𝑡 − 0.00012

Rock Modulus of elasticity 𝐸𝑓 MPa 22,410
Poisson’s ratio 𝜈𝑓 − 0.33
Mass density 𝜌𝑓 Kg/m3 2,483
Shear wave velocity 𝑉𝑆 m/s 1,939
Compressional wave velocity 𝑉𝑃 m/s 3,167

Water Mass density 𝜌𝑤 Kg/m3 1,000
compression wave velocity 𝐶𝑤 m/s 1,439
Winter water height 𝐿𝑊 m El. 268.21
Summer water height 𝐿𝑆 m El. 278.57
Normal water height 𝐿𝑁 m El. 290.00

dam-foundation interface is fully connected, thereby precluding any potential sliding. This assumption aligns with the
benchmark’s aim to streamline damage estimation using the continuum crack model. This deliberate choice enhances
consistency and simplifies further UQ.
In the context of the nonlinear simulations discussed within this paper, the foundation block beneath the dam spans

dimensions of 700mby 122m. Positioned 305m from the left edge of the foundation block, the damheel is located.Material
properties, as outlined in Table 1, are presumed homogeneous and isotropic for both the concrete and rock components.
Static load considerations encompass solely the weight of concrete and water. Participants were directed to incorporate
viscous damping based on their prevailing comprehension, modeling conventions, and chosen analytical methods, with
an expectation for thorough documentation.
The finite elementmodel representing the dam–water–foundation rock system incorporates three primary contributors

to energy dissipation: (1) material damping inherent in the concrete and rock; (2) radiation damping occurring in the
semi-unbounded foundation and fluid domains; and (3) dissipation of hydrodynamic wave energy at the reservoir
boundaries.73 The modeling of material damping in both concrete and rock introduces significant uncertainty. The lack
of experimental data on the energy-dissipating characteristics of concrete and rock under seismic loading conditions
renders phenomenological modeling impractical.74 Consequently, simplified viscous damping models are employed for
their mathematical tractability. Radiation damping emerges as a predominant source of energy dissipation in 2D dam
models, unless the foundation rock significantly exceeds the stiffness of the concrete.75 However, for 3D dam models,
this effect diminishes in significance.
In ICOLD workshop, for nonlinear simulations, participants were encouraged to explore the utilization of Rayleigh

viscous damping, with a recommended value of 2% applied to both the dam and foundation. Based on a compilation of
damping data from 32 concrete dams, obtained through forced vibration field tests and estimated from ambient vibra-
tion measurements,73 demonstrated that damping levels range between 1% and 5%, spanning from the 10th to the 90th
percentiles. Specifically, for Pine Flat Dam, damping was measured in the range of 2% to 3.5%.
Two distinct types of loads were furnished for the purpose of nonlinear transient simulations: a real ground motion

record and an artificially generated acceleration profile. The real acceleration data employed originates from the M 7.3
Kern County, California earthquake on July 21, 1952, recorded at the Lincoln School in Taft. Specifically, only the S69E
component of this ground motion was considered, with a peak ground acceleration (PGA) value of 0.18 g. The selec-
tion of the Taft ground motion, originally utilized in frequency-domain linear elastic analysis by Chopra et al.66–68 for
Pine Flat Dam, serves the purpose of facilitating model calibration. By allowing participants to compare their simulation

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1176 HARIRI-ARDEBILI

0 10 20 30 40 50 60
-0.2

0

0.2

0 10 20 30 40 50 60
-0.1

0

0.1

0 10 20 30 40 50 60
-0.1

0

0.1

0 1 2 3 4
0

0.5

1

(A) Real ground motion record

0 5 10 15
-1

0

1

0 5 10 15
-0.5

0

0.5

0 5 10 15

-0.4

-0.2

0

0 1 2 3 4
0

5

(B) Intensifying artifical acceleration

F IGURE 5 Acceleration, velocity, and displacement time histories, along with the response spectrum of the input excitation. Note that
the response spectrum is presented for the entire duration of real ground motion, whereas it is time-dependent in the case of IAA.

results, particularly in the linear elastic domain, with previous studies, the benchmark aims to mitigate potential bias in
probabilistic simulations arising from substantial errors in one participant’s predictions.
Given the assumption that this signal represents free surface conditions, its deconvolved counterpart is also included.

The deconvolution procedure assumes a vertically-propagating SH-wave within a uniform half-space characterized by the
samematerial properties as the foundation. Importantly, in the context of a free surface signal, themotion is assumed to be
recorded at the top of the foundation in the absence of the dam and reservoir. Consequently, the free-surface motions are
suitable for input inmassless foundationmodels, whereas the deconvolved signals are intended for employment inmassed
foundation models. The significance of deconvolving seismic waves using both frequency- and time-domain approaches,
along with a detailed implementation in the analysis of dam-foundation coupled systems, is extensively discussed in [76,
77]. Figure 5A illustrates the acceleration, velocity, displacement, and force time histories of the Taft ground motion,
accompanied by the 2% acceleration response spectrum.
The second signal constitutes an intensifying artificial acceleration (IAA) function. Such artificial accelerations are gen-

erated through unconstrained optimization techniques78 with a near-linear increase in peak acceleration and response
spectrum over time. This methodology can be likened to a dynamic pushover procedure, which assesses the seismic per-
formance of a dam across the spectrum from linear elastic behavior to nonlinear response, ultimately capturing theoretical
collapse. Figure 5B provides an illustration of the provided IAA to participants, showcasing its associated velocity, displace-
ment, and force time histories. Furthermore, acceleration response spectra are depicted at three distinct time points (5,
10, and 15 s). For further insights into the concept and application of IAA in dam engineering, refer to [79–83]. Validation
of IAA with respect to classical IDA has been discussed in [84, 85].

5 ANALYSIS AND SYNTHESIS OF RESULTS

This section presents a concise overview of the outcomes obtained from both real ground motion records and artificial
accelerations. The following scalar and vector quantities were analyzed: crest and heel displacements, relative crest dis-
placement compared to heel, net hydrodynamic pressure at the dam’s heel, ratio of crest-to-heel amplitude acceleration
spectrum, crack profile, area-based and length-based damage index (DI), and failure time (for IAA only). The data ana-
lyzed in this paper is derived from the findings presented in the following publications.86–104 Table 2 provides an overview
of key modeling variations observed across each study. These modeling variabilities encompass factors such as mesh type
(2D, 3D slice, full 3D), mesh size (which holds significance for DI assessment), time step utilized in the numerical simula-
tion, criteria for convergence, the specific foundation model and its associated boundary conditions for wave propagation

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1177

T
A
B
L
E

2
Su
m
m
ar
y
of
th
e
m
od
el
in
g
ch
oi
ce
sa
nd

va
ria
bi
lit
ie
sb
y
di
ffe
re
nt
pa
rt
ic
ip
an
ts
.

ID
M
es
h
ty
pe

M
es
h
si
ze

[m
]

Ti
m
e
st
ep

C
on
ve
rg
en
ce

cr
it
er
ia

Fo
un

da
ti
on

ty
pe

an
d
bo
un

da
ry

co
nd

it
io
ns

Fl
ui
d-
st
ru
ct
ur
e
in
te
ra
ct
io
n

C
on
cr
et
e
co
ns
ti
tu
ti
ve

(d
am

ag
e)

m
od
el

Ti
m
e
in
te
gr
at
io
n

So
ft
w
ar
e

11
2D

pl
an
e
st
ra
in

2
N
or
m
of
un
ba
la
nc
e
fo
rc
e

D
om

ai
n
re
du
ct
io
n
m
et
ho
d
w
ith

ab
so
rb
in
g
da
m
pi
ng

la
ye
rs
;N

o

fr
ee
-fi
el
d
BC

Fi
ni
te
vo
lu
m
e
m
es
h
fo
rf
lu
id
;

co
m
pr
es
si
bl
e;
ac
ou
st
ic
-s
tr
uc
tu
ra
l

co
up
lin
g

Im
pl
ic
it
(s
ol
id
);
ex
pl
ic
it

(fl
ui
d
co
up
lin
g)

Re
al
ES
SI

12
2D

pl
an
e
st
ra
in

2.
5–
5.
0

0.
01

1%
Pe
rf
ec
tly

m
at
ch
ed

la
ye
rs
;F
re
e-
fie
ld
BC

on
ly
on

th
e
bo
tto
m
of
fo
un
da
tio
n

In
co
m
pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

To
ta
ls
tr
ai
n-
ba
se
d
cr
ac
k
m
od
el
w
ith

fix
ed
-c
ra
ck

or
ie
nt
at
io
n
in
cl
ud
in
g

lin
ea
rt
en
si
le
cu
rv
e

Ex
pl
ic
it

D
ia
na

13
2D

pl
an
e
st
ra
in

10
0.
00
1

0.
00
01

V
is
co
us
-s
pr
in
g
bo
un
da
rie
sw

ith
m
as
se
d

fo
un
da
tio
n;
W
ith

fr
ee
-fi
el
d
BC

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

A
ni
so
tr
op
ic
da
m
ag
e
m
od
el
fo
r

pr
e-
pe
ak

an
d
or
th
ot
ro
pi
c
da
m
ag
e

fo
llo
w
in
g
Ra
nk
in
e
cr
ite
ria

fo
r

po
st
-p
ea
k;
Is
ot
ro
pi
c
da
m
ag
e
un
de
r

co
m
pr
es
si
on

st
re
ss
es
us
in
g

D
ru
ck
er
-P
ra
ge
rc
rit
er
ia

Im
pl
ic
it

C
od
e-
A
st
er

14
3D

sl
ic
e

2
0.
01

N
ew

to
n
Ra
ph
so
n
(F
or
ce
)

Fo
un
da
tio
n
bo
tto
m
is
ki
ne
m
at
ic
al
ly

co
ns
tr
ai
ne
d
(I
nf
in
ite

ra
di
at
io
n

bo
un
da
ry
);
N
o
fr
ee
-fi
el
d
BC

C
om

pr
es
si
bl
e
flu
id
;

au
co
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

D
ru
ck
er
-P
ra
ge
rc
on
st
itu
tiv
e
m
od
el

w
ith

pe
rf
ec
tp
la
st
ic
ity

Ex
pl
ic
it
(N
ew

m
ar
k)

A
ns
ys

15
2D

pl
an
e
st
ra
in

1
0.
01

N
ew

m
ar
k

Si
m
pl
e
su
pp
or
ts
;N

o
fr
ee
-fi
el
d
BC

In
co
m
pr
es
si
bl
e
flu
id
;F
SI
w
ith

lin
ks
-g
ap

m
et
ho
d

C
on
cr
et
e
m
at
er
ia
li
n
SA
P
w
ith

no
nl
in
ea
rp
ro
pe
rt
ie
s

Ex
pl
ic
it
(N
ew

m
ar
k
an
d

Li
ne
ar
)

Sa
p2
00
0

16
3D

sl
ic
e

1.5
0.
00
5

N
od
al
un
ba
la
nc
ed

fo
rc
es

V
is
co
us
-s
pr
in
g
ar
tif
ic
ia
lb
ou
nd
ar
ie
s

m
od
el
w
ith

m
as
se
d
fo
un
da
tio
n;
W
ith

fr
ee
-fi
el
d
BC

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
on
cr
et
e
da
m
ag
e
pl
as
tic
ity

Im
pl
ic
it
(H
H
T)

A
ba
qu
s

17
2D

pl
an
e
st
ra
in

3
0.
00
00
5

N
on
-r
ef
le
ct
in
g
vi
sc
ou
sb
ou
nd
ar
y
w
ith

fr
ee
-fi
el
d
co
nd
iti
on
sa
tt
he

la
te
ra
l

bo
un
da
rie
s

La
gr
an
gi
an

di
sp
la
ce
m
en
t-b
as
ed

flu
id

w
ith

ro
ta
tio
na
lp
en
al
ty
co
ns
tr
ai
nt
;

FS
I

D
is
cr
et
e
cr
ac
k
m
od
el
w
ith

a
bi
-li
ne
ar

so
fte
ni
ng

m
od
el
in
th
e
no
rm

al
an
d

sh
ea
rd
ire
ct
io
ns
an
d
Ro
ku
go

cu
rv
e

Ex
pl
ic
it
(c
en
tr
al

di
ffe
re
nc
e)

D
A
M
FA

2D

18
2D

pl
an
e
st
ra
in

3
0.
00
2

H
ilb
er
-H
ug
he
s-T
ay
lo
r(

𝛼
=

-0
.3
)

N
on
-r
ef
le
ct
in
g
bo
un
da
ry
co
nd
iti
on
s

us
in
g
pa
rt
ic
ul
ar
in
te
rf
ac
e
el
em

en
t

ch
ar
ac
te
riz
ed

by
a
ve
ry
lo
w
st
iff
ne
ss

an
d
da
m
pi
ng

co
ef
fic
ie
nt
;N

o

fr
ee
-fi
el
d
BC

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
ra
ck

w
id
th
by

m
ul
tip
ly
in
g
el
em

en
ta
l

cr
ac
k
st
ra
in
an
d
th
e
cr
ac
k
ba
nd
w
id
th

(s
et
to
0.
5
m
)

Im
pl
ic
it
(H
H
T)

D
ia
na

19
2D

pl
an
e
st
ra
in

15
0.
00
1

Fi
rs
tv
ib
ra
tio
n
fr
eq
ue
nc
y

Li
ne
ar
pa
ra
xi
al
ab
so
rb
in
g
bo
un
da
ry

el
em

en
ts
;W

ith
fr
ee
-fi
el
d
BC

Bo
th
ac
ou
st
ic
flu
id
an
d
W
es
te
rg
aa
rd
’s

ad
de
d
m
as
s;
FS
I

Eq
ui
va
le
nt
lin
ea
riz
at
io
n
te
ch
ni
qu
e

(te
ns
ile
-s
tr
ai
n
co
nt
ro
lle
d)

Im
pl
ic
it
(N
ew

m
ar
k)

C
od
e-
A
st
er

20
2D

pl
an
e
st
ra
in

1.5
0.
01

Pa
ck
ed

sc
he
m
e:
m
od
er
at
e

di
ss
ip
at
io
n

In
fin
ite

el
em

en
t;
W
ith

fr
ee
-fi
el
d
BC

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
on
cr
et
e
da
m
ag
ed

pl
as
tic
ity

Im
pl
ic
it
(D
ire
ct
)

A
ba
qu
s

(C
on
tin
ue
s)

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1178 HARIRI-ARDEBILI

T
A
B
L
E

2
(C
on
tin
ue
d)

ID
M
es
h
ty
pe

M
es
h
si
ze

[m
]

Ti
m
e
st
ep

C
on
ve
rg
en
ce

cr
it
er
ia

Fo
un

da
ti
on

ty
pe

an
d
bo
un

da
ry

co
nd

it
io
ns

Fl
ui
d-
st
ru
ct
ur
e
in
te
ra
ct
io
n

C
on
cr
et
e
co
ns
ti
tu
ti
ve

(d
am

ag
e)

m
od
el

Ti
m
e
in
te
gr
at
io
n

So
ft
w
ar
e

21
3D

sl
ic
e

1.0
-3
7.
0

0.
01

Re
la
tiv
e
co
nv
er
ge
nc
e

cr
ite
rio
n:
10
e-
6

Fi
xe
d
BC

fo
rm

as
sl
es
sm

od
el
;N

o

ab
so
rb
in
g
BC

A
dd
ed

m
as
sa
nd

ac
ou
st
ic
flu
id

N
on
lin
ea
rt
en
si
le
be
ha
vi
or
on
ly

Im
pl
ic
it
(N
ew

to
n)

C
od
e-
A
st
er

23
2D

pl
an
e
st
ra
in

0.
5–
4.
0

0.
00
5

5%
fo
rf
or
ce
an
d
m
om

en
t,

0.
01
%
fo
rm

as
sf
lo
w
ra
te
,5
%

fo
rd
is
pl
ac
em

en
t

V
is
co
us
-s
pr
in
g
bo
un
da
ry
co
nd
iti
on
;N

o

fr
ee
-fi
el
d
BC

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
ou
pl
ed

da
m
ag
e-
pl
as
tic
ity

m
od
el

ba
se
d
on

m
ic
ro
-p
la
ne

fo
rm

ul
at
io
n

Im
pl
ic
it

A
ns
ys

24
3D

sl
ic
e

1.4
0.
00
00
05
67

V
is
co
us
bo
un
da
rie
s(
Ly
sm

er
);

Fr
ee
-fi
el
d
el
em

en
ts
fo
rl
at
er
al

bo
un
da
rie
so
nl
y

C
om

pr
es
si
bl
e
flu
id
w
ith

m
ix
ed

di
sc
re
tiz
at
io
n
sc
he
m
e

El
as
tic

pe
rf
ec
tly

pl
as
tic

co
ns
tit
ut
iv
e
la
w

fo
rs
he
ar
w
ith

M
oh
r-
C
ou
lo
m
b

fa
ilu
re
cr
ite
rio
n;
el
as
tic

pe
rf
ec
tly

fr
ag
ile

co
ns
tit
ut
iv
e
la
w
fo
rt
en
si
le

Ex
pl
ic
it

Fl
ac
3D

26
2D

pl
an
e
st
ra
in

1
0.
01

En
er
gy
ba
se
d

In
fin
ite

el
em

en
ts
w
ith

fr
ee
-fi
el
d
BC

ba
se
d
on

Ly
sm

er

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
on
cr
et
e
da
m
ag
e
pl
as
tic
ity

Im
pl
ic
it

A
ba
qu
s

27
3D

sl
ic
e

3
0.
00
01

V
is
co
us
bo
un
da
ry
co
nd
iti
on
sp
ro
po
se
d

by
Ly
sm

er
;W

ith
fr
ee
-fi
el
d
BC

La
gr
an
gi
an

di
sp
la
ce
m
en
t-b
as
ed

flu
id

w
ith

on
ly
no
rm

al
st
iff
ne
ss
;F
SI

Bi
-li
ne
ar
so
fte
ni
ng

m
od
el
in
te
ns
io
n

an
d
sh
ea
rw

ith
Ro
ku
go

cu
rv
e

Ex
pl
ic
it
(C
en
tr
al

di
ffe
re
nc
e)

Pa
rm

ac
3D

28
2D

pl
an
e
st
ra
in

1.5
0.
01

Ly
sm

er
bo
un
da
ry
co
nd
iti
on
;W

ith

fr
ee
-fi
el
d
BC

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
on
cr
et
e
da
m
ag
e
pl
as
tic
ity

Im
pl
ic
it

A
ba
qu
s

29
2D

pl
an
e
st
ra
in

2.
5

0.
01

V
is
co
us
sp
rin

g
el
em

en
ts
;N

o
fr
ee
-fi
el
d

BC

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
on
cr
et
e
pl
as
tic
ity

w
ith

M
en
et
re
y-
W
ill
am

cr
ite
rio
n

Im
pl
ic
it

A
ns
ys

30
2D

pl
an
e
st
ra
in

3.
5

In
fin
ite

el
em

en
ts
;f
re
e
fie
ld

ro
ck
-c
ol
um

n
an
d
di
re
ct
fr
ee
fie
ld

fo
rc
es

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
on
cr
et
e
da
m
ag
e
pl
as
tic
ity

Im
pl
ic
it

A
ba
qu
s

31
2D

pl
an
e
st
ra
in

1.5
0.
00
1

V
is
co
el
as
tic

ar
tif
ic
ia
lb
ou
nd
ar
y;
N
o

fr
ee
-fi
el
d
BC

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
on
cr
et
e
da
m
ag
e
pl
as
tic
ity

Im
pl
ic
it

A
ba
qu
s

32
3d

sl
ic
e

10
0.
01

To
le
ra
nc
e
fo
rl
ar
ge
st
no
da
l

va
lu
e
av
ai
la
bl
e
in
th
e

sy
st
em

+
en
er
gy
no
rm

N
on
-r
ef
le
ct
in
g
BC

sw
ith

da
m
pe
r

el
em

en
ts

C
om

pr
es
si
bl
e
flu
id
-li
ke

st
ru
ct
ur
al

el
em

en
ts
;F
SI

C
on
cr
et
e
m
at
er
ia
la
cc
or
di
ng

to
EN

no
rm

s

W
ils
on

th
et
a

SO
Fi
ST
iK

33
Pl
an
e
St
ra
in

3
6e
-4
to
0.
01

Fo
rc
e
eq
ui
lib
riu

m
In
fin
ite

el
em

en
ts
;N

o
fr
ee
-fi
el
d
BC

C
om

pr
es
si
bl
e
flu
id
;

ac
ou
st
ic
-s
tr
uc
tu
re
in
te
ra
ct
io
n

C
on
cr
et
e
da
m
ag
e
pl
as
tic
ity

Im
pl
ic
it

A
ba
qu
s

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1179

(A) Detailed displacement time history at different locations and the median response.

0 20 40
0

20

40

60

0 20 40
0

20

40

60

0 20 40
0

10

20

30

(B) Standard deviation of displacement response at different locations.

0 5 10 15
0

0.1

0.2

0.3
Data
GEV
Lognormal
Gamma

0 50 100
0

0.02

0.04

0.06
Data
GEV
Lognormal
Gamma

0 5 10
0

20

40

60

80

100

(C) Distribution of different diaplcement quantities and their relationship.

F IGURE 6 Processing displacement response under real ground motion record.

and absorption (including considerations for the free-field), choice of fluid element and approach for FSI, the selected
concrete constitutive model and its response to dynamic loading-induced damage, the time integration scheme employed
for stability control, and finally, the software package employed for analysis.

5.1 Real ground motion

Figure 6 presents a summary of the results collected from 20 models, illustrating the displacement response. Prior to
statistical analysis, all results underwent initial data filtering, and any outliers were subsequently removed. Individual
displacement time histories for both the heel and crest, alongwith the relative crest displacement, are shown in Figure 6A.
This figure also highlights themedian values. The heel and crest displacements exhibit large oscillations, while the relative
displacement demonstrates a high-frequency nature. Fourmodels exhibit residual displacements at both the heel and crest
points, which may indicate either relative sliding between the dam and the foundation block or a rigid body motion of
the entire finite element model due to inadequate boundary conditions. This issue can often be resolved by introducing

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1180 HARIRI-ARDEBILI

a self-equilibrium loading condition, such as elastic boundary conditions, to one side of the foundation.1 Among these
four models, three employ viscous boundary conditions for the foundation, with one of them utilizing infinite elements.
Additionally, one simulation exhibits a temporal shift and employs infinite elements for the foundation, including the
rock-column and direct free field forces. In terms of the relative crest displacement history, two models exhibit residual
displacement, both of which involve infinite elements. Notably, the most substantial relative crest displacement arises in
a model utilizing a massless foundation approach with fixed boundary conditions. This case yields a maximum relative
crest displacement approximately 2.5 times higher than the median value. This observation aligns with prior findings that
distinguish between massed and massless foundation models.1
Figure 6B illustrates the standard deviation of displacement histories for both absolute and relative measurements.

The absolute displacements at both the heel and crest locations exhibit a nearly linear increasing pattern, punctuated by
intermittent spikes along this trajectory. Within the strong motion interval of the analysis (5–15 s), these displacements
demonstrate standard deviations ranging from 10 to 30 mm. The highest standard deviation is approximately 50 mm,
which can be compared to the previously mentioned median displacement of around 100 mm. In contrast, the pattern for
the standard deviation of the relative crest displacement differs. Here, the maximum standard deviation occurs around 10
s (coinciding with the peak ground motion), followed by a decrease to a smaller and relatively stable value for the rest of
the simulation.
All the previous comparisons have been conducted exclusively on “dynamic only” results, excluding the static response.

Figure 6C provides a comprehensive overview of the maximum static and maximum dynamic displacement responses
separately, along with three distributional models fitted to the data. The choice of using the generalized extreme value
(GEV), lognormal, and gamma3 distributions is arbitrary, as the goal here is not to determine themost suitable fit based on
goodness-of-fitmetrics. Given the limited scope of the database used in this study, any conclusion drawn from thesemodels
should not be regarded as definitive. Aswewill explore later in this section,with the exception of static displacement, these
three models yield similar outcomes for all other cases. For maximum static displacement, the GEV model appears more
reasonable, since the other two models result in very small mean or median values. Additionally, it is worth noting that
the lognormal distribution is commonly employed in the fields of structural and earthquake engineering to describe the
distribution of engineering demand parameters.105
For the maximum static displacement, the values vary from 0.2 mm to over 12 mm. The GEV model has a location

value of 3.4 mm, which is close to the crude mean of 4.6 mm. Concerning the maximum dynamic displacement, values
range from 13 to 112 mm, with a crude mean of 43 mm. All three distribution models exhibit striking similarity. The GEV
model’s𝜇 (location) value is 33mm,while the lognormalmodel’smedian is 39mm, and the logarithmic standard deviation
is 0.45. This value coincides with the recommendations from EMA P695106 regarding modeling uncertainty. EMA P695
characterizes modeling uncertainty through a matrix with two dimensions: (1) Representation of collapse characteristics
and (2) Accuracy and robustness of models4. Each dimension encompasses three levels: high, medium, and low. For
scenarios with a low representation of collapse and medium model accuracy (a condition that can be argued to apply in
this study), a dispersion of 0.5 is recommended, which is analogous to the logarithmic standard deviation in the lognormal
distribution. It is important to note that while this comparison is approximate, as FEMAP695 pertains directly to buildings
and not concrete dams, there are similarities between them, as a concrete gravity dam functions as a cantilever-type
short-period multi-degree-of-freedom system. They both share dynamic traits, responding to seismic forces through the
transmission of loads at their bases. Additionally, their seismic behavior is influenced by the first few vibration modes.
Lastly, Figure 6C (left) depicts the correlation between the maximum static and dynamic displacements. As seen, there

is minimal correlation observed across different models. The static displacement is primarily influenced by hydrostatic
pressure and the dam’s self-weight, along with the stiffness of the foundation. It is worth mentioning that in certain
instances, the hydrostatic pressure is applied solely to the dam, while in others, it is also applied to the upper face of
the foundation.
Figure 7 employs a similar approach to the acceleration response analysis as discussed for displacement in Figure 6.

Acceleration at the heel and crest points is depicted in Figure 7A (left and middle). Ideally, the acceleration at the heel
should be similar or close to the free-surface Taft record displayed in Figure 5A. However, the median of the recorded heel

3 The GEVmodel includes three parameters as shape, scale, and location factors. The lognormal model includes two parameters as log location and log
scale. The gamma model includes two parameters as shape and scale factors.
4 Representation of collapse characteristics refers to how completely and comprehensively the index archetype models capture the full range of design
parameters and associated structural collapse behavior that is envisioned within the archetype design space. Accuracy and robustness of models refers
to the degree to which nonlinear behaviors are directly simulated in the model, or otherwise accounted for in the assessment.

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1181

(A) Detailed acceleration time history at different location and the median response.

0 20 40
0

0.05

0.1

0 20 40
0

0.2

0.4

0.6

10-1 100 101 102
10-2

100

102

(B) Standard deviation of acceleration response at different locations.

0.1 0.2 0.3 0.4
0

2

4

6

8

10 Data
GEV
Lognormal
Gamma

0.5 1 1.5 2 2.5 3
0

0.5

1

1.5

2

2.5
Data
GEV
Lognormal
Gamma

0.1 0.2 0.3 0.4

1

1.5

2

2.5

(C) Distribution of different acceleration quantities and their relationship.

F IGURE 7 Processing acceleration response under real ground motion record.

acceleration exhibits a smaller value. This discrepancy can be attributed to themodeling uncertainty in the foundation and
its boundary conditions, including the free-field assumptions. As anticipated, the acceleration at the crest is amplified by
approximately 10 times. To better represent the heel and crest acceleration time histories, they are presented as amplitude
spectra in Figure 7A (right), showing the ratio of the acceleration amplitude spectrum. This ratio is obtained by computing
the Fourier spectrum of the acceleration time series.
For frequencies below 1 Hz, the median response ratio is around 2. Between the frequency range of 5 to 10 Hz, the ratio

fluctuates between 3 and 30, with an average value of 10. Beyond 10 Hz (up to about 50 Hz), the ratio gradually decays
linearly to 0.5. Among the individual models, three show significant spikes, such as 170, in the amplitude spectrum ratio.
Notably, all three models lack free-field boundary conditions.
The standard deviation of heel and crest acceleration time series follows a similar trend as displacement, peaking at the

heel with a value of 0.1 g, equal to the median acceleration, Figure 7B. The standard deviation of acceleration spectrum
amplitude averages around 0.8 for frequencies below 2 Hz, and then increases to an average value of 5 for 𝑓 > 2 Hz,
accompanied by fluctuations.

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1182 HARIRI-ARDEBILI

F IGURE 8 Processing hydrodynamic pressure response under real ground motion record.

0 0.05 0.1 0.15 0.2
0

2

4

6

8

10
Data
GEV
Gamma

5 10 15
0

5

10

15

F IGURE 9 Processing DI under real ground motion record.

The fitted distributions for maximum (absolute) heel and crest accelerations are shown in Figure 7C (left and middle).
The lognormal and gamma models show close agreement, while the GEV model exhibits a minor difference. Based on
the lognormal distribution, the median values for maximum heel and crest accelerations are 0.2 and 1.07 g, respectively.
The logarithmic standard deviation (dispersion) in both cases is approximately 0.3, which is lower than the displacement-
based dispersion. Finally, Figure 7C (right) illustrates the relationship between maximum acceleration at the heel and
crest. With the exception of two outliers, a positive correlation is evident. Further analysis of these two outliers indicates
that both suffer from spikes in their heel acceleration responses, which can be attributed to numerical integration issues
or the presence of superstructure (i.e., the dam) at the recording location.
Figure 8 employs a similar approach to analyze the hydrodynamic pressure, 𝑃ℎ𝑦𝑑, at the heel of the dam. Unlike the

displacement response, no instances of residual or shifted values are observed for this parameter. The ratio of themaximum
𝑃ℎ𝑦𝑑 to themedian value is approximately 2. The trend in standard deviation is comparable to that observed for the relative
crest displacement, with a peak value of around 0.15MPa. It is noteworthy that themaximum 𝑃ℎ𝑦𝑑 in themedian response
is even lower than 0.15 MPa. Lastly, Figure 7 (right) presents the distribution of the maximum 𝑃ℎ𝑦𝑑, which is consistent
across the three distribution models. The lognormal model yields a median value of 0.16 MPa and a logarithmic standard
deviation of 0.3, indicating less variation compared to the displacement response.
Figure 9 provides the DI for different simulations. Damage at the dam base (i.e., dam-foundation interface) is quantified

as the ratio of cracked base length to the total base length,78 𝐷𝐼 =
𝐿𝑐𝑟

𝐿𝑡𝑜𝑡

. As seen in this figure, some simulations do not report
any damage while some others report a damage as large as 18%. The median value from all simulations is about 4.5%. Five
out of 16 simulations report DI of 10% of higher. Four of these models (including the two top ones) are all based on no
free-field assumption for the foundation domain. Also, four of these models are associated with damage plasticity, and
four simulations use implicit time integration scheme. Tracking the tends of models with no or small damage also reveals
that they are different types of time integration scheme and concrete constitutive models; however, nearly all of them
are based on free-field boundary condition. Therefore, the most influential modeling aspect in DI variation is foundation
boundary condition and wave propagation assumption.

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1183

F IGURE 10 Processing displacement-based responses under IAA.

5.2 Intensifying Artificial Acceleration

While conventional modeling UQ using a single real ground motion record is inherently limited by the specific attributes
of that record, such as its frequency content and seismic intensity, the IAA offers a distinctive avenue for probing a broad
spectrum of intensity levels without incurring additional computational costs. IAA further facilitates the comprehensive
exploration of dam responses across varying seismic intensities in a continuous manner. In this context, 16 participants
undertook analyses employing IAA. Similar response quantities as of the real ground motion were recorded.
Figure 10 (left) illustrates the temporal evolution of the maximum absolute relative crest displacement. In this plot, the

displacement response terminated at the failure point. The determination of this failure time hinges on two primary crite-
ria: (1) when the simulation reaches a juncture where certain engineering demand parameters surpass defined thresholds,
for instance, when theDI at the base attains a value of one. Some participants used a displacement-based criteria, for exam-
ple, 0.01% of dam height as failure criteria (other references propose 0.03%107). (2) when numerical instability emerges,
triggering the automatic termination of the simulation due to unsatisfactory software internal criteria. Although the pre-
cision of establishing a failure point introduces an inherent uncertainty, it is important to acknowledge the absence of
a universally accepted metric to quantitatively define failure (or collapse) in concrete dams. Notably, the U.S. Bureau of
Reclamation’s definition of structural failure as the uncontrolled release of reservoir water presents a qualitative descrip-
tion unamenable to integration within UQ efforts. As illustrated in Figure 10 (left), failure is observed to manifest at
diverse levels of crest displacement and varied time intervals, ranging from 2 to 15 s, accompanied by failure displacements
spanning 30 to 600 mm, with one notable instance reaching 1200 mm.
In IAA, the time parameter exhibits a direct correlation with seismic intensity. This relationship is explained by [78,

79], demonstrating the potential to convert the “time” parameter into different intensity measures, such as PGA or first-
mode spectral acceleration. This study adopts the simplest approach: PGA is represented as PGA(𝑡) ≡ max{Abs(�̈�𝑔(𝜏) ∶

𝜏 ∈ [0, 𝑡])}. Unlike the scalar nature of PGA in real ground motions, it is a time-dependent parameter within IAA appli-
cation. Subsequently, the definition of displacement-based capacity curves ensues, depicting the relationship between
displacement and PGA, as shown in Figure 10 (middle). Increasing PGA yields a nonlinear displacement response, ulti-
mately leading to failure, as evidenced by the plateau in the capacity curve. Drawing inspiration from the IDAmethod,108
themedian, 16% and 84% fractiles are calculated for each individual capacity curve. The failure points are delineated by the
right vertical line, allowing for the tracking of the failure PGA value, PGA𝑓𝑎𝑖𝑙𝑢𝑟𝑒. It is paramount to reiterate that the vari-
ations observed in the individual capacity curves within Figure 10 (middle) are solely attributed to modeling uncertainty,
devoid of any reflection of ground motion RTR variability, as dispelled by [109].
Leveraging the data obtained from failure capacity points, it is possible to construct an empirical CDF and fit various

distributional models. In this study, it is assumed that these data points conform to a lognormal distribution, resulting in
the development of a failure fragility-type function. Two distinct methodologies are employed: (1) Method of Moments
(MM), which aligns the median and logarithmic standard deviation to match the moments of the data, and (2) Maximum
Likelihood Estimation (MLE), which optimizes parameters to maximize a “likelihood” function based on either failed
or safe data.110 The comparison of these failure fragility functions is presented in Figure 10 (right). Notably, the MM
and MLE methods exhibit remarkable closeness, yielding estimated medians of 0.39 and 0.37, and dispersions of 0.67
and 0.76, respectively. These dispersions surpass the RTR variability reported for the same concrete gravity dam, that is,
0.4.34

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1184 HARIRI-ARDEBILI

F IGURE 11 Processing damage indices under IAA.

Figure 11 illustrates capacity curves, fragility functions, and the distribution of the DI at the failure time, employing two
different versions of theDI. The first one, the length-basedDI discussed earlier, is complemented by the area-basedDI. The
area-based DI is quantified as the ratio of the cracked (cross-section) area to the total dam area, defined as 𝐷𝐼 =

𝐴𝑐𝑟

𝐴𝑡𝑜𝑡

. This
global index offers an inclusive perspective of damage across various locations without a specific crack path. Figure A1
compares the evolution of DI based on these two definitions as computed by various participants. In some cases, the
analysis is terminated at the failure time, while others continue to report the DI value for post-failure. Across all cases, it
can be observed that the length-based damage index (𝐷𝐼𝐿) consistently exceeds the area-based damage index (𝐷𝐼𝐴) at any
given time or equivalent intensity level. However, the specific ratio between the two indices varies based on the employed
damage modeling techniques.
Figure 11A presents the length-basedDI scenario, where themedian failure PGA stands at 0.33 g. Four fragility functions

are constructed corresponding to limit states (LS) of 1% (initiation), 5% (propagation), 10% (expansion), and 50% (critical)
of the length-based DI. All functions are fitted using theMLEmethod. The medians for these four LS values progressively
increase from 0.11, 0.17, 0.23 to 0.34 g, with dispersions varying at 0.8, 0.68, 0.67, and 0.83 respectively. These dispersion
values slightly surpass those of the displacement-based fragility functions and can be attributed to concrete modeling
complexities and the element size at the dam-foundation interface. Note that the crest displacement is a global response
and is less sensitive to mesh size and other local effects. The fitted distributions for length-based damage indices using the
GEV and lognormal distributions exhibit similarity, while the gamma distribution differs. The median of the lognormal
distribution is 23%, signifying that half of the models encounter 23% of the base crack at the time of failure.
In Figure 11B, the median failure PGA is 0.34 g for area-based DI. Three fragility functions correspond to LSs of 1%,

5%, and 10% of the area-based DI. All three functions exhibit a similar dispersion of approximately 0.95, considerably
larger than the displacement-based dispersion. Variability among results is exacerbated by the type of concrete constitutive
model andmesh size. A largermesh tends to create amore localized crack pattern, while a very finemesh leads to diffusive
(hair-like) crack propagation. The fitted distributions for area-based damage indices exhibit similarity, with the median of
the lognormal distribution at 8.5%. This indicates that half of the models fail before reaching 8.5% of overall body damage.
Figure A2 presents the evolution of crack propagation and failure modes of the dam across different times in the IAA

or equivalent PGA levels. In line with the variation observed in the DI and failure times, a substantial variability in crack

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1185

0 100 200
0

50

100

150

200

130 150 170 190
0

0.01

0.02

0.03

0 2 4 6
0

0.2

0.4

0.6
Data
GEV
Lognormal
Gamma

7 9 11 13
0

0.1

0.2

0.3

F IGURE 1 2 Processing modeling uncertainty results from the 12th ICOLD benchmark study.

propagation emerges, particularly concerning PGA (or time). The initiation of cracking is consistently observed at the heel
point and near the downstream neck discontinuity across nearly all cases. The base crack progresses horizontally towards
the dam’s toe, while the inclined propagation of the crack at the neck heads upstream. Unlike the base crack, a noticeable
pattern of multiple crack paths forms around the neck, with through cracks appearing before the base crack. Moreover, a
series of crack paths develop along the downstream face, primarily at the mid-height of the dam, propagating downwards
and towards the upstream face. Overall, the observed crack propagation in various cases aligns with findings from shake
table tests on gravity dams.111–113
Among the models displayed in Figure A2, eight models (#31 and #33 excluded) adopt a free-field foundation bound-

ary condition. However, specific conclusions regarding crack patterns based on foundation boundary conditions remain
elusive. Notably, five of these models rely on Abaqus software and concrete damage plasticity (#16, #26, #28, #31, and
#33). These models share a common FSI approach, albeit with variations in boundary conditions and mesh attributes.
Despite similar crack profiles, they exhibit differences in the time (or equivalent PGA) when the crack patterns manifest.
For example, the crack pattern observed at 𝑡 = 3 s (i.e., 0.17g) in model #16 is approximately replicated at 𝑡 =9 (for #26),
𝑡 =5 (for #28), 𝑡 =13 (for #31), and 𝑡 =5 (for #33).
Model #19 employs the Equivalent Linearization technique with large triangular elements (15 m), resulting in limited

accuracy in predicting crack paths and associated DI values. Models #17 and #27, despite sharing similar foundation
boundary conditions and concrete models, differ in their use of 2D plane strain and 3D slice model. While both models
exhibit a similar trend in crack pattern and progression, the extent of damage at the dam’s upper levels is greater in the
3D slice model, whereas the base crack at the dam-foundation interface is more pronounced in the 2D model.

6 DISCUSSION

Other ICOLD benchmark workshops have also investigated various case studies with the objective of identifying best
practices in seismic analysis of concrete dams. Notable examples include the 2013 workshop on seismic analysis of a
symmetry arch dam, the 2015 workshop on an asymmetry arch dam, and the 2017 workshop on seismic analysis of an
arch-gravity dam. This section aims to gather, refine, and quantify findings from these additional studies to facilitate
meaningful comparisons with the outcomes presented in this paper.

6.1 12th ICOLD benchmark results

The 12th ICOLD benchmark114 in 2013 revolved around a symmetric arch dam. Participants were provided with finite
element meshes of the dam and foundation, as well as material properties, a three-component ground motion record,
and suggested damping values. All simulations were conducted within the linear elastic analysis, without modeling any
concrete damage or joints. The primary focus of this benchmarkwas to perform dynamic FSI analyses. The 12 participants
submitted results, with 7 of them employing acoustic elements for fluid elements, while others adopted the added mass
approach. Notably, most participants utilized massless foundation model rather than free-field boundary conditions. The
collected data encompassed displacements and stresses at various locations. Figure 12 (left) visualizes the distribution of
static and dynamic (maximum and minimum) displacements along the crown cantilever in the stream direction. Data
points at the crest level were subjected to fitting with several distributional models. The computed dispersion value for
displacement was notably small, at 0.12, which aligns with the linear nature of the simulations and the limited modeling

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1186 HARIRI-ARDEBILI

variability presented by FSI. This figure also shows the fitted distributions to the maximum and minimum principal
stresses at the upstream point of the crest. It is important to note that this point was selected for uniformity of calculations,
despite not necessarily representing the highest stress value. The computed dispersion attributed to stress was measured
at 0.2.

6.2 13th ICOLD benchmark results

The 13th ICOLD benchmark115 conducted in 2017 focused on an asymmetric arch dam. Participants were furnished with
finite element meshes of the dam, foundation, reservoir information, material properties, three three-component ground
motion records, and proposed damping values. The dynamic properties distinct from static ones. The benchmark also
accounted for thermal gradients/analyses. While participants were advised to carry out linear simulations, the use of
dam-foundation joints was permissible. Participants were tasked with providing displacement results along the crown
cantilever, as well as the maximum and minimum principal stresses at 10 different locations: dam crest (left: A, center: B,
and right: C), dam upper part (left: D, right: E), abutment upper (left: F and right: G), abutment lower (left: H, right: I),
and dam lower part (J). Out of the total 11 participants, 5 of them employed compressible fluid elements in their analyses.
Upon processing the crest displacement results, a dispersion value of 0.25 was computed. However, the benchmark’s
complexity introduced numerous uncontrolled modeling factors, rendering the generalization of results challenging and
their reliability questionable.
Instead, Figure 13 illustrates the correlation matrix among the maximum principal stress values at 10 distinct locations.

The diagonal cells depict the histogram of ranks for each location, incorporating a non-parametric kernel-smoothing
distribution. Kernel distribution is employed due to the inability of a parametric distribution to adequately fit the data.
The lower triangular cells exhibit the discrete rank values for any two locations, denoted as 𝑙𝑖 and 𝑙𝑗 , alongside their local
regression using weighted linear least-squares and a 2nd-degree polynomial model. The upper triangular cells present the
Spearman’s linear correlation coefficient (R) in the format of 𝑅 = 1 −

6
∑

𝑑2

𝑛(𝑛2−1)
(where 𝑑 represents the difference between

the ranks of the two columns, and 𝑛 is the length of each column), along with 𝑝-values for assessing the hypothesis
of no correlation against the alternative hypothesis of a nonzero correlation (P). The correlation coefficient is further
distinguished by a distribution of colors, spanning from −0.41 (between locations D and F) to +0.93 (between locations
C and E). Consequently, excluding discussions of uncertainty in individual responses, this study highlights a substantial
correlation between the points at the crest and upper parts. As the data processing progresses towards the abutment and
lower part of the dam, the pattern and correlation of the results undergo notable changes (i.e., reduces in most cases).

6.3 14th ICOLD benchmark results

Finally, the 14th ICOLDbenchmark study116 has concentrated on the seismic analysis of an arch gravity dam. The provided
finite element model of the dam encompasses its foundation, with options for both free-field and non-free-field bound-
ary conditions. Concrete and rock material properties are provided, with the potential for nonlinearity being introduced
through dam-foundation interface joints. However, the stiffness values for these joints are not specified. Various types of
interface elements were employed, including frictional contact with a gap, cohesive law-based elements, spring elements,
joint elements, node-to-surface contact, and surface-to-surface contact. Themajority of simulations utilize an implicit time
integration scheme. Depending on the complexity of the modeling, different coupled system problems are formulated.
Foundation radiation boundary conditions encompass dashpot/Lysmer-Kuhlemeyer, perfectly matched layers, free-field,
and node-to-surface contact with a penalty.
While several parameters were provided to participants, the presence of numerous modeling variables beyond direct

control renders the UQ process challenging, particularly given the limited number of participants (ranging from 6 to 9
for different tasks). Three distinct seismic analyses are evaluated here, focusing on their relative variations rather than
providing generalized recommended values for future studies: Mod-I involves linear time-history analysis with simpli-
fied dynamic interactions (massless foundation), with the recommended use of generalized Westergaard’s added masses.
Mod-II entails nonlinear time-history analysis (applying uplift as external forces and allowing the interface elements to
model opening) with simplified dynamic interactions. Mod-III centers around linear time-history analysis with advanced
dynamic interactions (utilizing dynamic FSI and free-field foundation boundary conditions).
According to Figure 14, it becomes evident that utilizing Mod-III (linear analysis with advanced interaction models)

reduces the median and dispersion of the displacement response in comparison to Mod-I (linear analysis with simplified

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1187

0

10

20
(A)

0

5 (B)

0

5
(C)

0

5

10
(D)

0

5

10
(E)

0
1
2
3

(F)

0
1
2
3

(G)

0
2
4 (H)

-1
0
1
2 (I)

0 2 4 6
0

5

0 2 4 6 8 0 2 4 0 5 10 0 5 10 0 1 2 3 0 1 2 3 0 4 -2 0 2 0 2 4

(J)

F IGURE 13 Correlation among maximum principal stress values at 10 different locations collected from outcome collected from the
13th ICOLD benchmark study.

25 75 125
0

0.005

0.01

45 75 105 135
0

0.005

0.01

0.015

50 70 90
0

0.01

0.02

0.03

0.04 Data
GEV
Lognormal
Gamma

Mod-I Mod-II Mod-III

0

2

4

6

F IGURE 14 Processing modeling uncertainty results from the 14th ICOLD benchmark study.

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1188 HARIRI-ARDEBILI

interaction). Therefore, a bias term can be identified in simplified model compared to the advanced one. The employment
of the nonlinearmodel (Mod-II) increases themedian displacement in comparison to the linearmodel (Mod-I), but it does
not significantly impact the dispersion value. The logarithmic standard deviation (i.e., dispersion) for the threemodel types
(I, II, III) is measured at 0.39, 0.37, and 0.18, respectively. This figure also visualizes the results of tensile arch stresses at
the crest in the form of box-plots. The utilization of advanced interaction models (Mod-III) leads to a reduction in dam
stresses. Hence, it is apparent that the accurate modeling of the foundation-dam-reservoir coupled system, particularly
with respect to the interface joint, takes precedence. These results are pertinent to the specific ground motion record
provided for the study.

7 CONCLUSIONS, RECOMMENDATIONS, AND FUTUREWORK

Our comparative analysis, spanning a range of benchmark cases, emphasized the substantial impact of modeling-related
uncertainties on key dam response parameters: displacement, DI, stress distribution, crack propagation, and failure
modes. These modeling variations encompass factors such as mesh type, mesh size, time step used in numerical simula-
tion, convergence criteria, foundation model and associated boundary conditions for wave propagation and absorption,
fluid element selection and FSI approach, chosen concrete constitutive model and its damage response, time integration
scheme for stability control, and the software package for analysis.

7.1 Concluding remarks

Using various distributional models, including the lognormal, GEV, and gamma models, we have furthered our under-
standing. Notably, the lognormal distribution has allowed us to quantify the logarithmic standard deviation, referred to
as dispersion (𝛽𝑚𝑜𝑑). Its importance arises from its potential integration with other sources of uncertainty, such as ground
motion RTR variability (𝛽𝑅𝑇𝑅) and material randomness (𝛽𝑚𝑎𝑡), via SRSS combination.
By examining the dam response under the influence of a single component of the Taft groundmotion recordwith a PGA

of 0.18 g, the following high-level conclusions can be drawn. The computed 𝛽𝑚𝑜𝑑 (assuming a lognormal distribution) is
found to be 0.45, 0.30, 0.32, and 0.30 for the maximum dynamic crest displacement, maximum hydrodynamic pressure
at the heel, heel and crest maximum acceleration, respectively. The time-dependent standard deviation values exhibit
a nearly linear increasing pattern, occasionally punctuated by intermittent spikes along this trajectory. Conversely, the
pattern for the standard deviation of the relative crest displacement, absolute heel and crest acceleration, and heel hydro-
dynamic pressure differs. In this context, themaximum standard deviation is observed around 5–10 s, corresponding to the
strong motion part of the applied record, followed by a decrease to a smaller and relatively stable value for the remainder
of the simulation. The length-based DI results spanned from 0%–18%, with a median of 4.5%. Detailed discussions on the
contributions of different modeling uncertainties and assumptions are presented in the paper and are not reiterated here.
On the other hand, the IAA approach provided a unique opportunity to investigate the progression of failure in concrete

dams under various modeling strategies. This methodology has yielded a series of intensity-dependent capacity functions
for displacement and DI. Moreover, we introduced a novel fragility function that primarily reflects level III and IV mod-
eling variability. At the time of failure, the dispersion of PGA𝑓𝑎𝑖𝑙𝑢𝑟𝑒 was computed to be approximately 0.7, surpassing the
typical ground motion RTR validity of 0.4, highlighting the significance of modeling uncertainty in the failure analysis
of dams. Depending on the assumed limit state value, the dispersion of the DI may vary within the range of 0.67-0.83 for
the length-based criterion and around 0.95 for the area-based criterion. Such pronounced dispersion results from the low
level of confidence in damage analysis of concrete dams and the accurate documentation of the DI.

7.2 Recommendations

Based on the results of the four ICOLD benchmark studies examined in this paper (one in-depth and three dis-
cussed), coupled with insights from published materials on numerical simulations of concrete dams,1,117 several general
recommendations are proposed to foster confidence and credibility in modeling and computational simulations.
The incorporation of fluid-structure dynamic interaction is important for accurate hydrodynamic pressure determina-

tion on the dam’s upstream face. The utilization of compressible fluid elements is advised for advanced analyses. Utilizing
simplified models may inadequately replicate the actual hydrodynamic pressure acting on the dam face, particularly in
the presence of complex geometries and multi-component seismic excitation.

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1189

The dynamic interaction between the foundation and dam significantly influences safety assessment and potential
failure modes, particularly in the context of severe seismic scenarios. It is imperative to establish non-reflecting boundary
conditions for the foundation block to mitigate any manifestation of unrealistic response. Emphasizing the utilization of
free-field boundary conditions takes precedence, surpassing other considerations like system nonlinearity. This strategic
preference is grounded in recognizing the macro-scale impact of dam-foundation interaction across a broad spectrum of
seismic intensities, wherein system nonlinearity assumes a more micro-level significance, predominantly manifesting at
high seismic intensity levels.
In scenarios involving high seismic loading, introducing nonlinearity into the coupled system is imperative. Cohesive

crack models, capable of accounting for crack opening, closing, and sliding, should be employed at known crack or joint
locations. However, if the location of these features is not known beforehand, a continuum crackmodel should be applied
to the entire dam, allowing cracks to evolve based on the external load. Through a sequence of iterative simulations, dis-
crete cracks can subsequently replace smeared cracks at specific locations, providing amore direct control over the nature
of potential joint sliding and openings. It is also recommended to avoid a mixture of linear elastic and nonlinear regions
within the system. This can distort stress distribution and lead to unrealistic responses at the linear-nonlinear interface.
Factors such as mesh type and size, time step, convergence criteria, time integration scheme for stability control, damp-

ing model, and critical damping value all exert influence in seismic analyses of concrete dams, especially in the nonlinear
phase. The use of implicit schemes is favored for nonlinear analyses, despite their computational demands. Aspects such
as the dam’s stage construction, gradual reservoir impoundment, thermal stresses (especially in 3D models), and the
reservoir bottom wave reflection coefficient can impact initial stress states and structural results.

7.3 Limitations and future work

In this study, we explored modeling-related uncertainties in seismic analysis of concrete dams, examining a diverse range
of dam types andmodeling approaches through ICOLDbenchmarkworkshops. The variations in responses from different
modeling strategies underscore the complex interplay between choices and predictions, emphasizing the need for accurate
representation in dam behavior under seismic loads.
Limitations: Our study, while comprehensive, has some limitations. The assumed homogeneity in dam and founda-

tionmaterials might not fully capture the heterogeneity present, especially in older dams. Non-homogeneous foundations
with varying layer thicknesses pose challenges in wave propagation and deconvolution. For long canyon dams, spatial
variability in ground motions requires consideration for wave passage and incoherency. Seismic responses to multi-
component groundmotion differ from single-component motion, and the characteristics of input groundmotion, such as
long-duration or near-fault motion, directly influence structural responses.
Future Work: Reducing modeling uncertainties necessitates a robust validation and verification process for finite

element simulations. Future endeavors should concentrate on formulating a comprehensive roadmap for conduct-
ing structural analysis of concrete dams, incorporating a rigorous verification and validation framework. This effort
should be complemented by extensive laboratory experiments and field measurements. The verification process must
encompass code verification, solution verification, and the distinction between numerical errors and uncertainties. It
is crucial to emphasize that validation should always follow verification to mitigate the accumulation of errors within
computational models.
Practical Implications: The study’s outcomes carry practical implications for seismic analysis in dam engineer-

ing. Our quantification of modeling uncertainties underscores the need for rigorous practices and parameter selection
to ensure reliable outcomes. Engineers and practitioners should prioritize validation and verification in finite element
simulations, considering the complexities introduced by modeling, materials and ground motion characteristics.

ACKNOWLEDGMENTS
The author acknowledges the contributions of Dr Jerzy Salamon (US Bureau of Reclamation), Dr Richard Malm
(formerly in KTH Royal Institute of Technology), and Ms Giorgia Faggiani (Ricerca Sistema Energetico) in formulating
the ICOLD benchmark problem. Gratitude is also extended to all workshop participants whose invaluable insights and
contributions have greatly enriched our understanding of modeling variability in seismic analysis of concrete dams. The
author extends gratitude to Dr Siamak Sattar (National Institute of Standards and Technology) for general discussion on
modeling uncertainty quantification, and benchmark data processing. The presented materials and conclusions in this
paper solely represent the author’s studies and opinions and do not reflect the views or opinions of others in any way.

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



1190 HARIRI-ARDEBILI

CONFL ICT OF INTEREST STATEMENT
The author declares no potential conflict of interests.

DATA AVAILAB IL ITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

REFERENCES
1. Saouma VE, Hariri-Ardebili MA. Aging, Shaking, and Cracking of Infrastructures: From Mechanics to Concrete Dams and Nuclear

Structures. Springer; 2021.
2. Roache PJ. Verification and Validation in Computational Science and Engineering. Vol 895. Hermosa Albuquerque; 1998.
3. Roy CJ, Oberkampf WL. A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing.

Comput Meth Appl Mech Eng. 2011;200(25-28):2131-2144.
4. Ewing ME, Liechty BC, Black DL. A general methodology for uncertainty quantification in engineering analyses using a credible

probability box. J Verif Valid Uncertain Quantif. 2018;3(2):021003.
5. Ayyub BM, Chao R. Uncertainty modeling in civil engineering with structural and reliability applications. Uncertainty modeling and

analysis in civil engineering. CRC Press; 1998:1-8.
6. Bukenya P, Moyo P, Beushausen H, Oosthuizen C. Health monitoring of concrete dams: a literature review. J Civ Struct Health Monit.

2014;4(4):235-244.
7. Alves S, Hall J. System identification of a concrete arch dam and calibration of its finite element model. Earthquake Eng Struct Dyn.

2006;35(11):1321-1337.
8. Harris DW, Snorteland N, Dolen T, Travers F. Shaking table 2-D models of a concrete gravity dam. Earthquake Eng Struct Dyn.

2000;29(6):769-787.
9. Uchita Y, Shimpo T, Saouma V. Dynamic centrifuge tests of concrete dams. Earthquake Eng Struct Dyn. 2005;34:1467-1487.
10. Bolzon G, Sterpi D, Mazzà G, Frigerio A. Numerical Analysis of Dams: Proceedings of the 15th ICOLD International Benchmark Workshop.

Vol 91. Springer Nature; 2020.
11. Salamon J. Evaluation of Numerical Models and Input Parameters in the Analysis of Concrete Dams. Tech. Rep. DSO-19-13. US Bureau of

Reclamation; 2018.
12. Hariri-Ardebili MA, Seyed-Kolbadi S, Kianoush M. FEM-based parametric analysis of a typical gravity dam considering input excitation

mechanism. Soil Dyn Earthquake Eng. 2016;84:22-43.
13. Sevieri G. The Seismic Assessment of Existing Concrete Gravity Dams: FE Model Uncertainty Quantification and Reduction. PhD thesis.

University of Braunschweig and University of Florence; 2019.
14. Poul MK. Development of Viscoelastic PML and Deconvolution of Input Motion at Depth for Soil-Structure Interaction–Application to Dams.

PhD thesis. Drexel University; 2019.
15. Jalayer F, Iervolino I, Manfredi G. Structural modeling uncertainties and their influence on seismic assessment of existing RC structures.

Struct Saf. 2010;32(3):220-228.
16. Gokkaya BU, Baker JW, Deierlein GG. Quantifying the impacts of modeling uncertainties on the seismic drift demands and collapse risk

of buildings with implications on seismic design checks. Earthquake Eng Struct Dyn. 2016;45(10):1661-1683.
17. Haindl M, Burton H, Sattar S. Quantification of equivalent strut modeling uncertainty and its effects on the seismic performance of

masonry infilled reinforced concrete frames. J Earthquake Eng. 2023:1-21.
18. Li H, Li L, Zhou G, Xu L. Effects of various modeling uncertainty parameters on the seismic response and seismic fragility estimates of

the aging highway bridges. Bull Earthquake Eng. 2020;18:6337-6373.
19. O’ReillyGJ, SullivanTJ.Quantification ofmodelling uncertainty in existing ItalianRC frames.EarthquakeEngStructDyn. 2018;47(4):1054-

1074.
20. Sattar S, Liel AB,Martinelli P.Quantification ofModelingUncertainties Based on the BlindPredictionContest Submissions. ASCE; 2013:1997-

2008.
21. Parisse F, Cattari S, Marques R, et al. Benchmarking the seismic assessment of unreinforced masonry buildings from a blind prediction

test. Structures. 2021;31:982-1005.
22. Tomić I, Penna A, DeJong M, et al. Shake-table testing of a stone masonry building aggregate: overview of blind prediction study. Bull

Earthquake Eng. 2023:1-43.
23. VassiliouMF, BroccardoM, Cengiz C, et al. Shake table testing of a rocking podium: Results of a blind prediction contest. Earthquake Eng

Struct Dyn. 2021;50(4):1043-1062.
24. Folz B, Filiatrault A. Blind predictions of the seismic response of a woodframe house: an international benchmark study. Earthq spectra.

2004;20(3):825-851.
25. Bartoli G, Betti M, Biagini P, et al. Epistemic uncertainties in structural modeling: a blind benchmark for seismic assessment of slender

masonry towers. Bull Earthquake Eng. 2017;31(5):04017067.
26. Winkler RL. Uncertainty in probabilistic risk assessment. Reliab Eng Syst Saf. 1996;54(2-3):127-132.
27. Ellingwood B. Issues related to structural aging in probabilistic risk assessment of nuclear power plants. Reliab Eng Syst Saf. 1998;62(3):171-

183.

 10969845, 2024, 3, D
ow

nloaded from
 https://onlinelibrary.w

iley.com
/doi/10.1002/eqe.4064 by U

niversity O
f M

aryland, W
iley O

nline L
ibrary on [02/07/2024]. See the T

erm
s and C

onditions (https://onlinelibrary.w
iley.com

/term
s-and-conditions) on W

iley O
nline L

ibrary for rules of use; O
A

 articles are governed by the applicable C
reative C

om
m

ons L
icense



HARIRI-ARDEBILI 1191

28. Parry GW. The characterization of uncertainty in probabilistic risk assessments of complex systems. Reliab Eng Syst Saf. 1996;54(2-3):119-
126.

29. Helton JC, Johnson JD, Oberkampf WL, Sallaberry CJ. Representation of analysis results involving aleatory and epistemic uncertainty.
Int J Gen Syst. 2010;39(6):605-646.

30. Ang AHS. Extended Reliability Basis of Structural Design Under Uncertainties. Tech. Rep., SAE Technical Paper; 1970.
31. Der-Kiureghian A, Ditle