Seismic behaviour of rubble masonry: Shake table test and numerical modelling

The destruction of Amatrice and the surrounding villages in Central Italy after the 2016 seismic sequence was so impressive that engineers, authorities and local communities started sharing the common feeling that historical stone masonry buildings were too below current safety standards. The severe damage caused by the earthquakes led to a general distrust of traditional building techniques, leading to the conclusion that there is nothing to do but demolish and rebuild, perhaps with a false antique. Is there an alternative? Is there a way to combine safety and preservation of architectural heritage? This paper aims contributing to the understanding of the seismic behaviour of stone masonry by reproducing, through simulation on a shake table, the progressive loss of compactness of a real scale rubble masonry wall up to the ruinous collapse with the separation between the two external leaves. The laboratory simulation allowed to evaluate the decrease of the fundamental frequency with increasing damage and estimate the maximum displacement profile and the amount of cracking that the wall is able to sustain before failing. Eventually, two modelling strategies based on finite and discrete element methods were proposed and applied to verify the capability of simulating the out‐of‐plane seismic response and the failure mechanisms of rubble masonry.

F I G U R E 1 Wall sections and seismic damage after the 2016 Central Italy sequence adequate safety level. Inspections and damage surveys ( Figure 1) after the earthquake sequence showed that the poor performance was not a matter of insufficient in-plane shear or compressive strength. In some cases, the collapse occurred due to out-of-plane rigid-body overturning of the external walls lacking connections to the rest of the building, but in most cases the collapse was triggered by the loss of wall compactness with the external face detaching from the rest of the wall. Such fragmentation of stone masonry was observed also in other countries, such as Greece, 1 Iran, 2 Nepal, 3 New Zealand 4 and Pakistan. 5 Although collapse due to overturning has received considerable attention from the scientific community, which has defined reliable procedures currently incorporated in standards, 6 damage due to fragmentation has received little investigation in the literature. Nevertheless, failure due to separation of the two leaves of the wall precedes the onset of oscillatory motion and should therefore be considered as a prominent condition in the assessment procedure. The proneness of a rubble masonry wall to fail by fragmentation is difficult to assess, because it depends on the size, shape and bond pattern of stone units, as well as on mortar strength. Such features cannot be investigated from the outside only. Even the common flat-jack test fails to provide useful information about masonry monolithism. Collapse by wall disintegration began to be studied by Binda et al. 7 and references therein, who brought the vast range of masonry types and stone unit arrangements to the attention of the scientific community. A further contribution was made by Giuffrè 8 by means of a straightforward interpretation of masonry as a set of discrete blocks with frictional joints, while previously, masonry was modelled as a continuum homogeneous isotropic medium. A great effort was provided by Borri et al. 9 to establish the rules for evaluating the quality of stonework by means of an expert judgement based on the quality, size and shape of masonry units and geometric pattern of mortar joints. Nonetheless, because intervention techniques to prevent masonry disintegration are invasive and expensive, the decision to strengthen the masonry shall be based on quantitative assessments that are both accurate and feasible in current practice. An attempt to provide a mechanically based model of irregular masonry cross-sections was proposed by de Felice, 10 providing the capacity curve of wall sections surveyed in situ. However, a comprehensive and in-depth analysis of the failure due to fragmentation is still lacking.
The presence of stonework with poorly connected leaves has inspired several shake-table test programmes, emphasising the role of mortar as a bonding element between stone units. Meyer et al., 11 Elmenshawi et al. 12 and Giaretton et al. 13 studied single walls, with wall thickness ranging from 0.27 to 0.50 m and mortar being either absent or possessing up to a 2.6 MPa compressive strength. Costa et al. 14,15 and Candeias et al. 16 tested two C-plan wall assemblies without roof, with wall thickness ranging from 0.50 to 0.65 m and mortar compressive strength of 1.3 MPa. Sathiparan et al., 17 Ali et al. 18 and Wang et al. 19 tested single-room, single-story models with geometric scale ranging between 1:4 and 1:2 and mortar being either absent or up to a 4.0 MPa strong in compression. Benedetti et al., 20 Dolce et al., 21 Mazzon et al., 22 Magenes et al. 23 and Vintzileou et al. 24 investigated single-room, two-stories models with geometric scale ranging between 1:2 to approximately natural scale and mortar compressive strength between 0.7 and 4.6 MPa. Mendes et al. 25 and Senaldi et al. 26 tested multiple-stories models with geometric scale being 1:3 and 1:2, respectively, while mortar compressive strength was 2.5 and 1.8 MPa. Reconnaissance of previous experimental works showed that masonry fragmentation can only be replicated F I G U R E 2 Wall fair face of rubble masonry stonework in Central Italy Apennines in the laboratory if a number of conditions are met. The stone units shall be undressed, significantly smaller than wall thickness, arranged in a three-leaves pattern. The mortar must be weak, with mean compressive strength preferably below 1.5 MPa. The sample scale must be natural, otherwise even a poor mortar may be sufficient to guarantee monolithism. The upper boundary conditions should avoid the vertical arch effect and have a horizontal stiffness comparable to historical floors.
Aiming at complying with the above items, an experimental campaign on a shake table is presented in this paper to provide a clear understanding of the progressive lack of bond between the two leaves of stonework up to complete separation and to masonry failure under seismic action. In order to replicate the construction details typical in Central Italy, the specimen was built with stone units collected on site and with a very poor mortar, designed to simulate that of existing buildings.
In the second section, the experimental programme is presented, describing the construction of the prototype, the mechanical characterisation of the materials, the experimental setup, the sensors adopted for data acquisition and the ground motions selected to replicate the seismic sequence. In the third section, the experimental results are illustrated in terms of displacement profiles, progressive separation between the two leaves and overall dynamic behaviour, measuring the decrease in natural frequency with increasing damage. Finally, in the fourth section, two modelling strategies are presented, both based on the discretisation of individual blocks interacting with nonlinear joints, to assess the effectiveness and robustness of computational tools to predict the seismic failure of rubble stone masonry. As a further step towards the development of appropriate mitigation measures to prevent the failure by fragmentation and ensure the seismic protection of rubble masonry, the Reader is referred to a companion paper. 27 2 EXPERIMENTAL PROGRAMME

