Dynamic Numerical Simulations of Dry-Stone Retaining Walls: Identification of the Seismic Behaviour Factor

Abstract

Dry-stone retaining walls can be found worldwide and constitute critical assets of the built heritage for many sloped territories, holding cultural and economic value. Their design currently follows empirical rules, though the first steps towards a static safety assessment have recently been proposed in the scientific and engineering literature. However, the seismic design of these structures still lacks research studies. Therefore, this work conducts discrete element simulations to assess their dynamic behaviour. First, the approach is validated through existing scaled-down shaking table experiments, and it is found that the numerical simulations are conservative (i.e., on the safe side). Next, full-scale dry-stone retaining walls are subjected to harmonic excitations as an idealisation of earthquakes. Finally, based on a simplified limit-equilibrium analytical tool, their seismic behaviour factor is estimated for the first time in the literature, which falls within the proposed values of the European standards (Eurocode 8). This will allow engineers to adopt a validated behaviour factor in practice to assess and design dry-stone retaining walls with a pseudo-static approach.

1. Introduction

Masonry is one of the oldest building materials worldwide and is currently the material adopted in a significant part of the built heritage. Many different building techniques can be found within this general typology. Among them, dry-stone masonry consists of a clever assemblage of rubble (undressed) stones without mortar. This vernacular technique can be found in many countries, including several UNESCO sites (Lavaux’s Terraces, Machu Picchu, Douro’s Valley, the Great Wall of Zimbabwe). More recently, the dry-stone technique and associated know-how have been listed as a World Intangible Heritage by UNESCO, recognising the specificities of this still living technique and its contribution to high-quality landscapes and built heritage.

Dry-stone retaining walls (DSRWs), which constitute a high portion of this built heritage, are still used for agricultural or transportation purposes. One can also highlight their advantages, including their low embodied energy [1], permeability to water flow [2,3] and adaptability to different environmental sites. However, appropriate maintenance actions for the dry-stone heritage are missing due to insufficient knowledge concerning this building technique and, until recently, a lack of design standards [4,5].

This state of the art led many researchers to study the static behaviour of DSRWs, with different goals. While some studies aimed to reproduce the particular shapes, also called bulging phenomena, observed in the field [6,7,8,9], others carried out parametric analyses via the Discrete Element Method (DEM) [10,11]. Finally, some research groups produced useful analytical tools to assess the stability of existing structures [12,13,14]. The literature encompasses experimental [9,14,15,16,17], analytical [12,13,18] and numerical DEM [19,20,21,22,23] works, conducted either in 2D [6,10,14,24,25] or 3D [9,17,21,26]. Therefore, the design of DSRWs under static loading is now relatively easy to process, given that dedicated safety rules and criteria are provided [4,5,27]. If prejudice prevents broader use of this building technique, there is, however, still a lack of knowledge regarding the seismic assessment of DSRWs.

Contrarily to static designs, the seismic design of such assets has been investigated only recently, with experimental, analytical and numerical pseudo-static studies [18,22,27]. Still, most of the building standards worldwide take seismic hazards into account in the design procedures of structures [28]. The most complete (but also complex) approach used to assess the seismic capacity of structures is the performance-based method. Originally introduced by Newmark [29], it is still in use nowadays [30,31]. It computes the displacement of a retaining wall during an earthquake motion. It then compares the response (e.g., top displacement) with the maximum allowable displacement, from a structural point of view, to draw conclusions on the acceptability of the design. This method often requires several simulations with different seismic inputs on the same structure [30] to provide statistically significant results. Thus, a simplified (and more engineering-oriented) approach is also proposed by the Eurocode standards (EC8) for the seismic design of retaining walls [32], where the seismic assessment of DSRWs is conducted using a pseudo-static method. A seismic behaviour factor (also known as the seismic response modification factor), denoted as r, reduces the reference design acceleration to account for a larger dynamic capacity, related to a reserve of ductility given by the non-linear behaviour of the system. However, such a factor has not been explicitly evaluated for DSRWs. As the structural behaviour of DSRWs (the presence of weak block interfaces) is different from that of classical concrete retaining walls, a dedicated seismic behaviour factor r for DSRWs should be identified. Therefore, as its contribution to the current state of the art, the present work proposes a research path to characterise a relevant seismic behaviour factor r for DSRWs. It involved scaled-down dynamic experiments on DSRW mock-ups [33] to validate a numerical DEM model before simulating full-scale DSRWs. Herein, Discrete Element Modelling (DEM) is preferred to Finite Element Modelling, as large displacements are assumed to occur at interfaces between blocks [34]. Additionally, DEM has already shown high potential to reproduce the brittle failure of such structures [10,11,19,20].

