Analytical and direct FE methods for free-field foundation modeling in 3D analysis of concrete dams

Three-dimensional (3D) modeling is essential for reliable seismic assessment of concrete dams, as it captures canyon topography, monolith interactions, and joint nonlinearity that cannot be represented in two-dimensional (2D) analyses. This study systematically evaluates the seismic performance of a benchmark concrete gravity dam through comparative 3D and 2D finite element (FE) analyses, investigating two free-field foundation modeling approaches: (i) a direct FE method that explicitly accounts for canyon geometry, and (ii) an analytical approach employing simplified flat foundation assumptions. Using the domain reduction method, the nonlinear seismic analysis of the Pine Flat Dam implements effective earthquake forces as seismic input at the absorbing boundaries. The results reveal that while 2D analyses show methodological agreement (both relying on one-dimensional free-field assumptions), they fundamentally cannot reproduce 3D effects such as stress amplification from canyon geometry and monolith interactions. The direct FE method proves superior in predicting joint sliding and opening responses, as well as canyon amplification effects, whereas simplified approaches yield non-conservative estimates due to its failure to account for canyon-induced wave scattering. This study concludes that 2D modeling may be sufficient for preliminary analyses, but final design verification requires 3D simulation with direct FE free-field representation, particularly for dams located in irregular topographies.

1 Introduction

Precise evaluation of seismic response in concrete dams is critically essential, given that structural failure during seismic events may result in devastating impacts on downstream populations and civil infrastructure (Hariri-Ardebili, 2018). Historical events, such as the damage to the Koyna gravity dam in India (1967), the Hsingfengkiang buttress dam in China (1962), and the Pacoima arch dam in the United States (1974), highlight the critical need for rigorous seismic vulnerability assessments and effective mitigation measures (Wang et al., 2020; Akpinar et al., 2025). On the other hand, the seismic response of the dam is highly influenced by dynamic interactions with the foundation rock and the reservoir water. Irregular foundation topography can significantly alter seismic wave propagation, leading to scattering and refraction that result in nonuniform motion along the dam-foundation interface, as evidenced by recorded earthquake data (Nowak and Hall, 1990; Alves and Hall, 2006; Chopra, 2012). Recent studies, including those from ICOLD benchmark workshops, have emphasized the variability of response quantities by using different foundation modeling assumptions, and also the importance of free-field modeling in capturing these interactions (Hariri-Ardebili, 2024). Free-field modeling is particularly critical for systems with irregular canyon topography, where wave propagation patterns are complex and highly variable (Zhao et al., 1993). Despite its significance, implementing free-field boundary conditions in three-dimensional (3D) seismic analyses of concrete dams remains limited, especially for systems involving cross-stream vibrations and contraction joint behavior.

To facilitate dynamic analysis, the massless foundation method has been widely used for foundation modeling (Akköse and Şimşek, 2010; Wang et al., 2012; Mandal and Maity, 2018; Jin et al., 2019; Faggiani et al., 2021). Nevertheless, this approach has two major drawbacks: it overlooks the radiation damping caused by the unbounded foundation and does not consider the amplification of seismic waves due to topography. Fry and Matsumoto (Fry and Matsumoto, 2018) compared finite element (FE) modeling of concrete dams using the massless foundation method with earthquake records and experimental measurements, revealing a significant overestimation of surface motion and vibrations induced in the dam. Similarly (Abbasiverki et al., 2021), demonstrated that the massless method leads to unrealistic responses in the seismic behavior of concrete buttress dams under high-frequency excitations, particularly when topographic amplification is neglected. These findings highlight the need for more advanced modeling approaches that incorporate the mass of the foundation and accurately represent the semi-infinite behavior of the foundation domain.

Including the mass of the foundation in dynamic analyses requires appropriate boundary conditions to replicate the semi-infinite behavior of the foundation and incorporate relevant free-field motions. Neglecting free-field boundary conditions can result in wave leakage, where absorbing boundaries attenuate incoming waves propagating through the foundation domain (Nielsen, 2006). The general principle in soil-structure interaction analyses is that the free-field system must match the actual system in the region exterior to the absorbing boundaries (Wolf, 1989; Bielak et al., 2003). For foundations with irregular surfaces, such as those featuring canyons, the free-field system is inherently three-dimensional. Rigorous methods for modeling such systems often rely on perfectly matched layers (PML) (Basu and Chopra, 2004) as absorbing boundaries and the Domain Reduction Method (DRM) for the free-field system (Mohammadnezhad et al., 2019a; Zhang et al., 2019). DRM, first formulated by Bielak et al. (Bielak et al., 2003), employs a two-stage finite element framework to model seismic wave propagation through the application of effective earthquake forces. This theoretically robust approach, however, faces significant implementation challenges in commercial finite element environments, where combined support for DRM and PML boundary conditions remains exceptionally limited. Among mainstream commercial codes, only LS-DYNA currently provides a comprehensive implementation of this coupled methodology.

To address these limitations, Løkke and Chopra (Løkke and Chopra, 2018) proposed a modified DRM approach incorporating viscous boundary conditions, enabling the calculation of effective earthquake forces using the direct FEM and their application to non-reflective boundaries. While this method can be implemented in any FE code, it requires extensive data transfer and bookkeeping, particularly for 3D systems. Although effective earthquake forces can be automatically calculated with free-field boundary elements (Nielsen, 2014), only a few FE codes, such as Plaxis, support such elements. Mohammadnezhad et al. (Mohammadnezhad et al., 2019b) developed a free-field boundary condition method based on the direct FE approach, which eliminates the need for preliminary analysis or subroutines and can be implemented in most FE software. However, this method was limited to 2D foundations with uniform surfaces.

A simplified free-field foundation system assumes a flat foundation, ignoring the effects of canyon topography in the free-field analysis. The proposed approach can be analytically investigated using one-dimensional (1D) wave propagation principles (He et al., 2010). While its initial formulation and validation focused on viscous-spring boundaries in 2D foundations (Chuhan et al., 2009; Chen et al., 2012), Song et al. (Song et al., 2018) later generalized it to infinite absorbing elements. However, its accuracy has been evaluated only for foundations with uniform surfaces, leaving its applicability to more complex, irregular topographies uncertain. To clarify the overall modeling strategy and situate the present work within the broader context of dam-foundation analysis, Figure 1 provides a schematic flowchart outlining the evolution of modeling approaches and the specific position of the methods investigated in this study.

Figure 1

Figure 1. Evolution of dam-foundation modeling approaches from massless foundations to advanced radiation-damping methods.

In the design of concrete gravity dams, the structure is segmented into independent monoliths separated by vertical joints to account for thermal expansion and contraction. Under the assumption that the monoliths behave independently during an earthquake, seismic analysis of gravity dams is often conducted using 2D modeling of the tallest non-overflow monolith (Salamon, 2018; Salamon et al., 2019; Poul and Zerva, 2019). However, to accurately capture the seismic response of gravity dams, it is essential to employ 3D FE modeling that includes the dam, foundation, and reservoir. This comprehensive approach is necessary due to the significant influence of contraction joints and cross-stream vibrations on the dynamic response of the dam (Azmi and Paultre, 2002; Arici et al., 2014; Wang et al., 2017). Moreover, the shape of the canyon in which a gravity dam is situated significantly affects its seismic response (Bybordiani and Arıcı, 2017; Nikkhakian et al., 2020). Despite these complexities, the role of free-field foundation boundaries in concrete gravity dam performance has only been explored using two-dimensional finite element models (Chen et al., 2021; Enzell et al., 2021).

This study presents a systematic investigation into free-field foundation modeling approaches for the seismic analysis of concrete gravity dams, using the Pine Flat Dam as a case study (ICOLD workshop case study). The research evaluates two distinct methodologies: (i) the direct FE approach developed by Løkke and Chopra (Løkke and Chopra, 2018), which explicitly incorporates 3D canyon topography effects through rigorous FE discretization; and (ii) the analytical framework proposed by Song et al. (Song et al., 2018), based on simplified 1D wave propagation theory. The computational framework incorporates nonlinear finite element analysis of contraction joint behavior under combined opening–sliding mechanisms and cross-stream seismic excitation. An advanced seismic wave deconvolution technique is implemented to ensure accurate representation of input motions at the foundation boundary.

