Seismic rocking simulation of unreinforced masonry parapets and façades using the discrete element method

When subjected to high‐intensity out‐of‐plane seismic loading that results in the threshold value of static stability being exceeded, unreinforced masonry parapets and façades form mechanisms that exhibit rigid body rocking motion before collapse. The discrete element method was used to perform rocking simulations which were compared to experimental observations and analytical models to obtain relationships that enable the stiffness proportional Rayleigh damping parameter (β) to be estimated from the coefficient of restitution. The seismic performance of parapets and façades was compared using incremental dynamic analysis, and similar performance was observed between the compared models, confirming the validity of the developed expressions to obtain β. The gap size generated after failure between the façade and the return wall was parametrized and the outcome revealed the hazardous situation that arises when the connectivity between the uncollapsed façade and the adjacent return walls is damaged during an earthquake. Pulse‐like accelerations, generally associated with near‐fault earthquakes, caused collapse of parapets and façades at substantially lower accelerations when compared to earthquakes with oscillations that were evenly distributed in time. It was concluded that when assessing the seismic vulnerability of a rocking façade located near a fault, a capacity equal to that of a free‐standing block should be used, validating the assumption made in some seismic assessment guidelines of not considering the effect of return walls when assessing the façade capacity. For façades located far from a fault a range of earthquakes that represent the site conditions needs to be included in any dynamic capacity assessment.