The study first describes the numerical DEM modelling strategy (Section 2). Next, the DEM model is validated against scaled-down laboratory experiments (Section 3). Section 4 presents numerical simulations of the full-scale DSRWs. Their comparison with pseudo-static predictions gives an initial estimate of the behaviour factor r. The conclusions are drawn in the final section.

2. Numerical DEM Strategy for Dynamic Simulations

The numerical strategy uses the plane strain software UDEC (version 7.0) developed by the Itasca consulting group [35]. The numerical model consists of an assembly of masonry units retaining a homogeneous backfill. Each masonry unit is represented by a deformable homogeneous medium characterised by linear elastic isotropic parameters (Young’s modulus E and Poisson’s ratio ν). The backfill is modelled as an elastoplastic medium: its mechanical behaviour is driven by elastic parameters (Young’s modulus E and Poisson’s ratio ν) and a Mohr–Coulomb plasticity criterion (friction angle φ, cohesion C, tensile strength T and dilation angle ψ). Both the units and backfill are therefore meshed using finite elements of size Δz, which should respect the following inequality to ensure an appropriate propagation of elastic waves through the media:

Δz < 0.1 × λ_min = 0.1 × V_S × T_min = 0.1 × [G/ρ]^0.5 × T_min

(1)

where λ_min, V_S and T_min are the propagated waves’ smallest wavelength, shear velocity and minimum period, respectively. G and ρ are the medium’s shear modulus and unit weight, respectively. Equation (1) does not constrain Δz much (see Section 3 and Section 4 for a formal check). Therefore, herein, Δz is defined based on the masonry unit height so that at least six subcontact points are generated between each unit and the backfill, as this parameter appears to be critical to reasonably reproduce the actual force distribution at the interfaces [20,22,36,37].

Zero-thickness contact interfaces are assumed between the units and between the unit and the backfill. They follow an elastoplastic Mohr–Coulomb criterion.

In UDEC, equations of motion are solved explicitly via a central finite difference scheme. For each computational node of the finite element mesh (both units and backfill), and at each time t, the dynamic equations of motions are solved as follows, thus giving new positions and velocities at time t + Δt:

m_ia_i(t) = ΣF_i(t)
v_i(t + Δt/2) = v_i(t − Δt/2) + a_i(t) × Δt
u_i(t + Δt) = u_i(t) + v_i(t + Δt/2) × Δt

(2)

where Δt corresponds to the timestep of the simulations; F_i corresponds to all the applied loads, namely the external loads (gravity, dynamic shaking at the bottom of the model), contact interactions between the different deformable bodies, and internal loads given by classical continuum mechanics calculations [35]; m_i, a_i, v_i and u_i correspond to the associated mass, acceleration, velocity and displacement of the computational node, respectively. At the new simulation time t + Δt, the positions of each computational node give the strain of each finite element, which, in turn, updates the stress as a result of the constitutive law of the material. Contact forces are also recomputed at this stage according to the following equations:

ΔF_n = A_c × k_n × Δu_n = A_c × k_n × [u_n(t + Δt) − u_n(t)]
ΔF_s = A_c × k_s × Δu_s = A_c × ks × [u_s(t + Δt) − u_s(t)]

(3)

where ΔF_n and ΔF_s are the normal and shear contact force increments at a given subcontact, A_c corresponds to the subcontact area, k_n (respectively, k_s) is the normal (respectively shear) elastic contact stiffness, and u_n (respectively u_s) represents the normal (respectively shear) relative displacement between the two contact points. The constitutive law of the interface is then applied (the Mohr–Coulomb model in this case), which reads:

F_n(t + Δt) = F_n(t) + ΔF_n, F_n < T × A_c
F_s(t + Δt) = F_s(t) + ΔF_s, |F_s| < F_{s, max} = C × A_c − F_n × tan(φ)

(4)

Additional expressions exist to include the dilation effect but are not presented because hereafter, the dilation angle is always taken to be equal to zero [35]. Updating all the forces in the system allows the resolution of the equations of motion (Equation (2)) for the next simulation time t + Δt.

Finally, numerical damping is used to model the actual energy dissipation. In dynamics, Rayleigh damping is generally preferred [37,38], including mass- and stiffness-proportional contributions, because it gives constant damping across a relatively wide frequency range. However, only mass-proportional damping is considered here, leading to a reasonable timestep Δt [38]. Additionally, the input signals are all harmonic. Therefore, the under-evaluation of damping at higher frequencies when disregarding stiffness-proportional damping is insignificant for the current simulations.

In Section 3 and Section 4, the numerical models are characterised by three different boundaries: free, fixed-velocity and free-field boundaries (Figure 1). In particular, the velocity of the bottom computational nodes (the backfill and bottom wall units) are fixed to match the dynamic solicitation. Next, a free-field boundary is applied to the right part of the model (vertical boundary). An elastic 1D numerical model is first run parallel to the primary model. It has the same dynamic input signal, elastic properties and computational nodes as the backfill boundary. At each time t, the velocities produced from this 1D model are directly applied to the vertical boundary of the main simulation. In other words, it assumes that the free-field boundary is far enough from the masonry wall (i.e., no impact of reflected waves on the wall in the media close to the boundary).

Spurious large deformations are obtained for non-cohesive backfill free layers because there is zero confining pressure, and hence zero shear resistance, of the elements involved. Therefore, an artificial thin cohesive layer is added to the backfill to stabilise the behaviour of the shallow elements (Figure 1). Cohesion of only 2 kPa with a tiny depth avoids unrealistic deformation without significantly modifying the global response of the system.

The following two sections focus on the numerical modelling of the scaled-down experiments (Section 3) and full-scale walls (Section 4). The specifics of each model and their numerical parameters are described in the corresponding sections.

3. Validation of the Model on Scaled-Down Experiments

3.1. Experimental Results

The numerical strategy is validated against an existing experimental campaign [33]. The reader is referred to the original study for a detailed description, while the primary outcomes are gathered herein. Several scaled-down dry-joint retaining walls with different slenderness H/B values (constant width B and different heights H) were subjected to harmonic horizontal shaking (Figure 2). The walls were 400 mm long, 34 mm wide and between 50 mm and 140 mm high. The bottom course of the masonry units was fixed to the container. The walls were retaining a sandy backfill of the same height and length, with a depth of 400 mm. An acceleration input signal (S(t)) of increasing amplitude and constant frequency f (equal to 10 Hz, 30 Hz, 50 Hz or 70 Hz) was applied through a shaking table to the bottom of the mock-up, as follows:

S(t) = A × t × sin(ω × t)

(5)

where A is a constant that describes the amplitude rate, and ω is the angular frequency, defined as: ω = 2 × π × f. The primary campaign used a frequency f of 10 Hz. The horizontal displacement d_top of the second-to-last (top) course of bricks was monitored through a laser sensor (Figure 2). Section 3.3 uses d_top to compare the experimental and numerical outcomes. In reference [33], wall failure was declared when d_top exceeded 5% of the wall height H. Note that H refers to the wall height free to move throughout the shaking, disregarding the first fixed masonry course (see Figure 2 and Figure 3).

3.2. Numerical Parameters of the UDEC Simulations

The numerical models’ geometric data (external shape and unit dimensions) is drawn from the experiments (Figure 1).

Table 1. Properties used for the numerical simulations of the scaled-down experimental tests.

	Units	Backfill	Wall–Backfill Interface