Definition of the masonry wall type
As a first step in the design of the prototype, a preliminary survey of typical features of masonry walls in the municipalities of Amatrice and Accumoli and in the surrounding villages was carried out. The severe damage caused by the 2016 seismic sequence allowed for a detailed survey of stonework from both geometric and mechanical points of view, including typical dimensions of stone units, joint pattern, height and thickness of the walls, aspect ratio and boundary conditions. A total of 21 wall sections were surveyed: the average wall thickness is about 0.55 m and height to thickness ratio of about 7.5. The external masonry leaves systematically lack in connection because the stone units are generally small and undressed, with the internal nucleus made of small size, irregular and poorly cemented units ( Figure 1). As for bond pattern, horizontal joints are rarely continuous along the façade and vertical joints are mostly not regularly staggered ( Figure 2). The wall specimen was built with a thickness of 0.50 m, a height of 3.73 m and a width of 1.63 m ( Figure 3A and B), in agreement with common values surveyed in the area. The bond pattern simulates a typical rubble stone fair-face wall, with two external leaves, no bond stones and a central nucleus roughly filled with smaller elements. An ad hoc mortar was designed to reproduce that encountered on site during the post-earthquake inspections. The back of the wall was finished with plaster, to simulate the internal side of a residential building wall. The testing setup is meant to mimic a vertical spanning wall fixed at the base on a foundation beam and restrained at the top by a steel frame simulating a floor connection. Aiming at releasing vertical displacements while preventing horizontal ones, the top restraining steel beams are equipped with low-friction rubber rollers, as it is shown in Figure 4 and described in detail later in the paper.

F I G U R E 4 Test setup
The wall is built on a reinforced-concrete foundation beam fixed to the shake table by means of two steel profiles and bolts. The upper part of the wall is reinforced with glass fibre reinforced mesh (66 mm × 66 mm) placed in the last three bed joints for a height of 260 mm, using a natural hydraulic lime mortar with improved mechanical characteristics (f m = 6.20 MPa) provided with three vertical steel bars with diameter 20 mm ( Figure 3D), grouted for about 0.36 m depth in the wall, to prevent localised damage where the rollers constrain the lateral top displacement of the specimen. The total weight of the wall specimen is about 85 kN.

Mechanical characterisation of construction materials
During the on-site survey, mortar and stone samples were collected and subjected to experimental tests to determine physical, chemical and mechanical properties. An invaluable contribution to the significance of this project came from the small village of Collespada (part of the municipality of Accumoli), because its population donated the stone units belonging to  buildings collapsed during the 2016 sequence. These units were used for the specimen tested on the shake table. From  a lithological point of view, stone units are made of compact marlstone characterised by unit  . Cyclic direct shear tests performed on the fraction of non-cohesive mortar samples with particle-size less than or equal to 2 mm, carried out at different normal-stress levels, provide an average friction coefficient equal to 0.9 and a cohesion equal to 0.025 MPa. 28 The mortar used for the construction of the wall was prepared according to the defined mix design, specimens were sampled during construction and tested according to EN 1015-11 the day before the shake table experiments. The samples collected have unit weight w m = 16.9 kN/m 3 (CV = 3.3%, six tests), compressive strength f m = 0.87 MPa (CV = 55.9%, 10 tests) and flexural tensile strength f mf = 0.36 MPa (CV = 49.9%, four tests). It is worth pointing out a marked dependency of the mortar mechanical properties from the curing time, which varied because approximately ¾ of the specimen was built before CoViD-19 lockdown and the upper ¼ after it. In detail, the average compressive strength for the six samples tested after 4 months of curing is 1.24 MPa, while it is 0.31 MPa for the four specimens tested 1 month after being cast. Similarly, the average value of f mf calculated on the basis of three tests performed after 4 months of curing is 0.45 MPa and it decreases to 0.09 MPa when referring to the sole test performed after 30 days of curing. Experimental tests were also performed on premixed natural hydraulic lime plaster specimens, determining the following mechanical properties: w p = 14.7 kN/m 3 (three tests), f p = 5.55 (CV = 7.5%, six tests) and f pf = 2.79 MPa (three tests).

