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: f 0 = x 0 + c f r ̇ 0 (2) where x 0 represents boundary tractions and c f r ̇ 0 accounts for damping, with c f containing viscous coefficients and r ̇ 0 being nodal velocities from free-field analysis. The base boundary formulation simplifies to Equation 3: f 0 = 2 c f r ̇ I 0 (3) where 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 r ̇ 0 and boundary tractions 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.

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: f 0 = k f r 0 + c f r ̇ 0 + σ 0 n (4) with k f as system stiffness, r 0 = [ u v w ] T and r ̇ 0 = [ u ̇ v ̇ w ̇ ] T as boundary displacement/velocity vectors, σ 0 as free-field stress tensor, and 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 ( k f = 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:

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 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.

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 = 0 , when separation occurs  d ≤ − c 0 , p 0 e − 1 1 + d c 0 exp 1 + d c 0 − 1 , for contact  d > − c 0 . (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 τ surpasses the critical threshold τ crit = μ p , where μ 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 μ = 1.0 and a slip tolerance of 1 × 1 0 − 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).

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.

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.

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

Seismic input motion deconvolution using 1D FE column.

1:  Establish 1D finite element 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.

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.

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 ± 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.

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.

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.

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².

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 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 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.

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.

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.

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 ( ∼ 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.

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.