Volumetric properties
Unit weight γ (kN/m3) [33]	14.4	15.4	-
Elastic modulus E (MPa)	1000	10	-
Poisson ratio ν	0.2	0.3	-
Friction angle φ (°) [33]	-	45	-
Interface properties
Friction angle φ (°) [33]	32	-	32
Normal stiffness kn (Pa/m3)	2.6 × 1012	-	3.0 × 1010
Tangential stiffness ks (Pa/m3)	2.6 × 1012	-	3.0 × 1010
Numerical parameters
Timestep Δt (s)		≈1 × 10−6
Rayleigh viscous damping ξ		10% at 10Hz

gathers all the parameters considered in the simulations. The material properties considered in the numerical model are either collected from the original study [33] or from numerical studies that modelled the same material [21,22] when no experimental data were available. Apart from the small layer on the top of the backfill, the cohesion C (continuum and interface) is always equal to zero (dry sand backfill and dry contact wall assembly). Note that the timestep Δt is automatically computed by UDEC once all the numerical parameters are defined [35]. This led to 12 h of computation (on an HP computer with a speed of 2.3 GHz) for 1 s of simulation, while the experimental signals range from 15 to 30 s. According to Equation (1), the maximum allowable mesh size is 70 mm in the backfill and 800 mm in the wall, both being much larger than the actual numerical mesh size Δz, which is equal to 2 mm.

The viscous damping coefficient ξ used in the numerical simulations is based on the energy dissipation of the backfill. Preliminary simulations have allowed an evaluation of the cyclic strain of the backfill. By using computational nodes intentionally placed far from the wall (and the associated failure zones), a cyclic strain of 0.025% is found, which is the typical threshold where damping greater than 5% can be observed [39]. As two preliminary simulations using a damping ξ of 5% and 15% gave identical results (a relative difference less than 2%), a damping value ξ of 5% is considered thereafter. Since only the mass-proportional part of the Rayleigh damping is used, a multiplicative factor of 2 is applied to the chosen value to ensure appropriate effective damping (5%) at the critical frequency (Table 1) [35].

In the experiments, an absorbing foam placed at the back of the backfill was intended to reproduce a “free-field boundary condition”. Preliminary simulations directly modelling the foam with its properties (taken from [40]) revealed that modelling the right vertical boundary with either an absorbing foam or with a free-field boundary led to almost identical behaviour of the mock-up. The latter modelling strategy, being faster, has been considered in the rest of this section.

UDEC assumes a unique cross-section along the wall while the actual walls are 3D assemblages with four different cross-sections (Figure 3); the first three are equally distributed, while the fourth corresponds to the two ends of the wall. As the position of the headers (and, in general, the assemblage) in the height of the wall significantly influences the behaviour of dry-jointed structures [22,41,42], the four different cross-sections are expected to have different responses (Figure 3).

Figure 4 shows the horizontal displacement d_top of the second-to-last course obtained from the numerical simulations of the wall with a slenderness H/B = 2.7 using the first three cross-sections; the fourth one has the same response as the second one (Figure 3). Next, the global curve is computed, weighting the individual responses according to the proportions in Figure 3 [22]. For each d_top, one finds the corresponding times in each response and then weights them to obtain the global response time. This is illustrated in

Table 2. Computation of the weighted response illustrated with the threshold d_top equal to 5 mm.

Section No.	Proportion	Time for dtop = 5 mm (s)
1	7/23	18.4
2	7/23	21.5
3	7/23	19.9
4	2/23	21.5
Weighted average		20.1

for a single displacement d_top equal to 5 mm. As can be noticed in Figure 4, the weighted response is very close to the response of Cross-section 3 (a difference of less than 0.25 s for all points in the signal). Therefore, Cross-section 3 is assumed to represent the global response and is used alone in the following numerical simulations to save computational time.

3.3. Validation of the Modelling Strategy

As the first step, Figure 5 compares the numerical response d_top (only for Cross-section 3) with the experimental curves from five repeated tests on a wall with a slenderness H/B = 2.3. The obtained trends are globally similar (also found for all the studied configurations; not shown for brevity). In addition, the numerical response is encompassed by the experimental ones, highlighting the validity of the numerical strategy.

Wall (numerical or experimental) failure was declared when d_top exceeded 5% of the wall height H [33]. The acceleration at the instant of failure, taken from the dynamic signal, is denoted A_fail. Figure 6 and