Through a comprehensive comparison of these methodologies in both 2D and 3D modeling contexts, this paper provides new insights into the critical aspects of seismic wave propagation and foundation representation for dam seismic analysis. This comparison is particularly significant because a large body of published work in the past decade on concrete dams–including probabilistic analyses, incremental dynamic analysis, fragility curve development, comparisons of near- and far-field motions, short- and long-duration effects, and seismic risk assessments–has relied predominantly on either massless foundation models or simplified 1D wave propagation approaches. The current study offers a necessary benchmark on the validity of results obtained under such assumptions or simplifications, highlighting their potential implications for decision-making in dam engineering projects.

2 Dam-reservoir-foundation system

Hydrodynamic pressures during seismic events are primarily generated by the dynamic motion of the dam and foundation boundaries, while the interaction with the reservoir and foundation influences their distribution and magnitude. Simplified methods, such as the added mass approach, which neglect water compressibility, can result in unreliable seismic response predictions. More accurate models that consider water compressibility and dynamic interactions are essential for reliable seismic analysis of dams (Ameen et al., 2017; Salamon, 2011). The acoustic fluid method can be effectively employed to account for the compressibility of water accurately in the numerical modeling of reservoirs. This approach is particularly useful for capturing compressibility effects that are significant in various fluid–structure interaction scenarios (Goldgruber, 2015; Faggiani and Masarati, 2018; Mazzieri et al., 2022). The fluid-structure interaction in the reservoir can be modeled using the acoustic wave equation under the following assumptions: (i) fluid compressibility with irrotational flow, (ii) negligible viscous effects, and (iii) absence of mean flow or body forces. The acoustic formulation employs pressure as the sole degree of freedom for fluid elements. Within the system of dam-reservoir-foundation, boundary accelerations generate hydrodynamic pressures in the reservoir (Abbasiverki, 2023). The different boundary conditions applied to the reservoir are depicted in Figure 2.

Figure 2

Figure 2. Schematic of the dam-reservoir-foundation finite element model configuration, including reservoir interface boundary conditions (Abbasiverki, 2023).

The dynamic behavior of the coupled dam-reservoir-foundation system is modeled using finite elements with truncated boundary domains. This system can be expressed through the following matrix formulation (Equation 1), which accounts for the interaction between the structural and fluid domains:

$\begin{array}{ll}\hfill & \left[\begin{array}{cc}\hfill \mathbf{m}\hfill & \hfill \mathbf{0}\hfill \\ \hfill {\rho }_r\left({\mathbf{N}}_h^{\mathsf{T}}+{\mathbf{N}}_b^{\mathsf{T}}\right)\hfill & \hfill \mathbf{b}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\ddot{\mathbf{r}}}^t\hfill \\ \hfill {\ddot{\mathbf{p}}}_r^t\hfill \end{array}\right\}+\left[\begin{array}{cc}\hfill \mathbf{d}+{\mathbf{c}}_f\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathbf{h}+{\mathbf{c}}_r\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\dot{\mathbf{r}}}^t\hfill \\ \hfill {\dot{\mathbf{p}}}_r^t\hfill \end{array}\right\}\hfill \\ \hfill & +\left[\begin{array}{cc}\hfill \mathbf{k}-\left({\mathbf{N}}_h+{\mathbf{N}}_b\right)\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathbf{s}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\mathbf{r}}^t\hfill \\ \hfill {\mathbf{p}}_r^t\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill {\mathbf{f}}^0\hfill \\ \hfill \mathbf{0}\hfill \end{array}\right\}\hfill \end{array}(1)$

Within this framework, the structural properties of the dam and its supporting foundation are characterized by the mass $(\mathbf{m})$, damping $(\mathbf{d})$, and stiffness $(\mathbf{k})$ matrices. Conversely, the fluid domain is characterized by corresponding matrices ($\mathbf{b}$, $\mathbf{h}$, $\mathbf{s}$), representing its mass, damping, and stiffness properties, respectively. Energy absorption at artificial boundaries is incorporated through the damping matrices ${\mathbf{c}}_f$ (foundation) and ${\mathbf{c}}_r$ (reservoir). The system’s dynamic response is governed by two primary variables: ${\mathbf{r}}^t$ quantifying structural displacements and ${\mathbf{p}}_r^t$ defining fluid pressures. Interface phenomena between different media are modeled via coupling terms ${\mathbf{N}}_b$ (reservoir-foundation interaction) and ${\mathbf{N}}_h$ (dam-reservoir interaction). The effective earthquake forces ${\mathbf{f}}^0$ at the foundation boundaries serve as system excitation. Their computation is presented in the following section.

3 Free-field boundary condition

The computation of effective earthquake forces in this study relies on the fundamental principle of separating the seismic input from the local structural response. The assumption of linear elastic material behavior for the domain exterior to the absorbing boundaries is a fundamental component of the Domain Reduction Method (DRM) and similar seismic input techniques (Bielak et al., 2003; Løkke and Chopra, 2018). This approach conceptually separates the problem into two parts:

1. The determination of the free-field input motion, which is the wave field that would exist at the site in the absence of the structure (the dam). This is calculated assuming the seismic waves propagate through a linear elastic, semi-infinite medium (or horizontally layered media) from the seismic source.

2. The nonlinear structural response of the dam-foundation-reservoir system to this input motion.

The justification for this separation is twofold. First, the incident seismic wave field propagating from the bedrock base through the far-field geological structure is typically within the linear elastic range due to high confining pressures. Second, and most critically, the purpose of the free-field system is not to model the complex, nonlinear dam-foundation interaction, but to provide a physically correct input signal at the artificial boundaries of the local model. The nonlinear phenomena of primary interest—such as joint sliding and opening—are localized to the dam and its immediate near-field foundation. These occur in response to the applied free-field motion. Modeling the far-field as linear elastic is therefore a standard and valid simplification that allows for the efficient and accurate application of seismic forces while capturing the essential nonlinear behavior within the local region of interest.

It should be noted, however, that nonlinearity can develop in the near-field foundation. Previous research by Hariri-Ardebili (Hariri-Ardebili, 2014) has demonstrated the potential for foundation rock cracking and damage propagation under severe seismic loading, which can influence the overall structural response. A detailed investigation of this coupled dam-foundation nonlinearity is a recommended topic for future work.

The computation of effective earthquake forces employs two distinct free-field foundation models, both relying on the vertical wave propagation assumption¹: (i) an actual free-field system that includes a canyon, as shown in Figure 3a, and (ii) a simplified free-field system with a flat foundation surface (neglecting the canyon), as shown in Figure 3b. The former system is simulated through direct finite element computation, with the latter being treated analytically via a 1D wave propagation formulation.

Figure 3

Figure 3. Free-field foundation system (adapted from Løkke and Chopra, 2019): (a) With canyon, (b) Without canyon.

3.1 Direct FE method

The direct finite element method evaluates seismic forces at foundation boundaries through two components: free-field tractions and velocity-dependent damping. For lateral boundaries with infinite elements, the effective earthquake force combines these effects as shown in Equation 2:

${\mathbf{f}}^0={\mathbf{x}}^0+{\mathbf{c}}_f{\dot{\mathbf{r}}}^0(2)$

where ${\mathbf{x}}^0$ represents boundary tractions and ${\mathbf{c}}_f{\dot{\mathbf{r}}}^0$ accounts for damping, with ${\mathbf{c}}_f$ containing viscous coefficients and ${\dot{\mathbf{r}}}^0$ being nodal velocities from free-field analysis. The base boundary formulation simplifies to Equation 3:

${\mathbf{f}}^0=2{\mathbf{c}}_f{\dot{\mathbf{r}}}_I^0(3)$

where ${\dot{\mathbf{r}}}_I^0$ denotes upward-propagating waves, derived as half the outcrop motion (Poul and Zerva, 2018). The outcrop motion (theoretical free-surface displacement at equivalent depth) doubles the incident wave amplitude due to free-surface reflection requirements, obtained through surface motion deconvolution. The free-field velocity components ${\dot{\mathbf{r}}}^0$ and boundary tractions ${\mathbf{x}}^0$ are obtained through a preliminary analysis of the decoupled free-field system prior to the full dam-reservoir-foundation simulation. For two-dimensional foundations, this reduces to analyzing a single 1D finite element column. In three-dimensional configurations, the free-field system (Figure 4a) employs a simplified 2D representation supplemented by 1D corner columns, as shown in Figure 4b. The lateral boundary nodal forces computed via Equation 2 incorporate velocity and traction results from 1D column analyses along canyon boundaries, while the upstream and downstream cross-stream boundaries utilize solutions from the 2D system.