Test setup
The tested specimen is restrained at the top by a steel frame fixed to the shake table. The steel structure is made of two braced frames, at the two sides of the specimen, and two girder beams, bolted on top of the frames, on the front and the back side of the wall ( Figure 4). Rubber rollers, fixed to the girders, are tightened to the top of the wall, to restrain out-ofplane displacements while releasing vertical displacement and rotation. The setup has proven to be effective to simulate floor restraint 27 and is designed to allow the relative displacements between top beam and foundation of the wall as can be expected in historical constructions. The dead load of a typical roof was simulated by 15 steel plates (total weight of 15 kN) placed over the top of the wall and connected by vertical bars to the top beam. The plates and the top beam are secured to the crane bridge by fixing eyebolts at the end of each bar and hanging them at loose wire ropes to prevent their collapse at high intensity shakings, when the wall is expected to fail ruinously below.  Table 1, sorted by their date as they were ordered in the test sequence. For each signal, information is summarised in Table 1: date and moment magnitude (M W ) of the event, name and label of the record station, orientation of the horizontal component, peak ground acceleration (PGA), velocity (PGV) and displacement (PGD) in horizontal and vertical directions. The pseudo-acceleration response spectra of the strong motion records applied to the shake table are shown in Figure 5, together with the response spectra of the signals recorded on the foundation of the specimen during the tests carried out with SF = 0.6. Applied and recorded spectra are in very good agreement, except in the highest frequencies range. In addition to strong motion tests, white noise (WHN) tests AMT, Amatrice; CNE, Castel Sant'Angelo sul Nera; NRC, Norcia; PGA, peak ground acceleration; PGD, peak ground displacement; PGV, peak ground velocity.

F I G U R E 5
Pseudo-acceleration response spectra of the horizontal (left) and vertical (right) components of the input signals and of those recorded on the reinforced concrete foundation at scale factor equal to 0.6 with nominal PGA = 0.05 g were performed before each sequence with the same SF, to better investigate the evolution of the dynamic properties of the wall as the damage increased.

Instrumentation
Tests are carried out on a 4 × 4 m 2 , six degrees of freedom shake table. Four horizontal and four vertical hydraulic actuators govern the motion of the table, which has a payload capacity of 300 kN, and can reach ±3 g acceleration and a maximum displacement component equal to ±125 mm. A triaxial piezoelectric accelerometer with <2% accuracy and 0.01 g resolution is fixed to the reinforced-concrete foundation beam to record the three components of the acceleration of the high-resolution 3DVision motion system 29 is used to sample the displacements at pre-determined locations, by means of retro-reflecting spherical markers distributed on a (approximately) 40 × 40 cm 2 grid on both sides of the wall specimen. Additional markers are placed on the foundation, on the top beam and on the steel structure, for a total of 38+38 markers, as shown in Figure 4. Acquisition is performed at 200 Hz with nine near infrared digital cameras, placed around the shake table. Being a non-contact one, this system allows for displacement tracking eliminating some of the main drawbacks of traditional displacement monitoring system, as the construction of an instrumentation support setup or the risk of damage to the instruments during the test. Measurements on the wall are carried out for 20 s after the end of the input, to record displacements under free oscillations and perform dynamic identification.

Progressive damage and leaf separation
In Table 2, main information and results are collected, namely a sequential test #, the record and the SF. This information is also summarised in the test label. The maximum and minimum accelerations recorded on the reinforced concrete foundation beam are named as a h,max and a h,min for the horizontal component and a v,max and a v,min for the vertical component. Moreover, d A denotes the maximum (forward) out-of-plane displacement of the front side of the wall, calculated from the average of the horizontal displacements of the four markers A51, A52, A53 and A54 at z = 2.53 m height (see Figure 4), which exhibited the maximum displacements for SF = 0.6, 0.8. Similarly, d B is the minimum (backward) out-of-plane displacement of markers B51, B52, B53 and B54 on the back side (at the same height). Finally, f is the fundamental frequency of the wall, determined by means of multi-input/multi-output (MIMO) experimental modal analysis. The sequences with SF = 0.2, 0.4 and 0.6 were completed for all three records, whilst only NCR and CNE were applied in the last sequence with SF = 0.8. The wall exhibited no damage in the first two test sequences with SF = 0.2 and 0.4, during which the maximum absolute base acceleration recorded on the foundation beam was a max = max{a h,max ; − a h,min } = 0.26 g. After AMT02, the steel plates and the bars supporting the rollers were stiffened, as the behaviour of the top restraint was considered inappropriately deformable. The first signs of damage were observed during NRC06 F I G U R E 6 Progressive damage development (eye observed) (a h,max = 0.22 g) and included the expulsion of dust and small stone fragments from the front side of the wall, near the top left corner and the development of a vertical crack in the middle of the left side, revealing the onset of a local separation between the two wall leaves ( Figure 6A). Thin horizontal cracks were also detected on the plaster on the back. It is not surprising that the crack pattern developed even if the maximum acceleration did not monotonically increase, because damage progressively accumulated test after test. Moreover, the maximum base acceleration does not necessarily fully represent the intensity of an earthquake in terms of structural demand. Going forward, the cracks widened during CNE06 on both wall sides. Those on the front side extended and, at this stage, involved most of its height; furthermore, some stone units near the top left corner appeared detached from the surrounding ones ( Figure 6B).
The width and length of the vertical crack on the left side also increased, indicating the progress of leaf separation. A first partial collapse occurred during AMT06 (a h,max = 0.45 g), involving the stone units of the front leaf from about 1 m height to the top ( Figure 6C). Further stone units fell down from the same corner during NRC08 (a h,max = 0.29 g), whereas, despite a widening of the cracks on the plaster, no failure took place on the back side ( Figure 6D). The wall collapsed completely during CNE08 (a h,max = 0.42 g). The front leaf separated from the back one, suffered the formation of a horizontal crack along its entire width at about 3/4 of the height, and then collapsed with a two-block out-of-plane bending mechanism, associated with the fragmentation of the masonry ( Figure 7A-C). As a consequence, the back leaf remained isolated as a slender single-leaf wall, and, eventually, collapsed as well ( Figure 7D-F). The progressive separation of the two leaves of the wall was monitored by calculating the relative displacement (δ) between the markers near the corners on the front side and the corresponding ones on the back side. The relative displacement between three couples of markers on the left corner, namely A41−B41, A51−B51 and A61−B61 (A), as well as between three on the right corner, namely A44−B44, A54−B54 and A64−B64 (B) is portrayed in Figure 8. Positive values of δ indicate that the corresponding markers were moving apart, whereas negative values indicate that they were moving closer. The values recorded during all the strong motion tests are plotted in sequence to identify the leaf separation occurring during each test and its residual value progressively developing test after test. The first three strong motion tests with SF = 0.2 are omitted because no leaf separation took place and, for this reason, the time (t) on the x-axis starts at 120 s. The curves show the occurrence of a non-negligible leaf separation from NRC06 on the left side (residual displacements of 2-5 mm) and from CNE06 on the right side (2-6 mm), progressively increasing in the following tests. Large jumps of δ are associated with crack occurrence (due to the brittle nature of the mortar), whereas final values at test end are associated with residual crack width. Such residual values are larger on the left (2-5 mm after NRC06, 13-26 mm after CNE06, 35 mm after AMT06 and 58 mm after NRC08) than on the right (3-6 mm after CNE06, 12-22 mm after AMT06 and 28-40 mm after NRC08). The largest values of δ were recorded at z = 2.09 m (A41−B41 and A44−B44), consistently with the crack pattern displayed in Figure 6. Note that for A41−B41 and A51−B51 data are unavailable after AMT06 due to the collapse of the stone units from the corner.
The profiles of the out-of-plane displacements of the wall (d) were derived starting from the data of the markers placed on the front and back sides, averaging those at the same height as shown in Figure 9 for all the strong motion tests. They are calculated in the time instants when d A and d B were attained, respectively (each profile provides a picture of one specific time instant, those of the front and of the back do not coincide). In the first tests, when no significant damage developed, the profiles are basically linear with a rotation at the base of the wall.
On top, as expected, the restraining system prevented the out-of-plane overturning of the wall but allowed for some horizontal displacements (vertical displacements and rotations, instead, were free). Such top horizontal displacements increased from 3-4 mm in the first tests (SF = 0.2), to 7-9 mm (SF = 0.4) and to 8-13 mm (SF = 0.6, 0.8), with negligible differences between front and back sides. From AMT06, the onset of a failure mechanism is clearly detectable, as the profiles exhibit their larger displacements at mid-height, where the horizontal cracks formed. Masonry fragmentation took place after the attainment of d A and d B and, therefore, did not affect the portrayed displacement profiles until NRC08 whereas, for CNE08, a time instant was selected before collapse. Finally, it is worth noting that, especially in the last three tests, the displacements of the front side were larger than those of the back one, due to the plaster layer applied to the latter.