Table 3. Experimental and numerical acceleration at failure (A_fail) for different wall slenderness: the mean value is given when several experimental points are available. The relative difference is computed by taking the experimental value as a reference. Finally, the detected failure modes are also indicated.

Wall Slenderness (H/B)	Exp. Afail (m/s2)	Exp. Failure Mode	Num. Afail (m/s2)	Num. Failure Mode	Error (%)
H/B = 1.3	4.2	S	3.7	S	−11%
H/B = 1.7	4.6	S	3.4	S	−27%
H/B = 2.3	3.3	S/T	3.2	S	−5%
H/B = 2.7	3.2	T	2.6	T	−18%
H/B = 3	2.7	T	1.7	T	−35%
H/B = 3.3	2.0	T	1.4	T	−32%
				Mean	−21%

compare the outcomes, both in terms of acceleration at failure A_fail and the observed failure mode, between the experiments and simulations.

Firstly, it can be observed that the numerical simulations accurately reproduce the observed failure mechanisms. Secondly, the numerical simulations underestimate, on average, the actual experimental responses by 20%. However, since most numerical parameters have been taken by default, the agreement is considered to be good. Furthermore, in [33], the authors noted that the friction angle of the backfill sand could be as high as 55° (larger than the considered value of 45°).

Table 4. Effect of the backfill friction angle on the numerical acceleration at failure A_fail of the wall with a slenderness H/B = 2.7.

Backfill Friction Angle φ (°)	45°	50°	55°
Error compared with experiment (%)	−18%	−13%	−10%

shows the numerical underestimations of A_fail for the wall with a slenderness H/B = 2.7 for three different backfill friction angles (45°, 50° and 55°), which significantly affect the capacity of the walls.

In addition, in Figure 4 and Figure 5, one can also note that the numerical responses are softer than the experimental ones, indicating that the estimated material properties may need correction (e.g., units of Young’s modulus or joint stiffness). The comparison of the accelerations at failure becomes sensitive to the chosen threshold value. For example, considering a threshold of 10% (instead of the original 5%) to define the experimental and numerical failures (Figure 5) decreases the average error from 21% to 16%. Finally, UDEC does not catch 3D effects that are also partly responsible for the observed discrepancies. In particular, the most resistant cross-section (Figure 3) may play a more significant role than the weaker cross-sections, contrary to the assumptions of Table 2.

3.4. Effect of the Frequency

Figure 7 compares the numerical and experimental responses obtained for an input signal with f = 30 Hz (instead of 10 Hz). An excellent agreement (a difference of less than 1%) is also noticeable when computing the acceleration at failure A_fail (a threshold of 5%), similar to the case with f = 10 Hz (Table 3).

To conclude, the precise calibration of the numerical model to fit the experimental results is not sought here, since the number of parameters is significant. Nevertheless, it has been shown that the numerical strategy is accurate enough for modelling the influence of both the wall’s slenderness and the signal frequency on the wall’s response. Therefore, this strategy is used further to simulate larger-scale systems (Section 4). Finally, the modelling approach underestimates the actual capacity, which is conservative and particularly important, as design recommendations are proposed in the following section.

4. Full-Scale Numerical Simulations: Extraction of the Seismic Behaviour Factor

4.1. Description of the Case Study

The full-scale DSRWs considered in this section stem from the “C3s” DSRW taken from the literature [16] (Figure 8).

Table 5. Geometric parameters of the case study wall “C3s” and the DSRWs numerically simulated in Section 4 [16]. The reader is referred to Figure 9 to understand the geometries better.

Geometrical Parameters	Original Wall “C3s”	DSRW1	DSRW2	DSRW3	DSRW4
Height H (m)	2.5	2.5	2.5	2.5	2.5
Width at base B (m)	0.7	0.7	0.7	0.7	0.7
Width at top b (m)	0.55	0.55	0.55	0.55	0.7
Internal batter λm (%)	0	0	0	0	0
External batter λv (%)	6	6	6	6	0
Bed inclination α (º)	9.1	9.1	9.1	3.4	0
Backfill slope β (º)	32.6	0	32.6	0	0

gathers all the geometrical parameters of the DSRWs numerically tested throughout this section (Figure 9). According to the experiments, the masonry units are 0.1 m high and 0.35 m long [16];