INTRODUCTION
The rocking motion of a rigid body was first studied in the 1960s by Housner, 1 who established the equation of motion for a rocking block as being similar to that of an inverted pendulum. Following the findings of Housner, 1 authors such as Aslam, et al. 2 and Yim, et al. 3 observed that, for rather large displacements, minimal uncontrollable changes in the ground motion resulted in large changes in the response of the block. One of the most interesting features of rocking behavior when close to overturning was the inability to find clear trends via changing either the parameters of the strong motion input or the block parameters, as shown by Yim et al. 3 For example, scaling the intensity of a given earthquake record could change the rocking behavior of a block, without increasing the rocking response of the block. Yim et al. 3 also found that artificially created earthquakes with the same spectral density and time function imposed different behavior on a rocking block. Similarly, the addition of different vertical earthquake motions applied to rocking numerical simulations sometimes helped to stabilize the block and sometimes contributed to greater instability. As demonstrated by Maiorano, et al., 4 the size and slenderness ratio are good indicative parameters to assess the overturning potential of blocks when subjected to earthquakes, but exceptions were often found where stability did not systematically increase with increasing size and decreasing slenderness ratio. Also, reduction of the coefficient of restitution did not necessarily result in reduced response due to a consequence of greater energy dissipation, as would be intuitive. In order to understand at a fundamental level the different phenomena that could arise when blocks rock, Spanos and Koh, 5 Tso and Wong 6 and Wong and Tso 7 studied the effect of harmonic excitation on blocks and found steady states in their response. Such response phenomenon helped to explain the rocking behavior of blocks when subjected to more complex inputs (e.g., earthquake excitation) and reflected the complexity of the rocking phenomena, leading to further development of several models that aided in understanding of the problem. [8][9][10][11] Later research using incremental dynamic analyses (IDA) of 2D rocking blocks 12,13 further helped understanding of the rocking phenomena, although exceptions to recognizable trends were identified. It was established that blocks that rotated close to the point of overturning required very little additional energy to overturn, making the solution conditionally unstable, as illustrated by Yim et al. 3 via the use of potential energy versus rotation plots for two-sided rocking blocks and by Sorrentino et al. 14 for one-way bending walls. This conditional instability indicates that during an earthquake, once significant rocking of the block has developed, very small differences in the ground motion can greatly affect the response of the block. Due to the highly nonlinear nature of this problem, rocking blocks subjected to higher demand than the lowest theoretical level of excitation necessary to generate overturning resulted in so-called structural resurrection. Groups of safe and unsafe cases based on signal input parameters could always be identified within the IDA predictions, and sinusoidal signal pulse motion was shown to have a similar influence on rocking blocks as did near fault earthquakes. Therefore, depending on the combination of the frequency and amplitude of the pulse, the block might overturn or might not. A large number of simulations of ranging pulses have been plotted to illustrate the vulnerability of rocking blocks, with the plot described as an overturning envelop. 15 Rocking blocks subjected to ground motions receive a sequence of pulses having a combination of frequencies and accelerations, meaning that due to previous pulses within the earthquake record, the block pre-possesses a non-zero displacement and velocity at the arrival of a pulse, such that the block may overturn at a smaller intensity measure than that predicted for the single pulse case. 16 Ishiyama 17 developed approximate geometrical criteria to determine the minimum acceleration needed to start rocking and the minimum velocity that will produce overturning of slender blocks. The overturning results of rocking simulations for slender blocks that were produced by Ishiyama 17 were related to the peak ground velocity (PGV) ground motion, finding fair agreement. More recently, Sorrentino et al. 18 and Maiorano et al. 4 also demonstrated that the strong motion parameters that best agree with overturning of slender blocks were PGV and seismic wave energy flux density.
When subjected to an earthquake, unreinforced masonry (URM) façades fail by forming mechanisms that behave according to rigid body rocking mechanics before they collapse. Two-sided free-standing rocking block (hereafter referred to as parapet) mechanisms have been widely studied and their behavior has been forecasted mostly using linear, viscous single degree of freedom 19,20 oscillator (SDOF) models and through derivation of the formulation for an inverted pendulum as developed by Housner, 1 hereafter referred to as the Housner model. According to Makris and Konstantinidis, 21 SDOF models do not accurately represent the dynamic behavior of a rocking rock, and to assess URM rocking structures Lagomarsino 22 showed the differences between the Housner model and a nonlinear SDOF model where the elastic behavior before the rocking onset was included.
Rocking blocks as a representation of free-standing URM parapets and one-way bending walls have been widely studied by experimental testing and by using numerous models, 14,[22][23][24] but different configurations other than simple blocks have been observed after earthquakes and in most cases the boundary conditions do not relate to free-standing and one-way bending walls. Generally, a complete oscillation during rocking is not possible and instead one-sided rocking (1SR) occurs. 25 To the knowledge of the authors, no studies have yet been made that compare the seismic behavior of 1SR URM parts and how this behavior differs from free-standing rocking. Hence this lack of research is addressed here, placing special attention on damping factors in order to proceed with more complex configurations and boundary conditions in future studies. A 1SR configuration is used because it is the most common and most simple mechanism model of an overturning façade, which from hereon is simply referred to as a façade. The earthquake response of elevated façades and parapets is significantly influenced by the principal dynamic response of the building that they are connected to. For example, for the same earthquake and at the same location, the seismic input to a parapet located at the top of a multistorey building will be different to the seismic input for the same free-standing wall supported directly on the ground. [26][27][28] Recognizing the role of building amplification, in the reported study the difference between parapets and façades was studied subjected to application of the same (non-amplified) ground motion without the building filtering the signal. However, the influence of the dynamic behavior of the building in the rocking capacity of the façade was taken into account. Additionally, in some seismic assessment guidelines (e.g., NZSEE 29 ) no clear method to assess overturning façades is specified, leading to practitioners instead using the free-standing parapet methodology as a valid tool to assess façade vulnerability, but with insufficient supporting evidence. The study reported here was performed using a discrete element method (DEM) nonlinear model implemented in the software 3DEC 30 that allows for 3D simulation of nonuniform rocking blocks with multiple boundary condition options. DEM has been widely used for URM simulations and has been proven to accurately replicate collapse mechanisms and capture URM dynamic behavior. [31][32][33][34][35] For completeness, the final comparison between façade and parapet rocking behavior includes the Housner model properly calibrated as in previous studies, 18,36 offering a multi-model validating aspect to the study.