Seismic capacity and acceleration versus displacement behaviour
The acceleration versus displacement curves are plotted in Figure 10A to represent the global behaviour of the wall under earthquake inputs. 15,30,31 The displacements d A and the accelerations a h,max are considered for the positive portions of the curves and are associated with the forward movement of the front side, whereas d B and a h,min are taken for the negative parts and correspond to the backward movement of the back side. Three curves are plotted, one for each record (NRC, CNE, AMT), to highlight the sensitivity of the response to the earthquake base input. 32 For the first two sequences (SF = 0.2, 0.4), the relationship between base acceleration and out-of-plane displacement is basically linear and the three curves display a similar slope, consistently with the observed damage development ( Figure 6) and the displacement profiles ( Figure 9) described before. Then, the development of cracks is associated with a decrease of stiffness, such that the overall behaviour is nearly bilinear. Finally, and as expected, the displacements recorded at collapse are much larger than all the previous ones (the last value recorded before the marker fell on the table was considered for the plots, when necessary). As for the accelerations, instead, that attained in AMT06 (a h, max = 0.45 g), when the first partial collapses occurred, slightly overcomes that recorded in the last test (CNE08, a h,max = 0.42 g), and represents the seismic capacity of the wall in terms of peak base acceleration. It is worth noting that the response of the wall in the positive direction under AMT signal does not show the same reduction of stiffness at SF = 0.6 exhibited under NRC and CNE as well as in the negative direction. This behaviour may be due to the activation of a second mode of vibration, in addition to the first one, which governed the dynamic behaviour of the specimen in all the other tests.
The seismic performance exhibited by the specimen along the shake table investigation is not independent from the sequence of dynamic tests, due to the progressive damage accumulation. In order to derive an estimate of damage, the energy dissipated in the hysteretic cycles E hyst was calculated for each strong motion test starting from the acceleration and displacement time-histories recorded at each marker, and considering their tributary mass. 31,33,34 The result is shown in Figure 10B, and indicates that dissipated energy was relatively negligible in the first two test sequences (with SF = 0.2 and 0.4), grew in those with SF = 0.6 and exhibited a significant increase in the last two tests with SF = 0.8.