Table 6. Mechanical parameters of the full-scale DSRW simulations. Most of the data come from the reference experimental study [16,24,43], while the rest is taken on a default basis.

Mechanical Parameters	Units	Backfill [43]	Interface
Volumetric properties
Unit weight (kN/m3) [43]	20.0	14.9	-
Elastic modulus E (MPa)	10,000	92.7	-
Poisson ratio ν	0.3	0.22	-
Friction angle φ (°)	-	37.7	-
Interface properties
Friction angle φ (°) [43]	25.0	37.7	37.7
Normal stiffness kn (Pa/m3)	2.8 × 1012	-	2.9 × 1010
Tangential stiffness ks (Pa/m3)	2.8 × 1012	-	2.9 × 1010
Numerical parameters
Timestep Δt (s)		≈ 5 × 10−6
Rayleigh viscous damping ξ		10%

shows all the numerical parameters. Similar to Section 3, the viscous damping ξ is defined by preliminary simulations showing a cyclic strain of 0.002% in the backfill, thus justifying a value of viscous damping equal to 5% [39], again multiplied by 2 since only mass-proportional damping is considered [35].

According to Equation (1), the maximum allowable mesh size is 1.6 m for the backfill and 14.5 m for the units, which are again much larger than the actual numerical mesh size of 17 mm. As shown in Figure 9, the top block’s horizontal displacement d_top is monitored throughout the simulation and serves as the failure indicator (Section 4.2). Figure 9 also presents the numerical boundary conditions (Section 2), while the width W of the backfill is chosen to be equal to 5 m. Though this dimension may be considered small, Figure 10 illustrates that using a twice as large a width gives a similar wall response at the top when the wall is subjected to a harmonic input (see Section 4.2 for a description of the input signal). These simulations also validate the numerical procedure, analysing the deformation of the backfill at wall collapse (Figure 11). Though a slight detachment of the contact between the backfill and wall is noticeable at the very top of the wall, the backfill mostly follows the movement of the wall. For comparison, a similar amount of detachment is noticed for the simulations of the scaled-down experiments (Section 3). This continuous contact is sought to ensure an accurate transfer of forces at the wall–backfill interface, even up to the system’s collapse.

4.2. Dynamic Input Signal

As is common in earthquake engineering, the present study adopts horizontal harmonic signals (with different frequencies) as dynamic inputs to simplify the chaotic nature of actual ground motions: it allows a more straightforward (with fewer parameters) definition of the dynamic inputs [44,45,46,47].

The maximum magnitude of earthquakes expected in France does not exceed M = 7, the most significant event having had an estimated magnitude of M = 6.2 (Lambesc, 1909). According to the literature, the duration of earthquakes may be linked to their magnitude [48,49,50,51,52]; hence, to be as close as possible to the metropolitan French seismicity, the duration T of the input signal aims to reach T = 15s [45]. This appears to be slightly lower than the outcomes of the formulas proposed in [49,50,51] but keeps practical durations of the numerical simulations. In addition, no earthquakes with a longer duration T have been recorded in metropolitan France.

Finally, the intended input signal S₁(t) shape is shown in Figure 12a (analytical equations similar to Equation (4)). It is composed of an increasing and a decreasing part, idealising seismic accelerograms. Based on the analysis of several actual recorded earthquakes [53,54,55], the duration of the decreasing part is defined as twice the increasing one. One must note that the exact time of the signal peak A_peak (hereafter named peak acceleration) will not significantly influence the structural response, as the total energy of the signal is not affected by this parameter. The resistance of a given DSRW is defined as the maximum peak acceleration A_peak that leads to wall stability after applying an input signal S₁(t) of signal peak A_peak (Figure 12a). To obtain a first rough estimate of the resistance of a single DSRW, at least two 15 s simulations with different A_peak values are necessary: one leading to a stable wall and one to collapse.