RELATIONSHIP BETWEEN SP RAYLEIGH DAMPING AND COEFFICIENT OF RESTITUTION
DEM simulations performed with 3DEC follow a Rayleigh damping configuration either using mass proportional (MP) or stiffness proportional (SP) damping for dynamic analysis. Rayleigh damping makes use of the constants α (MP) and β (SP) to construct the damping matrix. The damping ratio ζ n for any circular frequency that results from the application of α and β can be found in Bathe and Wilson 37 as: Galvez et al. 38 studied the possibilities and challenges of an MP damping approach to simulate 3D DEM rocking blocks versus an SP damping approach, reaching the conclusion that SP damping provides greater benefits in terms of accuracy and parametric stability at the cost of higher computation time. For this reason, SP damping was further implemented and studied in the study reported here.
One of the most thoroughly studied models used to simulate the rocking behavior of free-standing blocks is based on the derivation of the Housner 1 formulation, which relies on the coefficient of restitution to account for the energy lost in every impact. The coefficient of restitution e is the ratio between angular velocities after and before each impact. Conservation of angular momentum before and after each impact gives the theoretical value, derived by Housner, of e H (see Equation (2), where = tan −1 ( ∕ℎ), b is the thickness and h are the height) valid for a homogeneous parallelepiped block.
Assigning of damping parameters in DEM can be a challenging process, especially when rocking blocks are to be simulated. DeJong 39 provided a procedure to assign values of β to free-standing rocking block DEM simulations, and Sorrentino et al. 36 proposed a reduction, based on experimental observations, of the value of e calculated from Equation (2). A mathematical correlation between e and β, delivered after a parametric analysis, would represent a powerful tool for DEM users that would allow them to avoid a lengthy and sometimes inaccurate process for the selection of damping parameters. Relating e and β has only been attempted by Peña et al., 40 but major drawbacks were found in the study, where the stiffness required for implementation of the relationship was extremely low, creating numerical instability in the model.
The damped motion of rocking bodies is heavily influenced by contact interface discretization, with relevant literature suggesting that between 3 and 16 sub-contact points are required to simulate the transition between the elastic and rocking phases, 12,38,[41][42][43][44][45] with agreement that a larger number of contact points delivers a smoother transition but substantially increases the computation time. The difference in the dynamic performance between blocks with different numbers of integration points was investigated by Galvez et al. 38 using SP damping, with more damped oscillations being observed when fewer contact points were implemented. This difference in damped behavior is a consequence of the decrease in the spring stiffness with an increasing number of contact points and therefore the rocking amplitude of the block increased with more contact points. Convergence in performance of the block was obtained when between 10 and 15 sub-contact points were implemented, but a significant change in computational time was observed when adding more contact points. A single time-history computation using 2 sub-contact points took 6 h, while the same computation with 3 contact points took 13 h. In order to study a model similar to that originally proposed by DeJong 39 and to perform a large number of simulations to obtain the correlations between e and β and develop the IDA presented in the reported study, a model with 2 sub-contact points was adopted. It is noted that researchers that perform simulations with large time demanding models usually use a 2 sub-contact point configuration. 33