Figure 4

Figure 4. Schematic representation of the direct finite element approach for computing effective earthquake forces, following the methodology of (Løkke and Chopra, 2018): (a) Three-dimensional free-field system showing uniform canyon excavation in the foundation rock half-space, (b) Two-dimensional simplification with auxiliary 1D corner columns.

3.2 Analytical method

The free-field modeling methodology developed by He et al. (He et al., 2010), based on 1D wave propagation theory for viscous-spring boundaries, determines effective earthquake forces via Equation 4:

${\mathbf{f}}^0={\mathbf{k}}_f{\mathbf{r}}^0+{\mathbf{c}}_f{\dot{\mathbf{r}}}^0+{\mathit{\sigma }}^0\mathbf{n}(4)$

with ${\mathbf{k}}_f$ as system stiffness, ${\mathbf{r}}^0={[uvw]}^T$ and ${\dot{\mathbf{r}}}^0={[\dot{u}\dot{v}\dot{w}]}^T$ as boundary displacement/velocity vectors, ${\mathit{\sigma }}^0$ as free-field stress tensor, and $\mathbf{n}$ as the boundary normal vector.

Song et al. (2018) extended this approach to infinite elements, incorporating viscous dampers that are uniformly embedded within the dynamic infinite elements. In this formulation, the spring coefficient is neglected (${\mathbf{k}}_f=\mathbf{0}$), and the effective earthquake forces at each artificial boundary surface are derived using 1D wave propagation theory.

The specific expressions for these effective earthquake forces in the 3D FE model are given in Equations 5–9 below, with all relevant parameters clearly annotated in the schematic representation shown in Figure 5:

Figure 5

Figure 5. Schematic of the 3D free-field foundation model according to the analytical method, showing orthogonal dampers and effective earthquake forces on all boundary surfaces.

Bottom surface

$\left\{\begin{array}{c}\hfill f_x^{z-}\hfill \\ \hfill f_y^{z-}\hfill \\ \hfill f_z^{z-}\hfill \end{array}\right\}=A_b\left\{\begin{array}{c}\hfill 2C_{\text{BT}}{\dot{u}}_0\left(t\right)\hfill \\ \hfill 2C_{\text{BT}}{\dot{v}}_0\left(t\right)\hfill \\ \hfill 2C_{\text{BN}}{\dot{w}}_0\left(t\right)\hfill \end{array}\right\}(5)$

Surface in X negative direction

$\left\{\begin{array}{c}\hfill f_x^{x-}\hfill \\ \hfill f_y^{x-}\hfill \\ \hfill f_z^{x-}\hfill \end{array}\right\}=A_b\left\{\begin{array}{c}\hfill C_{\text{BN}}\left[{\dot{u}}_0\left(t-\frac{h}{c_s}\right)+{\dot{u}}_0\left(t-\frac{2H-h}{c_s}\right)\right]+\frac{\lambda }{c_p}\left[{\dot{w}}_0\left(t-\frac{h}{c_p}\right)-{\dot{w}}_0\left(t-\frac{2H-h}{c_p}\right)\right]\hfill \\ \hfill C_{\text{BT}}\left[{\dot{v}}_0\left(t-\frac{h}{c_s}\right)+{\dot{v}}_0\left(t-\frac{2H-h}{c_s}\right)\right]\hfill \\ \hfill C_{\text{BT}}\left[{\dot{w}}_0\left(t-\frac{h}{c_p}\right)+{\dot{w}}_0\left(t-\frac{2H-h}{c_p}\right)\right]+\rho c_s\left[{\dot{u}}_0\left(t-\frac{h}{c_p}\right)-{\dot{u}}_0\left(t-\frac{2H-h}{c_p}\right)\right]\hfill \end{array}\right\}(6)$

Surface in X positive direction

$\left\{\begin{array}{c}\hfill f_x^{x+}\hfill \\ \hfill f_y^{x+}\hfill \\ \hfill f_z^{x+}\hfill \end{array}\right\}=A_b\left\{\begin{array}{c}\hfill C_{\text{BN}}\left[{\dot{u}}_0\left(t-\frac{h}{c_s}\right)+{\dot{u}}_0\left(t-\frac{2H-h}{c_s}\right)\right]-\frac{\lambda }{c_p}\left[{\dot{w}}_0\left(t-\frac{h}{c_p}\right)-{\dot{w}}_0\left(t-\frac{2H-h}{c_p}\right)\right]\hfill \\ \hfill C_{\text{BT}}\left[{\dot{v}}_0\left(t-\frac{h}{c_s}\right)+{\dot{v}}_0\left(t-\frac{2H-h}{c_s}\right)\right]\hfill \\ \hfill C_{\text{BT}}\left[{\dot{w}}_0\left(t-\frac{h}{c_p}\right)+{\dot{w}}_0\left(t-\frac{2H-h}{c_p}\right)\right]-\rho c_s\left[{\dot{u}}_0\left(t-\frac{h}{c_p}\right)-{\dot{u}}_0\left(t-\frac{2H-h}{c_p}\right)\right]\hfill \end{array}\right\}(7)$

Surface in Y negative direction

$\left\{\begin{array}{c}\hfill f_x^{y-}\hfill \\ \hfill f_y^{y-}\hfill \\ \hfill f_z^{y-}\hfill \end{array}\right\}=A_b\left\{\begin{array}{c}\hfill C_{\text{BT}}\left[{\dot{u}}_0\left(t-\frac{h}{c_s}\right)+{\dot{u}}_0\left(t-\frac{2H-h}{c_s}\right)\right]\hfill \\ \hfill C_{\text{BN}}\left[{\dot{v}}_0\left(t-\frac{h}{c_s}\right)+{\dot{v}}_0\left(t-\frac{2H-h}{c_s}\right)\right]+\frac{\lambda }{c_p}\left[{\dot{w}}_0\left(t-\frac{h}{c_p}\right)-{\dot{w}}_0\left(t-\frac{2H-h}{c_p}\right)\right]\hfill \\ \hfill C_{\text{BT}}\left[{\dot{w}}_0\left(t-\frac{h}{c_p}\right)+{\dot{w}}_0\left(t-\frac{2H-h}{c_p}\right)\right]+\rho c_s\left[{\dot{v}}_0\left(t-\frac{h}{c_s}\right)-{\dot{v}}_0\left(t-\frac{2H-h}{c_s}\right)\right]\hfill \end{array}\right\}(8)$

Surface in Y positive direction

$\left\{\begin{array}{c}\hfill f_x^{y+}\hfill \\ \hfill f_y^{y+}\hfill \\ \hfill f_z^{y+}\hfill \end{array}\right\}=A_b\left\{\begin{array}{c}\hfill C_{\text{BT}}\left[{\dot{u}}_0\left(t-\frac{h}{c_s}\right)+{\dot{u}}_0\left(t-\frac{2H-h}{c_s}\right)\right]\hfill \\ \hfill C_{\text{BN}}\left[{\dot{v}}_0\left(t-\frac{h}{c_s}\right)+{\dot{v}}_0\left(t-\frac{2H-h}{c_s}\right)\right]-\frac{\lambda }{c_p}\left[{\dot{w}}_0\left(t-\frac{h}{c_p}\right)-{\dot{w}}_0\left(t-\frac{2H-h}{c_p}\right)\right]\hfill \\ \hfill C_{\text{BT}}\left[{\dot{w}}_0\left(t-\frac{h}{c_p}\right)+{\dot{w}}_0\left(t-\frac{2H-h}{c_p}\right)\right]-\rho c_s\left[{\dot{v}}_0\left(t-\frac{h}{c_s}\right)-{\dot{v}}_0\left(t-\frac{2H-h}{c_s}\right)\right]\hfill \end{array}\right\}(9)$

4 Case study description

This study examines the Pine Flat concrete gravity dam as a case study. Situated on the King’s River near Fresno, California, USA, the dam was constructed in 1954. Its crest extends 560.80 m and comprises 36 monoliths, each 15.24 m wide, along with an additional monolith measuring 12.2 m in width. The tallest non-overflow section, Monolith 16, is analyzed in this work, with its cross-sectional geometry illustrated in Figure 6 (Rea et al., 1972).

Figure 6

Figure 6. Illustration of pine flat dam (Rea et al., 1972): (a) Downstream view, (b) Cross section geometry of the tallest monolith.