Envisioned as a less costly parameter, A_col is the maximum acceleration that still guarantees the wall’s stability after application of an input signal with only a single increasing ramp S₂(t), truncated to A_col (Figure 12b). The definition of the instant of collapse (Figure 12) is derived from a preliminary simulation of DSRW1 with f = 5 Hz (Figure 13). In particular, DSRW1 is found to be stable at simulation time t₁ (Figure 13a,b) but unstable a few instants afterwards, at time t₂ (Figure 13c,d). According to Figure 13e, a top wall displacement d_top equal to 20% of the wall height delimits these stable and collapse states well. Therefore, A_col is designated as the acceleration of S₂(t) when d_top exceeds 20% of the wall height.

Finally, based on DSRW1, the ratio between A_peak and A_col is estimated through one simulation with S₂(t) as the input signal and two simulations with S₁(t) (A_peak = 0.6 × A_col and A_peak = 0.7 × A_col, respectively). Note that the three increasing ramps have the same slope, resulting in slightly different durations T for S₁(t), though these are very close to the intended duration of 15 s (Figure 14a,b). As depicted in Figure 14, in the first case (A_peak = 0.6 × A_col), DSRW1 remains stable throughout the motion, while in the second case (A_peak = 0.7 × A_col), it collapses. Therefore, 60% is relevant to characterising the A_peak to A_col ratio. In Section 4.4, the numerical simulations apply only S₂(t) to the walls to identify A_col. Next, if we assume that it does not significantly depend on the system’s parameters (wall, frequency and signal ramp), the ratio 60% gives an estimate of A_peak. In particular, this strategy considerably reduces the computational time (by three to four times) compared with a direct evaluation of A_peak using only S₁(t) signals.

4.3. Different Configurations Tested

First, four geometries of DSRW are investigated numerically (Figure 9 and Table 5). They are subjected to input signals S₂(t) with a frequency of f = 5 Hz (Figure 12b). To keep approximately constant simulation durations (linked to the intended duration of the idealised seismic signal S₁(t), T = 15s), the collapse time t_col is always intended to equal 8.3 s. Therefore, as the walls have different intrinsic resistances, the respective slopes of S₂(t) need to be different. Herein, for each wall, the slope is defined so that its pseudo-static resistance is always reached at 2.5 s, which leads to approximately the same duration for the simulations (Section 4.4).

Next, DSRW1 is subjected to input signals S₂(t) of different natural frequencies (f = 3, 5 and 10 Hz). Frequencies larger than 10 Hz are not critical to DSRWs (Section 4.4), while frequencies smaller than 3 Hz are not representative of the frequency of earthquakes in metropolitan France. Finally, a vertical dynamic signal is added to the horizontal one S₂(t) for two additional simulations involving DSRW1. The vertical input signal has the same shape as the horizontal input S₂(t) and half the amplitude, according to the ratio proposed by Eurocode 8 [27,32]. In the first simulation, the positive peaks of the horizontal and vertical accelerations are in-phase, while in the second one, they are in opposite phases. In these cases, the slope of S₂(t) is recomputed, accounting for the vertical acceleration that updates DSRW1’s pseudo-static resistance. The outcomes are presented in Figure 15 and compared with the original case. The most critical case appears when the signals (horizontal and vertical) are in opposite phases; in other words, when the maximum (positive) horizontal load coincides with the minimum (negative) vertical load, i.e., with the minimum equivalent natural vertical gravity.

Table 7. Collapse acceleration A_col and time t_col for all simulations performed. Here, t_col is almost constant for all the simulations (with a range between 7.6 s and 9.5 s, and a coefficient of variation of 7%).

DSRW Name	DSRW1					DSRW2	DSRW3	DSRW4
Frequency (Hz)	3Hz	5Hz	10Hz	5Hz	5Hz	5Hz	5Hz	5Hz
Vertical acceleration	-	-	-	in-phase	opposition	-	-	-
tcol (s)	7.6	8.5	9.5	8.9	8.6	9.4	8.3	8.5
Acol (m/s2)	4.83	5.38	6.03	5.34	5.19	1.83	5.06	4.63

summarises the collapse acceleration A_col and time t_col for each simulation. As expected and intended, t_col is only slightly affected by variation in the parameters. On the contrary, the acceleration at collapse A_col is more affected: in particular, DSRW2, the only wall with a sloped backfill, has a much lower capacity, illustrating the negative effect of backfill slopes on the stability of DSRWs, similar to retaining walls in general.