Two-sided free rocking blocks (Parapets)
Free rocking of blocks with different geometries was simulated using a single discrete rigid element released from a displaced position, applying the value of β obtained by following the method developed by DeJong 39 for each simulation.
, where E is the Young's modulus and m is the mass of the block) for impact of the blocks against the ground, which should be critically damped as happens in nature. Once that ω crit is calculated, Equation (1) is used to obtain the corresponding β, keeping α as zero for the SP damping approach (further information can be found in Galvez, et al. 38 ). This value of β was parametrized keeping the rest of the input parameters unchanged. Stress along the nonlinear joint interface was calculated using the Mohr-Coulomb failure criterion with tension cut-off and a shear stress limit. The density of the blocks was ρ = 1800 kg/m 3 and Young's modulus E = 1800 MPa, which was used to calculate the interface joint stiffness j kn = E/h, and using elastic relations j ks can be calculated from j kn as ks = kn ∕2 ⋅ (1 + ), where v is the Poisson ratio which was taken in this study as 0.25. To allow rocking the base joint must open completely without cohesive resistance. Therefore, in all the rocking simulations reported herein the tensile strength and cohesion were set to 0, while to avoid slippage the friction angle was equal to 50 • . At the same time, the rocking equations of motion developed by Housner 1 were implemented in the software Matlab, including a three-branch moment-rotation law, and the DEM simulations were matched by varying the value of e. The most appropriate match was selected by obtaining the minimum weighted mean error of the time-history results of both models (the same procedure as in Shawa et al. 42 ). The values of β and e for the matched simulations can be seen in Figure 1, together with the trends for an increase of the applied damping. No changes in rocking behavior was expected when the length of the block (perpendicular to the rocking motion) was modified while keeping the height and thickness unchanged, because the polar moment of inertia does not vary. As part of the geometry parametrization, the length was changed and the results were included as smaller markers only in the chart for 0.1 m thickness ( Figure 1A), with the results remaining within the previously obtained trends. The trends obtained for the geometry variation followed a linear distribution that intercepted the horizontal axis at the theoretical value of e H (Equation (2)) on the horizontal axis as: where x is a fitting factor. The trend intercepting the horizontal axis suggests that when all external damping is turned to zero, the model is capable of simulating rocking behavior with damped motion equal to the performance of the Housner model with e H applied. However, this fortunate constraint for the regression is mainly theoretical because when zero viscous damping is selected, high-frequency undamped numerical noise and impact bouncing trigger unrealistic behavior. 38,[52][53][54][55] If x is plotted against the h/b ratio of the blocks, the power trends observed in Figure 2A can be expressed as: where y is, as x, a fitting factor. For any block dimensions, and requiring only that the value of b is known, y can be obtained by the relationship shown in Figure 2B, which will permit x to be obtained using Equation (4), and finally the relationship for β and e is obtained according to Equation (3). In the formulation developed by DeJong, 39 β is dependent on Young's modulus (E) and density (ρ). For this reason, the block stiffness and weight were also parametrized and the response was matched with the corresponding e as seen in Figure 3 with the trends for the relationship between β and e found to not change with varying E and ρ. The values obtained using the formulation given by Peña et al. 40   Obtained correlation comparison with existing literature Tomassetti et al. 56 previously correlated e with damping ratio in an analytical viscous damping model. Even though both Tomassetti et al. 56 and the model studied herein use viscous damping, and one of the models studied in Tomassetti et al. 56 also had the SP damping ratio implemented, the trends in parameter variation were different. While Tomassetti et al. 56 found that the main factor that resulted in varying of the distribution of e was the initial stiffness of the rocking block and that varying the block shape (h/b) had no influence on the e relationships, the results plotted in Figures 1 and 3 showed that blocks dimensions had a major impact on the trends, and that varying E for a fixed geometry had no influence on the relationship trend. There is a series of differences between both models that justify the difference in observed trend. First, the Tomassetti et al. 56 model was based on the solution of the piece-wise formulation of the equation of motion of a SDOF spring-mass model to compute the rocking behavior of a parapet. Tomassetti et al. 56 obtained the damping ratio distribution by calculating the system circular frequency at initial stiffness ( 1 ) and assigned a damping ratio at that frequency (ζ 1 ). When 1 and ζ 1 were calculated, the damping distribution was obtained as ζ n = ζ 1 · ∕ 1 and the relationship with e was made through the expression ζ 1 = −0.218 ⋅ 1 −0.195 ⋅ ln( ), where a 1 is the instability displacement ( ) multiplier to obtain the elastic displacement at peak force ( 1 ). However, the numerical model implemented in the reported study is a DEM, and the damping distribution with frequency was obtained by calculating and using Equation (1) to obtain β and the relationship between and ζ n . Secondly, the initial stiffness parameter that Tomassetti et al. 56 implemented for the spring of the model was 1 = 1 ∕ 1 ⋅ , whereas the stiffness parameters in the reported study correspond to the interface joint stiffness (j kn ) between the block and the floor and is based on j kn = E/h. For varying geometries, k 1 is affected by change in thickness (u ins and F 1 ) and height (F 1 ), while the model proposed here only varies with height. Therefore, for a SP damping configuration, the impact of geometry variation on damping is different. And finally, when Tomassetti et al. 56 relate ζ 1 with e, the relationship relies on a 1 , which is a known fixed number for parapets, while Equations (4) is affected by aspect ratio (h/b) and wall thickness b (see Figure 2B). In summary, the two models and the procedures used to assign damping are inherently different, and hence any attempt to compare the damping distribution with frequency or with e between the two models is deemed to give different results for each model.
Recently Vlachakis et al. 57 published a study expanding the research performed by Tomassetti et al. 56 by obtaining a correlation between e and ζ n to be implemented in numerical models. Similarly, to the model reported herein, Vlachakis et al. 57 included correlations to use viscous damping strategies to be applied on both two-sided and one-sided rocking problems. Damping was treated as a dashpot applied at the nodes on the interfaces, reducing the velocity of the block, and a procedure similar to that proposed herein was adopted to correlate e and ζ n . Vlachakis et al. 57 found that the geometrical relationship h/b greatly influenced the relationship between e and ζ n , similar to the findings of the study reported here. In the study performed by Vlachakis et al. 57 the relationship developed to correlate e and ζ n was applied to two numerical models: a rocking block discretised into an FEM deformable mesh for inclusion in the software Abaqus CAE, and to a rocking block discretized into an assembly of elements simulating the individual bricks when using 3DEC. However, the model presented herein used a single block with elastic deformations lumped at the base interface, which constitutes a considerable difference with respect to the study undertaken by Vlachakis et al. 57 Additionally, for the two-sided rocking problem Vlachakis et al. 57 used a 2D FEM model with a nonlinear interface between the rocking elements and the ground F I G U R E 4 Theoretical and experimental coefficient of restitution from various research and used two separate interfaces for ground and return walls for 1SR, which allowed for specification of different ζ n for the return wall and ground interfaces, not supported by the reported model. The bottom interface stiffness (k nV ) in the Vlachakis et al. 57 model was reported to have a large influence in the relationship between e and ζ n , while E was observed to not be influential in the reported relationship between e and β as seen in Figure 3. The comparison of parameters k nV and E is not straightforward due to the conceptually different understanding of the boundary contact interface. In addition, β is obtained from ζ n using Equation (1) where is needed, and is obtained from k nV according to Vlachakis et al. 57 and Multiple authors have studied the limitations of the value of e obtained with Equation (2) when compared to experimental tests but, with relation to URM blocks, the most relevant studies were performed by Sorrentino et al., 36 Costa et al., 58 Giaretton et al., 24 Degli Abbati and Lagomarsino, 59 Giresini et al. 60 and Giresini et al. 61 From these studies it was found that the experimental value of e always oscillated between 85% and 95% of the theoretical estimation of e (see Equation (2) and Figure 4). Such observations could appear contradictory to the improved formulation for the estimation of e as developed by Kalliontzis et al., 62 who increased the value of e from Equation (2) based on experiments reported by Muto et al., 63 Ogawa, 64 Priestley et al., 20 Aslam et al., 16 Lipscombe and Pellegrino, 10 Peña et al. 40 and Fielder et al. 65 (see Figure 4). However, none of the studies considered by Kalliontzis et al. 62 specifically addressed URM blocks. One of the possible differences between the URM tests included in Figure 4 and the tests described by Kalliontzis et al. 62 could be the reported damage at the rocking corners of URM blocks that could lead to a displaced rocking pivoting point O, but, as evident from Equation (2), reducing b at the interface would only make e larger. The difference of e between solid blocks and URM blocks could be linked to the rocking interface which was found to have a major influence on the rocking behavior. 66 The values of e obtained using Equation (3), as equivalents to the values of β obtained using the DeJong 39 formulation, were plotted in Figure 4 with three different Young's modulii and two different densities. For the first time to the knowledge of the authors, simulations using β can be displayed with confidence against other rocking research using different damping methods. The plots showed that the application of Equation (3) was reasonably within the experimental limits, but that for some h/b relationships the range of results was large.