Dynamic properties
The dynamic amplification (α) is analysed by means of the ratios between the maximum absolute acceleration recorded by the markers on the wall and that of the markers on the foundation. Until CNE06, the plots show a nearly linear increase of α with the height, and vary between 1.5 and 2.5 in the lower half of the wall and between 2 and 4 in the upper one, even if no clear relationships between dynamic amplification and intensity of the seismic input can be recognised. In the last three tests, instead, there is a clear deviation from the linear profile towards a vertical bending profile. The fundamental frequency of the wall is determined with an MIMO approach, as follows. First, the accelerations of the four markers on the foundation are considered as seismic inputs and are calculated by double derivation of their displacements. As for the output, the accelerations obtained from 48 markers on the wall, comprised between z = 1.28 m (row x2) to z = 3.44 m (row x7) and on both sides, are taken. No filters are applied to these accelerations. Then, the transfer function (Ω(f)) is calculated from each input marker to each output marker. Finally, all the transfer functions are averaged. The numerous measurement points available with the 3DVision system made it possible to average many transfer functions and derive a global dynamic characterisation of the wall that is not affected by possible local singularities of the response or by errors in acquired data (e.g. due to the presence of dust or the falling of a stone unit).
The transfer functions calculated in WHN tests are plotted in Figure 11A. The fundamental frequency f is determined as that corresponding to the main peak of Ω(f). The frequencies for all the shake table tests are shown in Figure 11B. At the beginning, for the undamaged wall, f resulted equal to 5.20 Hz (WHNa) and decreased to f = 1.20 Hz for CNE08. It is worth making some remarks on the variation of f. First, it results higher in WHN tests than in strong motion ones, due to the smaller acceleration and displacement amplitudes. This behaviour is consistent with that of a non-linear system whose frequency is higher for small oscillations. Second, the increase of the frequency after AMT02 (compare WHNb, f = 5.86 Hz, with WHNa, f = 5.20 Hz) was due to the modification of the restraint system. Finally, the non-linearity of the structural response and the different amplitude and energy content of the seismic inputs, actually affected the values of f, which even increased in some cases from one test to the following one. The overall trend is however well represented by the values of f averaged over each set of tests with the same SF (4.6, 4.7, 3.2 and 1.7 Hz), which starts decreasing after the end of the second sequence (with SF = 0.4), consistently with the trend shown by the dissipated energy in Figure 10B.

NUMERICAL MODELLING
As shown by the experimental results collected, the seismic behaviour of stonework can be strongly affected by the geometry of the stone units, their mutual arrangement and the bonding at the stone-mortar interface. The numerical approach to be adopted must therefore rely upon a micro-modelling technique, 35 which accounts separately for the stone units and the contacts, allowing the cracking of the latter, as well as the sliding and the complete separation of the units. Several computational strategies are theoretically eligible, but only some were actually used for the problem at hand. Lattice discrete particle model was used to simulate the in-plane static behavior of irregular stonework accounting for the interaction among coarse material heterogeneities. 36 Rigid units and nonlinear springs were used for modelling out-of-plane cracking of masonry, resorting to mass-proportional viscous damping, hysteretic energy dissipation and an implicit integration scheme. Although able to account for frictional-torsional wall interlocking and interaction between masonry leaves, 37 this approach was not utilised so far to account for full disintegration. Despite the investigated brickwork was not prone to disintegration, the applied element method was exploited to simulate whole-building shake-table tests. 38 Physical and hysteretic damping were implemented, within an implicit time-integration scheme. The non-smooth contact dynamics method was applied to account for masonry fragmentation of monumental constructions, 39 resorting to an implicit integration scheme and to friction alone for energy dissipation. In this section, two modelling approaches are proposed to simulate the dynamic behaviour of the wall both a priori, that is before the tests, and a posteriori, relying on a wealth of information collected during the experiments and on a more detailed description of geometry and boundary conditions. The first approach is based on a mixed finite element method-discrete element method (FEM-DEM) technique and considers elastic units interacting through nonlinear springs, resorting to an explicit integration scheme with stiffness-proportional viscous damping. It was originally developed in the context of short-duration dynamic problems, such as impacts and blasts. 40,41 The structure is discretised into blocks, representing units expanded to include half of the thickness of the mortar joints. Blocks are, in turn, meshed into finite elements and separated from each other by contact interfaces. Finite elements can be equipped with linear-elastic or non-linear constitutive laws, with blocks potentially separating in different fragments depending on the expected behaviour. Interfaces allow for cracking, separation and formation of new contacts. Contact stiffness response is governed by a penalty approach, hence there are no additional springs but interface properties depend on units geometric and material properties, as well as on a non-dimensional SF. Interfaces may fail either in tension or in shear. After shear failure, the interaction among the blocks is governed by friction, whereas after tensile failure no interaction in tension is possible. The equations governing the problem are solved according to an explicit strategy, which is well suited to analyse crack propagation and dynamic response. In the case at hand, the FEM-DEM model was implemented in LS-DYNA within an ANSYS environment. 42 The second approach is based on a classic DEM formulation, first conceived in the framework of rock mechanics. 43 It was efficiently implemented in the software UDEC and 3DEC 44 and widely adopted for the analysis of masonry structures under seismic or static loading. 45 The solution algorithm is based on the explicit integration of the equation of motion at each block centroid. Discrete blocks are modelled as rigid bodies separated by zero-thickness joints; they interact with each other by means of contact points located (in pairs) at each interface, in which elastic and strength properties of the assemblage are lumped and described by two nonlinear springs, reproducing the normal and the tangential response, respectively.
Although the adopted procedures were already used to study the behaviour of masonry structures subjected to horizontal actions, the novelty introduced in this paper relies in a systematic comparison between the two approaches, both stemming from the same constitutive parameters and the same input actions. Moreover, when dealing with DEM method, the influence of the discretisation adopted could be appraised thanks to a significant refinement of the mesh. This refinement represents a key aspect of the numerical study proposed as, unless in case of partial collapse occurrence, the internal structure of rubble stone masonry walls remains unknown and/or uncertain.