The concrete used in the dam has the following material properties: elastic modulus $E_c=30$ GPa, mass density ${\rho }_c=2400$ kg/m³, and Poisson’s ratio ${\nu }_c=0.2$. The analysis assumes a perfect bonding between the concrete structure and bedrock foundation, with identical material properties assigned to the rock ($E_f=30$ GPa, ${\rho }_f=2400$ kg/m³, ${\nu }_f=0.2$) (Salamon et al., 2019). The reservoir water is characterized by a bulk modulus $K_r=2.02$ GPa and density ${\rho }_r=1000$ kg/m³.

A material damping ratio of 5% was assigned to both the dam and the foundation. This value represents a conventional, conservative estimate for the inherent hysteretic energy dissipation within the concrete and rock. While recent field measurements suggest that the overall damping of dam-water-foundation systems typically ranges from 1% to 5% (Løkke and Chopra, 2019), the 5% value was selected for the material component to account for inherent energy dissipation in this study. For the time-domain analysis, Rayleigh damping is implemented using only a stiffness-proportional damping coefficient ($\beta =0.0046$ s). Mass-proportional damping is explicitly excluded $(\alpha =0)$ to prevent the generation of artificial damping forces proportional to absolute velocity, which could otherwise artificially restrain joint opening and sliding behavior—a critical aspect of this study (El-Aidi and Hall, 1989; Hall, 2006).

The stiffness-proportional coefficient $\beta$ was calibrated to achieve the target damping ratio $(\zeta =0.05)$ at the dam’s fundamental frequency. The calibration was based on the experimentally determined fundamental frequency from Rea et al. (Rea et al., 1972) ($f_1=3.47$ Hz). This approach was preferred because the linear FE model, which omits contraction joints, predicts a higher fundamental frequency (3.78 Hz) due to its overestimation of structural stiffness. The coefficient was calculated using the standard formula for a stiffness-proportional-only model as shown in Equation 10:

$\beta =\frac{2\zeta }{{\omega }_1}(10)$

where ${\omega }_1=2\pi f_1$ is the circular frequency. This calibration ensures the model possesses energy dissipation characteristics that are representative of the actual dam’s dynamic behavior.

4.1 Finite element model

For 3D FE modeling of Pine Flat Dam, the dam-foundation interface is contoured based on the downstream view (Figure 6a) and assumed constant in the upstream-downstream direction (Figure 7a). The stream-direction cross-section follows the configuration in Figure 6b, and the corresponding 2D FE model is constructed using this cross-section (Figure 7b). Figure 7c presents the complete 3D FE model, detailing both the dam structure and contraction joint arrangement. The discretization comprises 13,383 C3D8R solid elements, with the largest element dimension constrained to be smaller than $\frac{1}{8}$ of the minimum shear wavelength (Lysmer and Kuhlemeyer, 1969). The two-dimensional finite element representation (Figure 7b) is discretized with 201 CPE4R plane-strain elements.

Figure 7

Figure 7. "Illustration of the geometry of the Pine Flat dam: (a) 3D FE model, (b) 2D FE model, (c) Layout of the contraction joints.

Contraction joints between monoliths are explicitly represented in the model (Figure 7c). During seismic events, these joints may exhibit nonlinear behaviors such as opening, closing, and sliding. These interactions are simulated using a surface-based contact algorithm (Alembagheri and Ghaemian, 2016; Wang et al., 2017).

The normal behavior of joints is modeled through a contact formulation that permits no tensile strength, utilizing an exponential pressure-overclosure relationship described by Equation 11:

$p=\left\{\begin{array}{cc}0,\hfill & \text{when\hspace{0.17em}separation\hspace{0.17em}occurs\hspace{0.17em}}\left(d\le -c_0\right),\hfill \\ \frac{p_0}{e-1}\left[\left(1+\frac{d}{c_0}\right)\left(\exp \left(1+\frac{d}{c_0}\right)-1\right)\right],\hfill & \text{for\hspace{0.17em}contact\hspace{0.17em}}\left(d>-c_0\right).\hfill \end{array}\right.(11)$

In this formulation, $p$ represents the contact pressure (Pa), $d$ denotes the overclosure distance (m), $p_0$ is the reference contact pressure (Pa) at zero opening, and $c_0$ specifies the initial contact distance (m).

The normal contact interaction follows an exponential pressure-overclosure relationship where surfaces initiate pressure transmission when their separation distance in the contact normal direction decreases below the critical value $c_0$ (Figure 8a). The tangential contact behavior is governed by a penalty-based Coulomb friction formulation. Relative sliding initiates when the interfacial shear stress $\tau$ surpasses the critical threshold ${\tau }_{\text{crit}}=\mu p$, where $\mu$ denotes the friction coefficient. The formulation incorporates finite stick stiffness, permitting small elastic displacements during the sticking phase (Figure 8b), controlled by a specified slip tolerance. The analysis adopts $\mu =1.0$ and a slip tolerance of $1\times 10^{-4}$ to implicitly account for shear key effects at contraction joints. For normal contact behavior, the parameters $p_0=1$ MPa and $c_0=0.1$ mm were implemented (Nikkhakian et al., 2020).

Figure 8

Figure 8. Joint stress-displacement behavior, reproduced from Dassault Systemes (Systemes, 2022): (a) Normal direction, (b) Tangential direction.

The analyses consider a winter reservoir water level at EL. 268.21 m, corresponding to a depth of 94.48 m. Following established numerical wave propagation guidelines (Goldgruber, 2015), the reservoir domain extends to twice its depth to minimize artificial boundary reflections. This geometric configuration incorporates an upstream non-reflective boundary condition to properly absorb outgoing pressure waves. The interaction at the reservoir-foundation interface was modeled using a tie constraint, which rigidly connects the acoustic fluid nodes to the foundation solid elements. This approach creates a fully reflective boundary condition at the reservoir bottom. While some studies indicate that sediment absorption can influence the structural response, this effect is not included in the present analysis. This simplification is consistent with the findings of Løkke and Chopra (Løkke and Chopra, 2019), who demonstrated that for models which already include full dam-water-foundation interaction and the associated radiation damping, the additional energy dissipation from reservoir bottom sediments has a negligible influence on the overall system response. Therefore, following this approach, sediments and their associated wave reflection effects are not modeled. The reservoir water is modeled with 13,062 AC3D8 elements in three dimensions, contrasted with 224 AC2D4 elements in the two-dimensional representation.

The finite element model of the rock foundation (Figure 7a) consists of 13,420 C3D8R solid elements and 2,922 CIN3D8 infinite elements. The placement of absorbing boundaries was optimized via convergence analysis, requiring the rock foundation to extend laterally (cross-stream) and vertically to twice the dam height below its base. The final foundation dimensions are 264 m (stream), 264 m (cross-stream), and 120 m (depth).

The two-dimensional finite element model of the foundation is discretized into 290 CPE4R plane-strain solid elements and 49 CINPE4 infinite elements, as depicted in Figure 7b. Effective earthquake forces are computed using both the direct FE method and the analytical approach. Figure 9 demonstrates the 2D FE system, including a corresponding 1D corner column for auxiliary analyses, which facilitates the computation of effective earthquake forces via the direct FE method for the 3D foundation system.

Figure 9

Figure 9. FE model showing 2D system and companion 1D column for calculating effective earthquake forces through direct finite element method.

A critical practical consideration in dam seismic analysis is the computational cost associated with different modeling approaches. For the Pine Flat Dam case study, the complete nonlinear time-history analysis using the three-dimensional model–comprising a total of 42,787 elements (13,383 dam elements +13,420 foundation elements +2,922 infinite elements +13,062 reservoir elements)–required substantially greater computational resources than the two-dimensional model with 851 total elements (201 dam elements +290 foundation elements +49 infinite elements +224 reservoir elements). Based on the model complexity and established scaling relationships for nonlinear dynamic analysis, the 3D simulation demanded approximately 25–30 times more CPU time than its 2D counterpart. This significant computational expense is attributed to the increased degrees of freedom, contact nonlinearity at contraction joints, and smaller stable time increments required for accurate wave propagation in three dimensions.

4.2 Seismic excitation

The analysis employs recorded ground accelerations from the 1952 Kern County earthquake (Taft Lincoln School station) as seismic input (PEER, 2017). The direct finite element analysis of the coupled dam-reservoir-foundation system requires seismic excitation to be applied at the foundation base. While conventional approaches employ frequency-domain deconvolution techniques relying on one-dimensional shear wave propagation through layered media (Kramer, 1996; Hashash et al., 2009), this study adopts an innovative time-domain iterative methodology (Abbasiverki, 2023). The implemented technique models wave propagation through a one-dimensional finite element column, where P- and S-waves are systematically refined through successive iterations until the surface motion converges to the target acceleration record. The complete computational procedure is outlined in Algorithm 1.

Algorithm 1

Algorithm 1. Seismic input motion deconvolution using 1D FE column.

The time-domain deconvolution approach offers two key advantages: (i) consistency in material properties between the one-dimensional analysis column and the full three-dimensional finite element model, and (ii) compatibility with existing seismic analysis tools that assume vertically-propagating shear waves. In the present implementation, the deconvolution algorithm has been enhanced to process both compressional (P-wave) and shear (S-wave) seismic wave components. The deconvolution procedure is applied at a control point located on the abutment’s top surface (free ground surface of the rock foundation), considering all three orthogonal components (stream, vertical, and cross-stream) of the Taft earthquake record. Through an iterative refinement process, the solution converges when the discrepancy between the simulated free-surface motion and the target surface recording falls below 5%, as illustrated in Figure 10.

Figure 10

Figure 10. Time history of acceleration and the associated response spectrum for the Taft earthquake, showing measurements at the free surface (control point), the base (after deconvolution), and the column surface (after convolution).

5 Finite element analysis results

The computational methodology employs a sequential analysis framework. Initial eigenvalue analysis identifies the system’s fundamental dynamic properties, followed by a two-stage coupled analysis. The static phase establishes equilibrium under permanent loads (self-weight and hydrostatic pressure), while the subsequent dynamic phase performs nonlinear time-history simulation using an HHT-$\alpha$ algorithm with automatic time-step adjustment based on convergence monitoring.

5.1 Dynamic characterization and model validation

Validating the numerical model is fundamental to establishing confidence in the analysis results and ensuring the dam’s structural safety. This paper uses the steps recommended by Salamon and Hariri-Ardebili (Salamo et al., 2024) for verification and validation of the model.

An eigenvalue analysis was performed to evaluate the dynamic characteristics of both 3D and 2D numerical models, with validation against experimental tests conducted by Rea et al. (Rea et al., 1972). The measured natural frequencies were 3.47 Hz, 4.13 Hz, and 5.4 Hz for the first, second, and fourth bending modes, respectively. The 3D FE model predicted corresponding values of 3.78 Hz (+9%), 4.17 Hz (+1%), and 6.04 Hz (+12%), while the 2D FE model yielded 3.39 Hz (−3%), 4.01 Hz (−3%), and 5.23 Hz (−3%) for these modes. Figure 11 compares the crest mode shapes derived from the 3D FE model and experiments.

Figure 11

Figure 11. Crest mode shapes derived from the 3D FE model and experiments by Rea et al. (Rea et al., 1972).

The mode shapes of the 3D FE model are shown in Figure 12, demonstrating good agreement with experimental crest mode shapes (Figure 11), though some eigenvectors required sign adjustment (a common numerical convention due to the arbitrary $\pm$1 scaling of eigenvectors in eigenvalue analysis). The systematic frequency overestimation in the 3D model suggests greater structural stiffness than the actual dam, likely due to the omission of contraction joints in the linear eigenvalue analysis.

Figure 12

Figure 12. Mode shapes of the 3D FE model (left to right: modes one–4).

Notably, the third bending mode was not captured experimentally because the vibration generators were positioned on monolith no. 16 – coinciding with the inflection point of this mode. The 2D FE results (Figure 13) show greater deviation from experimental data, particularly for the third bending mode, which exhibits exaggerated curvature. This discrepancy highlights the 3D model’s superior ability to represent actual boundary conditions and structural constraints.

Figure 13

Figure 13. Vertical mode shape at centerline of block 16 derived from 3D FE model, 2D FE model and experiments by Rea et al. (Rea et al., 1972).

While the 2D model shows reasonable frequency predictions (within 3% of experimental values), its simplified geometry leads to less accurate mode shape representation. These findings collectively demonstrate that a 3D FE model is essential for proper investigation of the dam’s dynamic behavior, as it captures critical 3D effects that significantly influence both modal characteristics and structural response.

5.2 Comparison of the canyon response

Figure 14 presents peak ground acceleration (PGA) variations along the canyon’s free surface in three orthogonal directions, contrasting direct FE with analytical solutions for free-field foundation analysis. The direct FE method yields PGA ranges of 1.49–2.05 m/s² (stream), 1.48–2.27 m/s² (cross-stream), and 0.94–1.23 m/s² (vertical). In contrast, the analytical method produces corresponding ranges of 1.34–2.09 m/s², 0.97–1.69 m/s², and 0.92–1.32 m/s².

Figure 14

Figure 14. PGA variations across the canyon’s free surface computed through direct FE and analytical approaches for free-field foundation modeling.

The analysis reveals substantial discrepancies in cross-stream direction results, with the analytical method systematically producing lower ground motion estimates (0.97–1.69 m/s²) compared to the direct FE approach (1.48–2.27 m/s²). This underestimation occurs because the analytical method’s simplified 1D column representation fails to capture the complex wave propagation effects induced by the canyon geometry. In contrast, the direct FE method accurately accounts for these 3D topographic effects by modeling the complete free-field system, including the actual canyon configuration.

The 2D FE model was excited by the streamwise component of the Taft ground motion. Figure 15a shows the variations in PGA along the foundation surface for both direct FE and analytical free-field modeling approaches. The direct FE method yields a uniform PGA of 1.94 m/s² along the surface, consistent with the model’s uniform foundation geometry. In contrast, the analytical method overestimates PGA near side boundaries due to its inability to properly account for damping effects in free-field motions at these locations. This is confirmed by Figure 15b, which demonstrates that neglecting Rayleigh damping produces a uniform PGA distribution. Comparatively, the 3D FE model predicts a stream-direction PGA of 1.6 m/s² at monolith 16 (Figure 14), notably lower than the 2D result. This reduction reflects the significant influence of Pine Flat Dam canyon’s irregular topography on seismic wave propagation, which attenuates surface motions in the canyon center.

Figure 15

Figure 15. Comparison of PGA along the foundation free surface in 2D FE analysis: Direct FE and analytical approaches for free-field foundation modeling, evaluated with and without Rayleigh damping.

Figure 16 shows the stream-direction velocity contours at $t$ = 0.15 s, capturing the first notable divergence between the 2D flat and 3D canyon models. The 2D simulation exhibits classical layered wave propagation with coherent reflections consistent with Snell’s law, whereas the 3D model displays complex scattering patterns caused by canyon topography. The irregular geometry disrupts wavefront coherence through multi-directional reflections and diffraction, leading to localized velocity amplifications (“hotspots”) near canyon walls. This contrast highlights how topographic complexity fundamentally modifies seismic wave propagation, with implications for site-specific ground motion predictions. The complete temporal evolution ($t$ = 0.04–10 s) is provided in Supplementary Appendix A.

Figure 16

Figure 16. Comparison of stream-direction velocity contours in 2D and 3D models at $t$ = 0.15 s.

5.3 Comparison of the dam response

Figure 17a compares the joint sliding responses obtained from direct FE and analytical approaches, illustrating the spatial distribution and temporal behavior of the sliding displacements. The spatial variation along the dam crest in the upstream and downstream directions is presented, along with the time-history response at the critical 25th joint.

Figure 17

Figure 17. Comparison of joint sliding and opening responses between direct FE and analytical approaches for free-field foundation modeling: (a) Joint sliding, (b) Joint opening.

The direct FE analysis reveals a maximum downstream sliding displacement of approximately 10 mm, with this peak displacement occurring at the 25th joint. In contrast, the analytical method demonstrates consistent underestimation of sliding magnitudes, ranging from 20% to 50%, across the central crest region (−90 m to +90 m). This systematic underestimation pattern correlates well with the PGA discrepancies previously illustrated in Figure 14. The link between inaccurate input motion and non-conservative structural response is most critically demonstrated at the 25th joint, which experiences the maximum sliding displacement. While the analytical method yields a 5% lower stream-direction PGA at the foundation free-field, this error is dynamically amplified by the dam’s inertia and its interaction with the reservoir and foundation. Consequently, the underestimation of the peak sliding displacement at the dam crest is significantly larger, reaching 17% (downstream) and 28% (upstream). This amplification at the most critical location highlights the vital importance of accurate free-field input for predicting ultimate limit states.

Figure 17b illustrates notable differences in joint opening behavior when comparing direct FE and analytical approaches. The spatial distribution highlights the maximum opening at the 27th joint (proximal to the right abutment), and the corresponding time-history response at this location is also presented. The direct FE method yields a maximum opening of 69.90 mm at the 27th joint, contrasting sharply with the analytical method’s prediction of 43.80 mm - representing a significant 37.3% underestimation.

Across the dam structure, the analytical method systematically underestimates joint openings, with approximately 60% of joints showing 20%–50% lower values compared to direct FE results. This discrepancy is particularly pronounced in the cross-stream direction, where prediction errors exceed those observed for stream-direction responses. The pattern of underestimation correlates strongly with the canyon free-surface motion shown in Figure 14. This correlation indicates that the analytical approach fails to properly capture critical cross-stream dam-foundation interaction mechanisms.

The analytical method’s limitations in capturing joint opening behavior prove more severe than its shortcomings in predicting sliding responses. This performance gap highlights fundamental challenges in modeling complex 3D dam-foundation interactions using simplified analytical approaches, particularly for cross-stream motions that significantly influence joint opening behavior.

Figure 18 presents orthogonal acceleration response spectra components in the center of the crest and near the left abutment. The maximum acceleration responses obtained from the direct FE method reach 48.50 m/s² (stream) and 138.10 m/s² (cross-stream), as shown in Figure 18a. The corresponding values for the analytical method are 30% and 54% lower, respectively. It is worth noting that high-frequency accelerations with frequencies higher than 10 Hz are generated due to the impact load from the opening/closing of the joint. As seen in Figure 18b, such high-frequency accelerations do not exist near the left abutment, where joint opening is insignificant. Figure 18b also demonstrates that the analytical method yields maximum response spectra values (stream: $-10\%$; cross-stream: $-30\%$) lower than the direct FE results. This systematic underestimation is especially significant for cross-stream responses and at the dam center location, correlating with areas of increased joint opening.

Figure 18

Figure 18. Comparison of acceleration response spectra between direct FE and analytical approaches for free-field foundation modeling: (a) crest center, (b) Near the left abutment.

The envelope of maximum principal stresses over the upstream and downstream faces of the 3D dam model is presented in Figures 19a,c for direct FE and analytical approaches, respectively. The direct FE results (Figure 19a) reveal peak stresses of 1.8–4.6 MPa concentrated in monoliths 18–20, particularly at the heel and downstream regions. While these stresses remain below the threshold for significant nonlinear behavior in Pine Flat Dam under the Taft earthquake loading, the analytical method (Figure 19c) shows consistently lower stress magnitudes (1.8–3 MPa) localized in monoliths 16–17, demonstrating a systematic underestimation compared to the direct FE solution.

Figure 19

Figure 19. Envelope of maximum principal stresses in 3D and 2D free-field modeling of the rock foundation: (a) 3D Direct FE (left: downstream face, right: upstream face), (b) 2D Direct FE, (c) 3D Analytical (left: downstream face, right: upstream face), (d) 2D Analytical.

Figures 19b,c display the maximum envelope of principal stresses in the 2D model for the same free-field modeling techniques. The stress distributions from both methods are nearly identical, with high stress concentrations at the downstream side and heel–a pattern consistent with the 3D FE model results. However, the maximum stress in the 2D model ($\sim$1 MPa) is lower than in the 3D model. This discrepancy occurs because the simplified 2D analysis neglects several critical factors, including monolith interactions, contraction joints, cross-stream vibrations, and canyon shape effects. This fundamental difference originates from a reversal in the PGA trend between the foundation and the dam crest. While Figure 14 showed that 3D canyon effects cause wave scattering and attenuation at the foundation level (leading to a lower foundation PGA compared to the 2D model), this scattered wave energy is subsequently amplified by the dam’s own complex 3D inertial response and monolith interactions.

This significant amplification is unambiguously captured at the dam crest, as shown in Figure 20. The 3D FE model exhibits a peak ground acceleration (PGA) of 15.34 m/s²—86% higher than the 2D model’s PGA of 8.23 m/s². Furthermore, the response spectra reveal that for frequencies higher than 3.62 Hz, the 3D model produces consistently larger spectral accelerations. This demonstrates that the simplified 2D model fails to capture the critical dynamic phenomena that generate significantly higher inertial forces. These elevated forces at the crest are a primary driver of the more severe stress concentrations observed in Figure 19 and the pronounced nonlinear joint behavior detailed in Figure 17.

Figure 20

Figure 20. Comparison of stream-direction acceleration time-history and corresponding acceleration response spectra at the dam crest center, obtained from 2D and 3D finite element models.

These results highlight a crucial consideration for engineering practice: the choice between 2D and 3D modeling should be guided by the project’s design phase. While the computational efficiency of 2D analysis may be suitable for preliminary design stages, its inability to capture the amplified responses and complex interactions demonstrated here makes it inadequate for the final design of dams in seismically active regions or complex topographies. For such critical applications, the comprehensive insight provided by 3D modeling is indispensable for ensuring structural safety and reliability.

Figure 21 compares acceleration time histories and response spectra between methods in 2D FE model. The results show peak and response accelerations of 8.8 m/s² and 49.5 m/s², respectively, with negligible differences between the two approaches. While both methods produce comparable results in the 2D FE model, the analytical approach significantly underestimates dam responses in the 3D analysis. This discrepancy arises from fundamental differences in how each method handles foundation modeling: in the 2D case, both approaches use equivalent 1D free-field analysis to determine earthquake forces, while in the 3D approach, the direct FE method properly accounts for the complete free-field system, including irregular foundation topography. The analytical method, by contrast, employs simplified foundation assumptions that fail to capture these critical geometric complexities. Consequently, for reliable seismic assessment of concrete gravity dams, a comprehensive 3D FE model incorporating direct free-field foundation modeling becomes imperative.

Figure 21

Figure 21. Acceleration time history and corresponding response spectra at crest of the dam in 2D FE model for free-field foundation modeling using direct FE and analytical methods.

6 Conclusion

This study presented a comprehensive investigation of the seismic response of concrete gravity dams through comparative analysis using both 3D and 2D nonlinear finite element models. The research specifically evaluates two free-field foundation modeling approaches: (i) the direct FE method, which accurately represents canyon effects, and (ii) the analytical method, which employs simplified free-field assumptions. The analysis incorporates critical factors, including dam-reservoir-foundation interaction, nonlinear joint behavior, and cross-stream seismic excitation. Three key findings emerge from this investigation:

First, regarding free rock surface motion, the analytical method demonstrates significant limitations. In 3D FE analysis, the simplified 1D column representation fails to capture complex wave propagation effects induced by canyon geometry. For 2D FE models with uniform foundation surfaces, the analytical method conversely overestimates motions near non-reflecting boundaries due to inadequate damping representation.

Second, the foundation modeling approach substantially impacts seismic response predictions. For the Pine Flat Dam’s 3D FE analysis, the analytical method systematically underestimates critical responses, particularly at locations experiencing significant joint sliding and opening. This non-conservative bias underscores the importance of precise free-field modeling. In contrast, 2D analyses show comparable results between methods since both employ equivalent 1D free-field analysis.

Third, comparative analysis reveals fundamental limitations of 2D modeling approaches. While computationally efficient, 2D models cannot capture essential three-dimensional effects, including canyon-induced stress amplification, monolith interactions, and true free-field response in irregular topography. These findings suggest that 2D analysis may serve as a preliminary assessment, but the final design requires 3D modeling with a direct FE free-field representation. A key insight from this study is that 2D analysis is adequate only when the foundation topography is sufficiently uniform to be accurately represented by a one-dimensional wave-propagation model, as evidenced by the nearly identical dam responses from both free-field methods in the 2D simulation.

Based on these findings, it is advised that 2D analysis be confined to the preliminary design phase for sites with simple, planar foundation topography. For the final design, the results underscore the necessity of a 3D model. This is particularly crucial for dams situated in narrow, V-shaped, or U-shaped canyons, where the analysis confirms that pronounced curvature generates complex arching actions that a 2D model fundamentally cannot represent. For such topographically complex sites, a 3D analysis is concluded to be indispensable.

Consequently, the findings mandate a phased modeling strategy. The efficiency of 2D analysis suits preliminary design, while its failure to capture critical 3D mechanisms renders it inadequate for final design. For dams in high-seismic or topographically complex sites, the computational investment in a 3D model is not a choice but a necessity for ensuring safety.

Future studies should investigate how the divergence between analytical and direct FE solutions varies as a function of canyon geometry, canyon symmetry, foundation rock properties and their interaction with the dam, as well as the characteristics of seismic input–specifically intensity, duration, and frequency content. Additionally, foundation heterogeneity is a critical factor that limits the applicability and accuracy of both modeling approaches, and warrants further examination in terms of its influence on model correlation and predictive capability.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

RA: Conceptualization, Formal Analysis, Investigation, Methodology, Resources, Validation, Visualization, Writing – original draft. RK: Conceptualization, Validation, Writing – review and editing, Methodology, Resources. MH-A: Conceptualization, Validation, Writing – review and editing, Project administration, Visualization.

Funding

The author(s) declared that financial support was not received for this work and/or its publication.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author MH-A declared that they were an editorial board member of Frontiers at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbuil.2025.1737466/full#supplementary-material

Supplementary_Appendix_A | Temporal evolution of stream-direction velocity contours comparing 2D and 3D finite element models of the Pine Flat Dam foundation at time instances t = 0.04, 0.15, 0.26, 0.5, 1, 5, and 10 seconds.

Footnotes

¹It should be noted that this assumption is violated for nearby seismic sources or when wave propagation from the source to the site is inclined, in which case the procedure is not applicable (Mohammadnezhad et al., 2025).

References

Abbasiverki, R. (2023). Numerical modelling considerations for analysis of concrete hydraulic structures subjected to high-frequency seismic loads. Stockholm, Sweden: Kungliga Tekniska högskolan.

Google Scholar

Abbasiverki, R., Malm, R., Ansell, A., and Nordström, E. (2021). Nonlinear behaviour of concrete buttress dams under high-frequency excitations taking into account topographical amplifications. Shock Vib. 2021, 4944682. doi:10.1155/2021/4944682

CrossRef Full Text | Google Scholar

Akköse, M., and Şimşek, E. (2010). Non-linear seismic response of concrete gravity dams to near-fault ground motions including dam-water-sediment-foundation interaction. Appl. Math. Model. 34, 3685–3700. doi:10.1016/j.apm.2010.03.019

CrossRef Full Text | Google Scholar

Akpinar, U., Arici, Y., and Binici, B. (2025). Post-earthquake effects on the seismic performance of concrete gravity dams. Struct. Infrastructure Eng. 21, 10–23. doi:10.1080/15732479.2023.2180522

CrossRef Full Text | Google Scholar

Alembagheri, M., and Ghaemian, M. (2016). Seismic performance evaluation of a jointed arch dam. Struct. Infrastructure Eng. 12, 256–274. doi:10.1080/15732479.2015.1009124

CrossRef Full Text | Google Scholar

Alves, S., and Hall, J. (2006). System identification of a concrete arch dam and calibration of its finite element model. Earthq. Engineering and Structural Dynamics 35, 1321–1337. doi:10.1002/eqe.575

CrossRef Full Text | Google Scholar

Ameen, A. M. S., Ibrahim, Z., and Othman, F. (2017). Three-dimensional seismic response analysis for a rockfill dam. J. Comput. Theor. Nanosci. 14, 6003–6013. doi:10.1166/jctn.2017.7048

CrossRef Full Text | Google Scholar

Arici, Y., Binici, B., and Aldemir, A. (2014). Comparison of the expected damage patterns from two- and three-dimensional nonlinear dynamic analyses of a roller compacted concrete dam. Struct. Infrastructure Eng. 10, 305–315. doi:10.1080/15732479.2012.753921

CrossRef Full Text | Google Scholar

Azmi, M., and Paultre, P. (2002). Three-dimensional analysis of concrete dams including contraction joint non-linearity. Eng. Struct. 24, 757–771. doi:10.1016/s0141-0296(02)00005-6

CrossRef Full Text | Google Scholar

Basu, U., and Chopra, A. (2004). Perfectly matched layers for transient elastodynamics of unbounded domains. Int. J. Numer. Methods Eng. 59, 1039–1074. doi:10.1002/nme.896

CrossRef Full Text | Google Scholar

Bielak, J., Loukakis, K., Hisada, Y., and Yoshimura, C. (2003). Domain reduction method for three-dimensional earthquake modeling in localized regions, part I: theory. Bull. Seismological Soc. Am. 93, 817–824. doi:10.1785/0120010251

CrossRef Full Text | Google Scholar

Bybordiani, M., and Arıcı, Y. (2017). The use of 3d modeling for the prediction of the seismic demands on the gravity dams. Earthq. Eng. and Struct. Dyn. 46, 1769–1789. doi:10.1002/eqe.2880

CrossRef Full Text | Google Scholar

Chen, D., Du, C., Yuan, J., and Hong, Y. (2012). An investigation into the influence of damping on the earthquake response analysis of a high arch dam. J. Earthq. Eng. 16, 329–349. doi:10.1080/13632469.2011.638697

CrossRef Full Text | Google Scholar

Chen, D., Hou, C., and Wang, F. (2021). Influences on the seismic response of the gravity dam-foundation-reservoir system with different boundary and input models. Shock Vib. 2021, 6660145. doi:10.1155/2021/6660145

CrossRef Full Text | Google Scholar

Chopra, A. (2012). Earthquake analysis of arch dams: factors to be considered. J. Struct. Eng. 138, 205–214. doi:10.1061/(asce)st.1943-541x.0000431

CrossRef Full Text | Google Scholar

Chuhan, Z., Jianwen, P., and Jinting, W. (2009). Influence of seismic input mechanisms and radiation damping on arch dam response. Soil Dyn. Earthq. Eng. 29, 1282–1293. doi:10.1016/j.soildyn.2009.03.003

CrossRef Full Text | Google Scholar

El-Aidi, B., and Hall, J. F. (1989). Non-linear earthquake response of concrete gravity dams part 1: modelling. Earthq. Engineering and Structural Dynamics 18, 837–851. doi:10.1002/eqe.4290180607

CrossRef Full Text | Google Scholar

Enzell, J., Malm, R., Abbasiverki, R., and Ahmed, L. (2021). “Non-linear behavior of a concrete gravity dam during seismic excitation: a case study of the pine flat dam,” in Numerical analysis of dams: proceedings of the 15th ICOLD International Benchmark workshop 15 (Springer), 99–112.

Google Scholar

Faggiani, G., and Masarati, P. (2018). “Seismic behaviour of monticello dam: a comparison between numerical results and experimental measures,” in Validation of dynamic analyses of dams and their equipment (Boca Raton, FL: CRC Press), 535–551.

CrossRef Full Text | Google Scholar

Faggiani, G., Masarati, P., and Frigerio, A. (2021). “Assessment of the dynamic response of pine flat concrete gravity dam. fem simulation of dam-foundation interaction,” in Numerical analysis of dams: proceedings of the 15th ICOLD International Benchmark workshop 15 (Springer), 113–127.

Google Scholar

Fry, J.-J., and Matsumoto, N. (2018). “Validation of dynamic analyses of dams and their equipment,” in Contributions to the international symposium on the qualification of dynamic analyses of dams and their equipments, 31 August-2 September 2016, saint-malo. France: CRC Press.

Google Scholar

Goldgruber, M. (2015). Nonlinear seismic modelling of concrete dams, Ph.D. thesis. Graz, Austria: Technical University of Graz.

Google Scholar

Hall, J. (2006). Problems encountered from the use (or misuse) of rayleigh damping. Earthq. Eng. Struct. Dyn. 35, 525–545. doi:10.1002/eqe.541

CrossRef Full Text | Google Scholar

Hariri-Ardebili, M. A. (2014). Impact of foundation nonlinearity on the crack propagation of high concrete dams. Soil Mech. Found. Eng. 51, 72–82. doi:10.1007/s11204-014-9257-9

CrossRef Full Text | Google Scholar

Hariri-Ardebili, M. A. (2018). Risk, reliability, resilience (R3) and beyond in dam engineering: a state-of-the-art review. Int. Journal Disaster Risk Reduction 31, 806–831. doi:10.1016/j.ijdrr.2018.07.024

CrossRef Full Text | Google Scholar

Hariri-Ardebili, M. A. (2024). Quantifying modeling uncertainties in seismic analysis of dams: insights from an international benchmark study. Earthq. Eng. and Struct. Dyn. doi:10.1002/eqe.4064

CrossRef Full Text | Google Scholar

Hashash, Y., Groholski, D., Phillips, C., and Park, D. (2009). Deepsoil v5, user manual and tutorial.

Google Scholar

He, J.-T., Ma, H.-F., Zhang, B.-Y., and Chen, H.-Q. (2010). Method and realization of seismic motion input of viscous-spring boundary, Shuili Xuebao. J. Hydraulic Eng. 41, 960–969.

Google Scholar

Jin, A.-Y., Pan, J.-W., Wang, J.-T., and Zhang, C. (2019). Effect of foundation models on seismic response of arch dams. Eng. Struct. 188, 578–590. doi:10.1016/j.engstruct.2019.03.048

CrossRef Full Text | Google Scholar

Kramer, S. (1996). Geotechnical earthquake engineering, prentice-hall civil engineering and engineering mechanics series. Upper Saddle River, NJ, USA: Prentice Hall.

Google Scholar

Løkke, A., and Chopra, A. K. (2018). Direct finite element method for nonlinear earthquake analysis of 3-dimensional semi-unbounded dam–water–foundation rock systems. Earthq. Eng. and Struct. Dyn. 47, 1309–1328. doi:10.1002/eqe.3019

CrossRef Full Text | Google Scholar

Løkke, A., and Chopra, A. K. (2019). Direct finite element method for nonlinear earthquake analysis of concrete dams: simplification, modelling, and practical application. Earthq. Eng. and Struct. Dyn. 48, 818–842. doi:10.1002/eqe.3150

CrossRef Full Text | Google Scholar

Lysmer, J., and Kuhlemeyer, R. (1969). Finite dynamic model for infinite media. J. Eng. Mech. Div. 95, 859–878. doi:10.1061/jmcea3.0001144

CrossRef Full Text | Google Scholar

Mandal, K. K., and Maity, D. (2018). Transient response of concrete gravity dam considering dam-reservoir-foundation interaction. J. Earthq. Eng. 22, 211–233. doi:10.1080/13632469.2016.1217804

CrossRef Full Text | Google Scholar

Mazzieri, I., Muhr, M., Stupazzini, M., and Wohlmuth, B. (2022). Elasto-acoustic modeling and simulation for the seismic response of structures: the case of the tahtalı dam in the 2020 i^.zmir earthquake. J. Comput. Phys. 466, 111411. doi:10.1016/j.jcp.2022.111411

CrossRef Full Text | Google Scholar

Mohammadnezhad, H., Zafarani, H., and Ghaemian, M. (2019a). Domain reduction method for seismic analysis of dam-foundation-fault system. Sci. Iran. 26, 145–156. doi:10.24200/sci.2018.20696

CrossRef Full Text | Google Scholar

Mohammadnezhad, H., Ghaemian, M., and Noorzad, A. (2019b). Seismic analysis of dam-foundation-reservoir system including the effects of foundation mass and radiation damping. Earthq. Eng. Eng. Vib. 18, 203–218. doi:10.1007/s11803-019-0499-4

CrossRef Full Text | Google Scholar

Mohammadnezhad, H., Bidaki, S. Z., Hariri-Ardebili, M. A., and Noorbakhsh-Saleh, M. (2025). Seismic response of geostructures to obliquely incident pulse-type near-fault p and sv waves. Eng. Struct. 330, 119917. doi:10.1016/j.engstruct.2025.119917

CrossRef Full Text | Google Scholar

Nielsen, A. H. (2006). “Absorbing boundary conditions for seismic analysis in abaqus,” in ABAQUS users’ conference, 359–376.

Google Scholar

Nielsen, A. H. (2014). Towards a complete framework for seismic analysis in abaqus. Proc. Institution Civil Engineers-Engineering Computational Mechanics 167, 3–12. doi:10.1680/eacm.12.00004

CrossRef Full Text | Google Scholar

Nikkhakian, B., Alembagheri, M., and Shayan, R. (2020). Parametric investigation of canyon shape effects on the seismic response of 3d concrete gravity dam model. Geotechnical Geol. Eng. 38, 6755–6771. doi:10.1007/s10706-020-01467-3

CrossRef Full Text | Google Scholar

Nowak, P., and Hall, J. (1990). Arch dam response to nonuniform seismic input. J. Engineering Mechanics 116, 125–139. doi:10.1061/(asce)0733-9399(1990)116:1(125)

CrossRef Full Text | Google Scholar

PEER (2017). Ground motion database. Available online at: http://ngawest2.berkeley.edu/.

Google Scholar

Poul, M. K., and Zerva, A. (2018). Nonlinear dynamic response of concrete gravity dams considering the deconvolution process. Soil Dyn. Earthq. Eng. 109, 324–338. doi:10.1016/j.soildyn.2018.03.025

CrossRef Full Text | Google Scholar

Poul, M. K., and Zerva, A. (2019). Comparative evaluation of foundation modeling for ssi analyses using two different abc approaches: applications to dams. Eng. Struct. 200, 109725. doi:10.1016/j.engstruct.2019.109725

CrossRef Full Text | Google Scholar

Rea, D., Liaw, C., and Chopra, A. K. (1972). Dynamic properties of pine flat dam. Berkeley, CA: Earthquake Engineering Research Center, University of California.

Google Scholar

Salamon, J. W., and Hariri-Ardebili, M. A. (2024). Verification, validation, and uncertainty quantification (vvuq) in structural analysis of concrete dams. Front. Built Environ. 10, 1452415. doi:10.3389/fbuil.2024.1452415

CrossRef Full Text | Google Scholar

Salamon, J. (2011). Seismic induced loads on spillway gates, phase I - literature review, technical report DSO-11-06. Denver, CO: U.S. Bureau of Reclamation, Denver, Colorado.

Google Scholar

Salamon, J. (2018). Evaluation of numerical models and input parameters in the analysis of concrete dams, technical report DSO-19-13. Denver, CO: U.S. Bureau of Reclamation, Denver, Colorado.

Google Scholar

Salamon, J., Wood, C., Hariri-Ardebili, M. A., Malm, R., and Faggiani, G. (2019). “Seismic analysis of pine flat concrete dam: formulation and synthesis of results,” in ICOLD international benchmark workshop on numerical analysis of dams, ICOLD-BW (Springer), 3–97.

Google Scholar

Song, Z., Wang, F., Liu, Y., and Su, C. (2018). Infinite element static-dynamic unified artificial boundary. Shock Vib. 2018, 7828267. doi:10.1155/2018/7828267

CrossRef Full Text | Google Scholar

Systémes, D. (2022). Abaqus 2022 theory manual, dassault systémes Simulia Corp. Providence, RI, USA.

Google Scholar

Wang, H., Feng, M., and Yang, H. (2012). Seismic nonlinear analyses of a concrete gravity dam with 3d full dam model. Bull. Earthq. Eng. 10, 1959–1977. doi:10.1007/s10518-012-9377-4

CrossRef Full Text | Google Scholar

Wang, G., Wang, Y., Lu, W., Yu, M., and Wang, C. (2017). Deterministic 3d seismic damage analysis of guandi concrete gravity dam: a case study. Eng. Struct. 148, 263–276. doi:10.1016/j.engstruct.2017.06.060

CrossRef Full Text | Google Scholar

Wang, G., Wang, Y., Lu, W., Yan, P., and Chen, M. (2020). Earthquake direction effects on seismic performance of concrete gravity dams to mainshock–aftershock sequences. J. Earthq. Eng. 24, 1134–1155. doi:10.1080/13632469.2018.1453423

CrossRef Full Text | Google Scholar

Wolf, J. P. (1989). Soil-structure-interaction analysis in time domain. Nucl. Engineering Design 111, 381–393. doi:10.1016/0029-5493(89)90249-5

CrossRef Full Text | Google Scholar

Zhang, W., Seylabi, E. E., and Taciroglu, E. (2019). An abaqus toolbox for soil-structure interaction analysis. Comput. Geotechnics 114, 103143. doi:10.1016/j.compgeo.2019.103143

CrossRef Full Text | Google Scholar

Zhao, C., and Valliappan, S. (1993). Seismic wave scattering effects under different canyon topographic and geological conditions. Soil Dyn. Earthq. Eng. 12, 129–143. doi:10.1016/0267-7261(93)90040-x

CrossRef Full Text | Google Scholar