One-sided free rocking blocks (Façades)
Parapets and façades can be approximated by a rocking block having an equivalent overturning pushover capacity. A façade is restricted to rock only around one side (1SR), having different dynamic rocking behavior and boundary F I G U R E 5 Parametric study of damping for URM façades conditions than a parapet. The rocking behavior of façades using a single discrete rigid element as the rocking block was studied in order to develop a relationship between e and β for façades. The 1SR configuration was first considered for a generic block by Hogan, 9 followed by Sorrentino et al., 25 who studied the earthquake behavior of 1SR URM blocks using a model based on the Housner 1 rocking equations of motion. Lacking enough experimental information to calibrate energy dissipation, an extensive number of parametric simulations led to a greater understanding of the nature of the problem. Later, an experimental campaign performed on 1SR free rocking led to improved comprehension regarding the fundamental relevance of energy damping in such problems. 36 Classical approaches 1,9,67 were considered to no longer be valid, and an adapted formulation to account for the distinct damped motion caused by multiple impacts of the block against the return walls was developed as: Sorrentino et al. 36 observed during the experimental campaign that the experimental coefficient of restitution, shown in Equation (6), 1 of the third bounce after impacting with the foundation and the return walls was generally the closest to e 1SR in Equation (5). In Equation (6), |θ n | is the maximum absolute rotation at the n th bounce.
Further study of the experimental data from Sorrentino et al. 36 showed that, in addition to the third bounce being similar to e 1SR , the first impact e exp was between 64% and 72% of e 1SR and the second impact was between 74% and 86% of e 1SR . Such findings allow the experimentally correct solution to be identified and subsequently compared against numerical simulations with varying parameters by analyzing the e exp of the three first bounces. Therefore, a parametric analysis was performed to explore a possible estimation of β for 1SR walls by finding the results that best matched the aforementioned experimental percentages of e 1SR . A range of β was applied to four batches of free rocking façade simulations with varying block geometry [height (h) and thickness (b)], Young's modulus (E = 1800 and 4050 MPa, see Figure 5), friction angle equal to 30 • (value obtained from Sorrentino et al., 36 where no sliding was observed and based on Navier, 68 Rankine, 69 Colombo, 70 and Meyer et al. 71 ) and different gap sizes between the return walls and the façade. Simulations with high values of β close to the value calculated using ω crit resulted in overdamped behavior that wrongly restricted the bouncing of the block after the first impact ( Figures 6A, B), not allowing the façade to impact the return walls for a second time. A lower value of β can be obtained by adapting the formulation given by DeJong 39 for vertically stacked blocks (ω Vcrit = ω crit · √ 2) by considering the 1SR problem as a series of horizontally stacked blocks. When substituting in the ω crit formulation, height would be thickness and vice versa, resulting in = √ 2⋅ ⋅ℎ 2 ⋅ (further information can be found in Galvez et al. 38 As seen in Figure 6C, some overlap between the façade and the return walls was only observed for β < 0.006. The values of β used as part of the parametric study can be seen in Figure 5A for  Figure 5A represent the values where overlap between the façade and the return walls (as per Figure 6C) was observed. As seen in Figure 6, these simulations showed a more realistic rebound behavior. The DEM simulation results using β that best matched analytical simulations with experimental percentages of e 1SR as previously explained are circled in Figure 5A, including the trendline of the best matching results for each h/b block combination. As expected, the best fitting β was found to always be closer to the β values obtained using ω Hcrit . The best fitting trends for different h/b of façades with b = 0.11, 0.30 m, and E = 1800, 4050 MPa, as observed in Figure 5B, were studied to obtain a relationship for estimating β 1SR as follows: wherein b has units of m, and E of MPa. As E reduces and b increases the energy damping needed to perform a realistic 1SR simulation increases. All the trends studied herein were experimentally validated in Galvez et al. 38 and their consistency was tested in Galvez et al. 72 The gap between return walls and façades that influences damping of the rocking behavior was also studied by performing DEM simulations with varying gap sizes (0, 0.3, 1, 3, 4, and 6 mm) as seen in Figure 7. Equation (6) was used to identify the level of energy dissipation in each simulation by obtaining an e gap for each simulation, and afterwards was compared to the e 0 value corresponding to the coefficient of restitution with zero gap between the return walls and the façade. The