A priori numerical predictions
Stemming from the same set of information, the expected seismic behaviour of the wall was first blindly predicted with both models before the execution of the tests. At this stage, a preliminary ambient vibration identification test had already been performed in the laboratory, which provided a first natural frequency of 3.7 Hz for the wall unrestrained at the top. The geometry of the 2D model was designed to reproduce in a relatively rough way the effective pattern of the specimen while mechanical parameters were assigned on the basis of literature values (strength) and dynamic identification (stiffness). In both numerical models, failure is expected to occur at the joints, thus units are considered either elastic in the FEM-DEM model (linear elastic, isotropic, homogeneous behaviour with Young's modulus and Poisson's coefficient equal to 1400 MPa and 0.20, respectively) or rigid in the DEM model, while nonlinearities are lumped in the joints. The assigned stiffness of the joints in the DEM model are k n = 8.2 GPa/m and k s = 4.1 GPa/m for normal and shear spring, respectively, targeting the measured first natural frequency of the wall.
Albeit with different accuracy, both discretisations account for the presence of two external leaves, for which block density was assumed equal to 1948 kg/m 3 and an inner core, made of smaller stone units in which a 15% lower value of unit mass is considered to account for the presence of voids. In the case of the FEM-DEM mesh the bond pattern is crudely simulated with two external leaves made of rectangular units, of varying height and width, and a nucleus made of smaller square units. 46 The DEM geometry is more refined, but still assumed to be representative of a tentative geometry, ideally assigned to reproduce a 'typology' rather than a surveyed section. The reinforced-concrete foundation, the top beam and the steel plates were explicitly modelled with their actual size and mass. Joint behaviour was described by means of a Mohr-Coulomb strength criterion with friction angle ϕ = 35 • , cohesion c = 0.06 MPa and tensile strength f t = 0.05 MPa either in the FEM-DEM or in the DEM model. In the FEM-DEM model, a 10% stiffness-proportional damping ratio was introduced, as also suggested by Meyer et al., 11 while a 2% mass-proportional damping at 3.85 Hz was considered in the DEM model.
After the activation of self weight, the analyses were performed subjecting the wall to the same seismic signal to be assigned to the shake table, conveniently cut to a 10 s duration. Each run, characterised by a progressively increasing SF, was performed after 2 s of free vibrations, starting from the deformed/damaged configuration at the end of the previous run. The FEM-DEM mesh is formed by 206 blocks, with a total of 6200 finite elements and 27,334 nodes, while the DEM model is made up of 151 blocks. Computer time for the FEM-DEM analyses is about 2.5 that of DEM simulations.  Figure 12A, D and G, respectively. These pictograms also report the location of A and B points on the front and back sides, respectively, whose displacement was monitored either experimentally (these points correspond to the location of two markers on the left, A51 and B51, and right sides, A54 and B54) and numerically (two points at the same height of the fifth markers row), calculating the relative value plotted on the right ( Figure 12C, F and I). The progressive leaf separation (δ) is plotted starting from SF = 0.4 up to failure. This choice is motivated by the fact that prior to that series of signals the behaviour of the wall remained essentially in the elastic range. The experimental response, as commented in the previous section, exhibits a different evolution of the crack pattern on the two sides of the wall, reflected in the fact that, looking at the first column of markers, collapse is attained at the AMT06 signal, while referring to the fourth column of markers, closer to the right side of the wall, a higher seismic capacity is attained prior to failure (occurred during the NRC08 event).
It can be observed that during the experiments, the two faces of the wall experienced very little relative displacement during the first series of shakings, with some pounding and inner core crushing during NRC04 and CNE04 tests, leading to some oscillation of the δ time history as reported in Figure 12C. The FEM-DEM model is not significantly sensitive to these small excitations and this behaviour might be ascribed to the relatively regular pattern that characterises the mesh, to which also a prevalence of shear sliding mechanism on the horizontal planes could be related. The DEM model, on the contrary, despite a shift to the AMT04 signal, shows comparable small δ perturbations, most probably because of the more irregular discretisation. Looking at the time histories reported in Figures 12F and I it is possible to observe how the two F I G U R E 1 3 Discrete element method (DEM) geometric micro-modelling discretisation numerical models essentially agree in terms of relative cumulative displacement under the CNE06 signal, which is also in line with the experimental evidence on the right side of the specimen. Differently from the FEM-DEM model, which shows an additional capacity and fails under NRC08 signal, the DEM model rapidly attains collapse during the AMT06 test, and this resembles what experienced on the left corner of the wall. The overall behaviour depicted by the results of the FEM-DEM simulations, however, seems to be consistent with the A54-B54 relative displacement time history. Figure 12B, E and H report the failure configuration, that is the last converged step in case of numerical simulations, which are both characterised by separation of the two leaves and, in case of the FEM-DEM simulations, by a more marked out-of-plane flexural mechanism involving the external front leaf.
From an engineering perspective, both simulations, performed with a relatively poor set of input information, provided a realistic representation of the behaviour of the wall. The significance of these analyses is two-fold: on the one hand, predictive simulations are useful, in the context of experimental tests, to set convenient testing programmes and seismic actions to be applied. On the other hand, in real engineering practice problems, it is also common to rely on few mechanical parameters and to have a rough or at worst completely unknown internal structure representation, which could only be approximated resorting to standard section typologies. In this sense, both the models proved reliable in providing satisfactory results in terms of seismic capacity and an overall good agreement in terms of damage accumulation (here expressed through leaves detaching).