4.4. Comparison of Pseudo-Static Predictions and the Estimated Behaviour Factor r

Table 8. Comparison of the collapse acceleration A_col, the estimated peak acceleration A_peak and the pseudo-static resistance A_pseudo, along with the estimated behaviour factor r for every single numerical case in this study. When the vertical acceleration is accounted for, the displayed accelerations only show the horizontal components, the vertical ones being always half of the latter.

DSRW Geometry	DSRW1					DSRW2	DSRW3	DSRW4
Frequency (Hz)	3Hz	5Hz	10Hz	5Hz	5Hz	5Hz	5Hz	5Hz
Vertical acceleration	-	-	-	in-phase	opposition	-	-	-
Acol (m/s2)	4.83	5.38	6.03	5.34	5.19	1.83	5.06	4.63
Apeak = 0.6 × Acol (m/s2)	2.90	3.23	3.62	3.21	3.11	1.10	3.04	2.78
Apseudo (m/s2)	1.58	1.58	1.58	1.50	1.50	0.49	1.36	1.52
r = Apeak/Apseudo	1.83	2.04	2.29	2.14	2.07	2.26	2.00	2.04

estimates the peak acceleration A_peak for each DSRW, considering the ratio of 60% found in Section 4.2. In addition, Table 8 also shows the pseudo-static resistance A_pseudo of each wall computed using an analytical approach based on the limit-equilibrium of a Coulomb’s wedge of soil [18,27]. The properties used in the analytical model are the same as for the numerical simulations (Table 5 and Table 6). Furthermore, the stones are assumed to be perfect, which leads to a high interface stiffness in the numerical simulations (Table 6) and no internal rotation of the stones in the analytical approach (see the details in [56]).

Finally, Table 8 estimates the behaviour factor r for the investigated DRSWs as the ratio between the peak acceleration A_peak, and the pseudo-static resistance A_pseudo. On average, the behaviour factor r equals 2.08, with a minimum value of 1.83. Therefore, to fit the thresholds (1, 1.5 or 2) recommended by Eurocode 8 for retaining walls [32], the authors propose r = 1.5, which is conservative, as it is intended to be used for practical design/assessment. Furthermore, the latter can readily be used in conjunction with a pseudo-static approach to give fast practical recommendations on-site. This is of particular interest when many DSRWs need to be assessed/designed, requiring fast methods where performance-based approaches are less competitive, as they require more computational time and often need several expensive time-history analysis to give reliable results.

5. Conclusions

This work presents dynamic numerical simulations involving dry-stone retaining walls (DSRWs). Firstly, a discrete element modelling (DEM) tool is validated against a scaled-down laboratory mock-up subjected to horizontal harmonic shaking. The outcomes show good agreement regarding the collapse mechanism and the shape of the response while predicting a slightly lower failure capacity, thus being conservative. It is also noted that the agreement holds for various slenderness and frequency values.

Next, time-history dynamic simulations involving full-scale DSRWs are carried out. The associated data are taken from the literature. The dynamic input signals are harmonic, aiming to model seismic motions in a simplified and fully controlled manner. The signals’ properties are based on the most critical seismicity of French metropolitan areas. However, the outcomes also hold for any country with a similar (or lower) seismicity. The effects of the frequency and the vertical acceleration are investigated. It is found that the frequency of the input signals is the most critical characteristic of the motion, and it noticeably modifies the maximum acceleration supported by the DSRW.

Finally, the behaviour factor r is deduced for each simulation, comparing the pseudo-static resistance with that obtained from the time-history dynamic simulations. An average factor r = 2 is found, with a minimum value (for low frequencies) of 1.8. Therefore—and this is the novelty of the present study—the authors advise a conservative value of r = 1.5 for DSRWs, chosen among the recommended values by the Eurocode 8 standard. This value can readily be used in a practical framework to assist engineers and practitioners in designing new (or assessing existing) DSRWs following a classical pseudo-static approach for fast analyses.

Future work will involve sets of artificial and recorded seismic motions to confirm the recommended value for the seismic behaviour factor r.