COMPARISON OF PARAPETS AND FAÇADES
After understanding the role of damping in the DEM simulation of rocking behavior of façades and parapets the seismic performance of four blocks with different aspect ratios was studied. Each block was modelled as an URM parapet and façade in order to compare the dynamic capacity between each collapse mechanism that, while different, most times are assessed as equal (i.e., no differentiation is made in the NZSEE 29 seismic assessment guidelines). Ideally, a large range of earthquakes would be included to study the seismic performance of the blocks. Because of the highly demanding computation time required in DEM, only four earthquakes with the widest range possible of seismic features (see Table 1) were selected to perform IDA on parapets and façades. The records BUCAR0, 1ST280, and RRS228 were obtained from Sorrentino et al., 18 and the HVSC1 earthquake record was selected from the Canterbury earthquakes in Christchurch 2011. 73 This selection was made by considering a range of earthquake moment magnitude (M W 6.3-7.5) and distance from the surface projection to the source (0.1-150 km), as well as a range of station soil type as defined by Eurocode 8 74 Figure 8A), and the sine cycle (SC, Figures 8A-D). 77 Such pulse shapes are known to induce different responses on rocking blocks 77 and are well represented in the selected earthquake records. BUCAR0 and RRS228 are pulse-like earthquakes that have a large pulse near the beginning of the time-history. These earthquakes are generally typical of near-fault regions within 20 km from the source, 78,79 but it has been proven that several exceptions can be found in historical records. 80,81 On the contrary, 1ST280 and HVSC1 have evenly distributed motion along time.
The four blocks were modelled with E = 1800 MPa, ρ = 1800 kg/m 3 , friction angle equal to 30 • and conhesionless contacts. The dimensions of the blocks chosen to have a distributed range of aspect ratios (h/b) were h = 1.50, 3.00, 6.44, and 8.00 m and b = 0.170, 0.250, 0.375, and 0.375 m respectively. The one and two-sided rocking models used in Shawa et al. 42 and Sorrentino et al. 36 were also included in the reported study for completeness and to increase confidence in the results. Other than dimensions and density of the block, the one-sided rocking model required the nondimensional displacement parameters Δ 1 , Δ 2 and Δ u . Such parameters control the three-branch moment-rotation law 42 and were calibrated to match the pushover analysis performed with the DEM model of each block geometry (see Figure 9). A three-branch law (either moment-rotation or acceleration-displacement) was proposed in several analytical works (e.g. Doherty et al., 82 Derakhshan et al., 83 Sorrentino et al., 84 Godio and Beyer, 85 Tomassetti et al., 56 Padalu et al., 86 Baraldi et al. 87 ), with such a law not found in standards, although in the Commentary to the Italian Building code 88 the simplified two-branch law F I G U R E 8 Pulse shape sequences (blue: original acceleration; black: low-pass filtered acceleration; red: seismic pulse shape).