A posteriori numerical simulations
More detailed numerical analyses were carried out after the experimental tests. Both models were updated in terms of mechanical properties. However, in the FEM-DEM model, the block geometry was kept the same of a priori simulations to avoid too small FEs that would involve a very short integration time step. On the contrary, in the DEM simulations, the mesh was refined to reproduce the actual geometry of the blocks and their arrangement. 10,47 To this end, both lateral sections of the wall were modelled,to deliver a representation of the two lateral sides of the real structure and, at the same time, to highlight the geometry-related effects for the same loading and boundary conditions. The process of definition of the geometry for the left (A) and right (E) side of the specimen is presented in Figure 13. In detail, stone units were first identified on the front side of the wall (C), then labelled with a number denoting the course and ascending from bottom to top and a letter indicating whether a unit is placed on the left (L) or right (R) side. When  Figure 4 and adopted throughout the previous sections.
In both models, the stiffness and the block density are left unchanged. Conversely, strength parameters were increased (cohesion c is now equal to 0.18 MPa and tensile strength f t = 0.15 MPa). This choice stems from the mechanical characterisation performed on mortar samples taken during construction. The seismic input is modelled assigning to the foundation and top beam the velocity time histories recorded on the shake table (acceleration input is adopted in case of FEM-DEM simulations). Accounting for the seismic signals experimentally recorded during the tests allows to implicitly consider the real boundary conditions at the top beam, which could play an important role in the response of the wall and the slight variation between the input signal assigned to the table and the recorded ones, as reported in Figure 5.
These actions are obtained by deriving and then averaging displacement histories recorded by markers AC1-BC1 (left side) and AC4-BC4 (right side) for the foundation and markers AU1-BU1 and AU4-BU4 for the top beam. The derivation process is carried out applying the Savitzky-Golay smoothing filter for data noise reduction. 48 The same seismic signal, in terms of accelerations, adopted for the right side of the wall in the DEM simulations is adopted for the FEM-DEM analyses. The time necessary to perform the a posteriori simulations is increased by about 15% compared to a priori simulations.
A complete representation of the numerical-experimental matching for all the performed tests is shown in Figure 14 in terms of displacement profiles. Differently from the profiles plotted in Figure 9, the ones reported here are not the result of an average of the four markers at the same height, but are that recorded on the left side of the specimen, to be more immediately correlated with the DEM simulation. The most remarkable aspect emerging from this comparison is that while DEM simulations are in fairly good agreement with the experiments, the FEM-DEM model shows an overall lower value of the out-of-plane displacement, related to the fact that the model tends to exhibit a quasi-elastic behaviour during the first seismic sequences, undergoing some shear sliding on the horizontal joints and then a sudden and brittle loss of equilibrium with out-of-plane failure mechanism under NRC08 seismic signal. Figure 15 summarises the most representative results in terms of displacement on the front (A marker) and back (B marker) side of the wall at the sixth marker row height. The experimental time histories, plotted in gray, are compared with those of two selected points on the FEM-DEM (Figures 15A and D) for the NRC06 and AMT06 signals, respectively, while, having performed two different DEM analyses for the left and right sides, two plots for each signal are reported in this case ( Figure 15B, C, E and F). The results are plotted for the SF = 0.6 seismic sequence, the last completed for all three records. It can be observed that the FEM-DEM model under the NRC06 signal ( Figure 15A) provides an overestimation of the displacement while correctly replicating the absence of separation between the two leaves. In the DEM model ( Figure 15B), a rigid sliding of both leaves without mutual separation is responsible for the shift of the numerical curves in the negative portion of the displacement graph. While for the right side, apart from the shift, the DEM model correctly reproduces the experimental time histories, in the case of the left side ( Figure 15C), the model does not replicate the separation that is instead experimentally measured in the second half of this record. Looking at the AMT06 graphs, the FEM-DEM model results show, differently from the experiment, a delayed trigger of the leaf separation process but, on F I G U R E 1 5 Numerical and experimental displacement histories at markers A51, B51 (left side) and A54, B54 (right side) average, the displacement trend on front and back side seems to be captured. The DEM numerical model captures the local failure mechanism occurring on the left portion of the wall front side, with just a limited delay ( Figure 15F). The right side simulation is in fair agreement with the experimental 3DVision measures too ( Figure 15E) though, on average, the numerical model provides a 40% higher estimate of the relative displacement between the two selected points.
The update of the models from a mechanical and (for the DEM model) from a geometrical point of view proved to be effective in improving the descriptive capability of the numerical simulations essentially in terms of profiles and, especially for seismic signals closer to collapse, in terms of leaf separation. With respect to this aspect, the wealth of information that can be extracted from DEM simulations is further compared to experimental outcomes.
The cracking process that the DEM model undergoes during the sequence is portrayed in Figure 16. The combination of thickness and shade changes marks the evolution from early crack formation to significant separation (>10 mm), which is in very good agreement with the experimental damage accumulation observed by eye ( Figure 6). On the left side, in fact, the first cracks are visible after the NRC06 test, propagating along the height of the wall and widening until the partial collapse of the front-left edge under the AMT06 test. On the right side, crack formation in the numerical model starts a bit earlier than what experimentally measured, but the final picture obtained from the simulations very faithfully reproduces the detected crack pattern at the end of the last completed test before collapse (NRC08).
In order to represent the different behaviour as observed on the left and right sides (and consistently reproduced via the DEM model), leaf separation profiles are plotted in Figures 17B and F, at the time instant of maximum relative displacement between front and back sides along the first and the fourth marker columns and along the wall sections of the model for left and right sides, respectively. Numerically obtained relative displacements are extracted along the height of the wall on n pairs of points on the opposite sides of the specimen (being n the maximum number of rows in the discretisation, i.e. 34 for the left profile and 33 for the right one) while the dashed red lines are obtained interpolating the relative displacements between the eight pairs of equally high markers, whose position is also reported in Figure 17A and E. On the left side, in spite of a slight underestimation of the actual leaf separation, the model is able to account for the formation of the failure mechanism, with the height of the central hinge almost coincident with the experimental one. On the right side, the localised increment of the relative displacement, for a height between 2.0 and 2.7 m, is due to the sliding of the blocks 18A through 22A ( Figure 17G). In general, along the height of the right side of the wall, until the localised failure, the numerical model shows larger relative displacements, as already anticipated in Figure 15F) for a specific pair F I G U R E 1 6 Progressive crack opening at the left and the right sides in the discrete element method (DEM) model of markers. This behaviour is coherent with the first cracks of the numerical model being visible much before than the first cracks of the specimen.
As a closing remark of this section, the laboratory conditions at the end of AMT06 and NRC08 tests for left and right sides, respectively ( Figure 18A and D), are compared with the failure mechanisms numerically obtained ( Figure 18B, C and E). The numerical collapse for the left side occurs under the AMT06 excitation and a compatible failure is also experimentally found. On the contrary, the numerical model of the right side fails during the NRC08 test because of a much more localised instability involving a small group of stone units (18A, 19A and 20A) without an experimental correspondence, though the crack pattern along the wall side is well reproduced by the model. In this case, the FEM-DEM model succeeded in capturing the separation between the two leaves.