F I G U R E 9
Pushover analyses of each wall geometry and the three-branch linearisation to calibrate Δ 1 , Δ 2 and Δ u . discussed by Lagomarsino 22 is included. A more detailed four-branch analytical law can be found in Ferreira et al. 89 Each return wall of the façade modelled using DEM included four cohesionless blocks positioned at the back of the façade, with the ground motion applied to the base and with the amplified motions input to each block positioned above the base, as described in Galvez et al. 72 IDA results are shown in Figure 10, for the different geometric configurations and for the four earthquakes studied, in terms of acceleration and normalized maximum displacement of each simulation (d/b). It can be seen that regardless of the parapet and façade dimensions, higher acceleration was needed to induce large rocking displacements when subjected to 1ST280 and HVSC1 than when subjected to RRS228 and BUCAR0, due to the façades being more vulnerable to pulse-like earthquakes. For all dimensions, the DEM and the analytical approaches performed similarly for the parapet simulations, whereas for all façade simulations except for the 8 m high parapet, greater resistance was attained for the analytical model than for the DEM model, probably because inclusion of the building dynamic behavior in the return walls of the DEM simulations made the façade more vulnerable. 72 When the parapets and façades reported in Figure 10 were subjected to the BUCAR0 earthquake the presence of the sharp seismic pulse shape in the earthquake record dictated that the peak capacity of the rocking element was governed by F I G U R E 1 0 IDA of façade and parapet using DEM and analytical models. the pulse alone, resulting in a similar capacity for the DEM and the analytical models, except for the 3.00 m high façade that was simulated with the analytical model, which sustained approximately twice the acceleration for displacements close to instability. When the simulated blocks were 1.50 and 3.00 m tall, façades sustained a higher intensity of earthquake scaling than did parapets, probably due to the low building amplification of the earthquakes. Simulated blocks with 6.44 and 8.00 m height showed similar capacity between façades and parapets, coinciding with strong building amplification that worsened the façade performance and resonance of rocking periods for the parapets, 90 except for HVSC1, where parapets were clearly stronger than façades. Structural resurrections 91 were observed in both DEM and analytical models but the results after the first observed collapsed simulation were not included in Figure 10 for clarity.
Based on the observed IDA analyses (Figure 10), a comparison study was performed between the DEM façade simulations (to account for the inclusion of building behaviour from motion of the return walls) and an average of the analytical model and DEM parapet results. An average of the results was used for parapets as the most sensible choice to compare against façades because it was shown that both models were capable of simulating parapet dynamic behaviour and the results were found to be similar ( Figure 10).
All results from the study are shown in Figure 11 by means of comparing the normalised displacements of parapets and façades for pairs of simulations subjected to the same acceleration record. A diagonal line was included in the plots to F I G U R E 1 1 Comparison between facade and parapet seismic performance using DEM and analytical model (▼= negative direction, • = positive direction) represent equal performance of parapet and façades. Any pair of simulations located above the diagonal line indicate a weaker performance of the façade with respect to the parapet, and the opposite applies for pairs of simulations located below the diagonal line. The simulations of BUCAR0 and 1ST280 showed weaker performance of parapets than façades, while the results from simulations using RRS228 and HVSC1 were more scattered, with weaker façade performance being slightly predominant for RRS228 and strongly predominant for HVSC1. Both RRS228 and HVSC1 had an MH pulse shape (see Figure 8) that, regardless of the earthquake being pulse-like or not, caused the motion combination of the façade rebounding back after impact with the return wall, together with a reversal in ground motion, resulting in the one-sided façade being more vulnerable. This phenomena is probably one of the factors why façade overturning was amongst the most common failures observed after the Canterbury earthquake (HVSC1). 92,93

CONCLUDING REMARKS
Discrete element method (DEM) models that incorporated the Rayleigh stiffness proportional (SP) damping parameter (β) were studied to analyse the seismic performance of parapets and façades. As widely recommended in the literature, any study on rocking blocks using Rayleigh damping should use the SP damping approach rather than mass proportional damping that has been shown to lead to non-accurate results. The performance of DEM free rocking simulations was compared to experimental observations and well-known analytical models to obtain relationships allowing β to be estimated from the coefficient of restitution (e). When studying parapet rocking motion, the analytical simulations with an estimated e were taken as a reference to match the DEM simulations. The parapet dimensions were found to have a major impact on the trends obtained when relating β and e, whereas varying Young's modulus and density for a fixed geometry had no influence on the relationship trend. For the first time, simulations using β for a DEM model were displayed with confidence against experimental observations and other free standing rocking numerical research. A relationship was also developed for façades based on experimental observations and parametric analysis for different wall geometries and Young's moduli. The values of β obtained for façades were considerably lower than those for parapets, in the order of 10. 1 The gap size between the façade and return walls was studied with relation to energy dissipation of the one-sided rocking motion, and it was found that in the majority of cases for different gap sizes the rebound motion of the façade was less damped when there was a gap than when not (registering coefficients of restitution up to 1.8 times larger than when there was no gap), illustrating the hazardous situations that arise when the façade bond with the adjacent return walls is damaged during an earthquake, and the urgent need to repair the connection before the arrival of aftershocks. However, the level of damped motion did not follow a trend or an ordered correlation with gap size.
The values of β obtained through the relationships developed for façades and parapets were implemented into the DEM model to perform incremental dynamic analysis using four earthquakes and to compare the predictions with analytical models. Similar performance was observed between the compared parapet models, confirming the validity of the developed expressions to obtain β. The DEM model for the façades included the building dynamic behavior in the return walls that, as expected, lowered the performance of the façades in most cases when compared to the analytical model. Two out of the four earthquake records displayed pulse-like accelerations, which lowered the performance of parapets and façades to between 1/3 and 1/12 when compared to the performance derived for earthquake records having oscillations that are evenly distributed in time. The dynamic capacity obtained with pulse-like earthquakes was close to the static capacity. It was observed that for earthquake records exhibiting the "Mexican hat" pulse shape the capacity of façades reduced due to such pulses generating the motion combination of the façade rebounding from the return wall together with a reversal in ground motion. For other studied pulse shapes, the capacity of façades was higher than that of parapets, as expected. The abovementioned observations led to the conclusion that when assessing the vulnerability of a façade located in a near-field area, at least two types of pulse-like records need to be included in the analysis, with one record having a "Mexican Hat" shape, and a capacity equal to that of a two-sided free standing rocking block should be used. This conclusion confirms the criteria appearing in some seismic assessment guidelines, where it is recommended to analyse façades as parapets without considering the transverse walls impact effect. Elevated seismic capacity might be obtained when assessing unreinforced masonry buildings located outside the near-field region, by performing dynamic rocking analyses including records that represents the hazard of the area. The simulations presented herein were performed using the software 3DEC and an analytical model where the rocking equations of motion developed by Housner 1 were implemented. Therefore, the trends obtained were validated and can be expected to be observed when using other models. The reported study opens multiple possibilities for simulating the collapse mechanisms of URM building element having more complex geometries and boundary conditions, whilst having confidence in the stiffness proportional damping parameters applied.
Future research could focus on increasing the number of records that include pulse-like earthquakes to explore more trends that are dependent on this type of excitations. Any upcoming study that investigates the seismic behavior of onesided rocking façades should include as an input the dynamic response of the building interacting with the façade. This building response could be modelled by means of the input motion of the return walls, because it has been shown that the input of the return walls has a detrimental effect on the façade response. Similarly, when studying the behavior of parapets located at the top of multistorey URM buildings, the response of the building needs to be included as widely demonstrated by existing literature. Exploring the introduction of algorithms such as the Hilber-Hughes-Taylor, that control the high frequency noise more efficiently, is recommended.

A C K N O W L E D G M E N T S
This project was (partially) supported by QuakeCoRE, a New Zealand Tertiary Education Commission-funded Centre. This is QuakeCoRE publication number 708. The authors also gratefully acknowledge the software 3DEC provided by Itasca consulting group under the Itasca Educational Partnership program. Dr Omar Al Shawa is gratefully acknowledged for kindly sharing the raw data to perform the simulations reported herein.
Open Access Funding provided by Universita degli Studi di Roma La Sapienza within the CRUI-CARE Agreement

D ATA AVA I L A B I L I T Y S TAT E M E N T
Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.