CONCLUSIONS
Fragmentation by leaf separation is the most critical collapse mechanism of rubble stone masonry. With the aim to better understand the onset and evolution of masonry disintegration, a shake table test campaign was carried out on a masonry wall replicating the typical conditions of a historical building (real scale, medium to small size units, absence of bondstones and weak mortar). The specimen was built with stone units collected from buildings collapsed after the 2016 Central Italy earthquake, arranged in two separate leaves with inner nucleus and laid with poor quality mortar. The testing setup simulated a vertical spanning wall fixed at the base and horizontally restrained at the top by a steel frame simulating floor connection. The dead load of a typical roof was simulated by steel plates placed over the top of the wall. Three input signals recorded during the 2016 Central Italy sequence in three different locations, namely NRC, CNE and AMT, were applied to the shake table in both the horizontal and the vertical directions, with increasing amplitude SF (SF = 0.2, 0.4, 0.6 and 0.8).
The first damage was detected during the third sequence (SF = 0.6), NRC06 (PGA = 0.22 g), which caused a vertical crack in the middle of the left side with the onset of separation between the two leaves, and two horizontal cracks on the back side. AMT06 (PGA = 0.45 g) induced the failure of the left corner of the front side and the widening of the horizontal cracks on the back surface. Under NRC08 (PGA = 0.29 g) a progressive failure of the right edge of front leaf occurred, along with flexural mechanisms on the front and on the back leaves. Under CNE08 (PGA = 0.42 g) global collapse occurred as a two-block, out-of-plane flexural mechanism of each of the two leaves on its own side.
Although testing a single wall specimen does not allow definitive conclusions to be drawn, the experimental investigation provided useful information on the seismic behaviour of rubble masonry walls and succeeded in reproducing the actual failure by fragmentation that typically occurs in the case of severe earthquakes, including: -the progression of masonry leaf separation, from the initial cracking, to the subsequent debonding of stone units of the outer face, until the overall collapse; F I G U R E 1 8 Experimental and numerical damage patterns for NRC08 (right side A, B, C) AMT06 (left side D, E). AMT, Amatrice; NRC, Norcia -the intrinsic brittle response of the base-acceleration versus displacement curve up to the partial failure and the subsequent inelastic behaviour, due to the progressive separation between the two leaves; -the reduction in the natural frequency of the wall as the damage increases.
Two different approaches, based on a 2D model of the lateral section of the wall, resorting either to mixed FEM-DEM or DEM, were used to predict the response of the wall and to back-analyse its experimental behaviour.
A limited but consistent (among the two models) set of input parameters were adopted, namely calibrated stiffness properties and reasonably chosen strength values, which were then adjusted in the postdictions. The adopted geometric description aims at representing with reasonable accuracy the internal structure of the wall that however, in principle, is not known in detail when dealing with real full-scale engineering problems involving existing structures. The update of the geometric description that was adopted in the case of DEM analysis aimed at emphasising the effect of the bond pattern, showing a different response of the left and right sides of the wall, as also experimentally observed.
The performed numerical analyses highlighted: -the reliability of discrete models for a clear understanding of the mechanical behaviour and a good replication of both the seismic capacity and the failure mechanisms; -the strong influence of the geometry of masonry units in the out-of-plane capacity of rubble masonry walls enhanced by the weakness of the mortar, and the consequent suitability of discrete models over equivalent continuous ones.
Overall, the results provide a benchmark for further investigation and indicate possible directions for the development of assessment techniques to be used in engineering practice.

A C K N O W L E D G E M E N T S
This work was carried out within the Research Projects 'SICURA Sustainable technologies for the seismic protection of the cultural heritage' (2018-2020) and 'SISMI Technologies for the safety improvement and the reconstruction of historical centres in seismic prone areas' (2018-2019) funded by Lazio Region, and it was partially funded by the 'Dipartimento di Protezione Civile -Consorzio RELUIS'. The intermunicipal operative office of Accumoli and Amatrice is kindly acknowledged for allowing the collection of the stone units from the village of Collespada. The opinions expressed in this publication are those of the authors and are not necessarily endorsed by the funding bodies.
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
The data that support the findings of this study are available from the corresponding author upon reasonable request.