Fragility curves for light damage of clay masonry walls subjected to seismic vibrations

The probability of light damage to unprepared, unreinforced masonry structures exposed to induced seismicity due to gas extraction in the north of the Netherlands is still under investigation. Repeated light seismic excitations caused by frequent, light and nearby earthquakes have been linked to economical losses and societal unrest in particular, with extensive damage claims. Moreover, the damaging potential of the seismic events has been related to the condition of the structure, especially if damage corresponding to settlement causes is already present. A comprehensive testing campaign oriented towards the initiation and progression of light damage of replicated clay brick masonry has been conducted at Delft University of Technology. Based on these tests, calibrated finite element models have been produced. This article uses the calibrated non-linear time-history models to simulate the effect of earthquake ground motion on a variety of initial conditions, wall geometry, material properties, and number, type and intensity of earthquakes. The results are then used to regress a relationship between damage and these parameters. This is subsequently employed to run a MonteCarlo simulation and produce fragility curves where the probability of exceeding specific damage values for various initial damage levels is presented against the seismic hazard. The vulnerability or fragility curves show that visible damage, with cracks wider than 0.1 mm, appears, with a 10% exceedance probability, at 13 mm/s of peak ground velocity; but, if the masonry had already undergone some light, yet imperceptible damage, a PGV of 6 mm/s was sufficient to aggravate it into visible cracks. To attain a 1% probability of exceeding light damage however, for which the masonry would need more invasive repair, it was observed that PGVs larger than 15 mm/s were required. These fragility curves were finally compared to graphs from other authors and found to capture well the variability in the range assigned to light damage.


Introduction
Light earthquakes, induced by gas extraction activities in the north of the Netherlands (NAM et al. 2015;Crowley et al. 2019b), have triggered extensive studies into the safety and damage vulnerability of unreinforced masonry structures (van Elk et al. 2019). These structures, built throughout the region and in the rest of the country, were never prepared to be exposed to seismic activity and its resulting vibrations. Consequently, the unreinforced fired-clay brick masonry, ubiquitous to masonry veneers or outer cavity wall leaves and also present in older (double-wythe) structural walls, has exhibited damage, principally, in the form of cracks opening in the plane of the walls (Staalduinen et al. 2018). The earthquakes, frequent, but light in nature, with vibration levels in the order of 5 to 20 mm/s (Noorlandt et al. 2018), have led to mostly aesthetic or cosmetic damage, depicted by cracks which are considered herein to range between 0.1 and 5 mm in width. While wider cracks and more severe damage have been reported, these are outside the scope of this study, laying the focus instead, on lighter and more common expressions of damage due to processes acting in the plane of the walls. Out-of-plane failures are more commonly associated with the ultimate limit state (ULS) or near-collapse states of the structures, while light damage is linked to in-plane actions (Staalduinen et al. 2018).
To estimate the impact of seismic activity on a regional level and to predict its future consequences, it is necessary to quantify the probability of damage to buildings. Damage states, as defined by the EMS (Grünthal et al. 1998;Vent et al. 2011), categorise the condition of a structure on a scale of zero to five where DS0 corresponds to no damage, DS1 to aesthetic damage, and ultimately, DS5 to total collapse. Many studies have sought to characterise the behaviour of masonry structures for the states DS3 to DS5 as these are relevant for quantifying the safety of the structures (Esposito et al. , 2019Messali and Rots 2018;Tomassetti et al. 2017;Graziotti et al. 2017); yet, the aesthetic/cosmetic damage state (DS1), or light damage, more difficult to evaluate due to the increased influence of the variability in the structures and materials, and the interaction of other loading processes, has proven to be especially relevant in the context of induced seismicity where an operator is responsible for the damage with associated economic losses and social unrest (Vlek 2018).
Several studies have looked at the development of light damage in masonry structures (Mojsilović and Salmanpour 2016;Mignan et al. 2015;Didier et al. 2018; García-Macías and Ubertini 2019; Kita et al. 2020). For example, Mojsilović and Salmanpour 2016 looked at crack formation in masonry walls, and Mignan et al. 2015 considered light damage of masonry structures as part of their risk assessment. Similarly, Didier et al. 2018 quantified damage based on cracks on masonry panels subjected to in-plane actions. Nevertheless, the quantification of (light, cosmetic) damage and its physical relationship to earthquake vibrations, especially in a probabilistic manner, has yet to be explored. Some studies available [e.g. 33, 36] do incorporate DS1 or DS2 as a probabilistic outcome, but only as an extrapolation or inference from the computation of other damage states or limits related to structural components. Conversely, this paper aims to present fragility curves that determine the probability of light damage for various earthquake vibration levels expressed as peak ground velocities (PGV). In this fashion, the probability of light damage state (DS1) for fired-clay-brick masonry is characterised.
To accurately determine the probability of damage, several fundamental aspects about the masonry need to be known; in this paper, these are obtained from earlier studies. Firstly, the characterisation of the material properties and their variability is procured from the extensive research of Jafari et al. 2017; here, the expected variations of, for instance, the tensile strength or the Young's Modulus of the masonry, are investigated from replicated material tests based on samples from real structures in the Groningen region. Secondly, to determine the probability of damage, numerous and accurate models of structures need to be evaluated; these models call for a representative view of reality. Experiments conducted on masonry walls, monitoring the initiation and propagation of cracks, were used in various studies to develop calibrated computational models tuned to replicate the behaviour of the walls not only in terms of stiffness, strength, and hysteresis but also of the intensity and pattern of cracks (Korswagen et al. 2019a(Korswagen et al. , b, 2020a. These models are expanded herein to simulate the behaviour of masonry walls within the constraints of a structure interacting with the soil and the vibrations of earthquakes; the latter component having been also evaluated in a parallel study (Longo et al. 2020). The deformation of the soil, during earthquake vibrations and in prior potentially damaging processes, was seen to lead to increased light damage and earlier damage detection (Korswagen et al. 2019b(Korswagen et al. , 2020cStaalduinen et al. 2018). Thirdly, to effectively quantify damage, a continuous scalar parameter has been proposed to express the intensity of the varying crack patterns; the parameter, Ψ (Psi), tailored to evaluate bare masonry walls, is based on the ease of detectability and repair of the cracks and comprises the number of cracks, crack width and crack length (Korswagen et al. 2019a;. Moreover, it helps in the validation between experiments and models and with the assessment of (the progression of) damage.
These prior studies have dealt with the deterministic, yet vital background for the probabilistic assessment performed in this paper. This data and tools are employed to further the insight into light damage. Using physical modelling to characterise probabilistic behaviour can be a powerful tool (Fiore et al. 2014;D. Rapone et al. 2020). Figure 1 summarises the background studies and depicts the remaining three steps treated in this paper. The extension or extrapolation with the existing calibrated model is explored in section two and used to create a large pool of results where the relationships between damage and material strength, soil type, earthquake intensity, etc. can be observed. Such a regression model, investigating the physical implications of the established relationships, is developed in section three and used in section four to, on the basis of a Monte Carlo simulation, compute the pursued fragility curves. Finally, the last section compares the fragility curves with other attempts in literature to assess light damage and summarises the main conclusions found herein.

Background on validation experiments
The quasi-static, cyclic experiments performed in related studies (Korswagen et al. 2019a(Korswagen et al. , b, 2020a are briefly summarised in this section: Five nominally identical, single-wythe walls of fired-clay brick masonry approximately 3.1 m wide and 2.7 m tall with an asymmetrically-placed window opening were tested quasi-statically in-plane by enforcing a controlled displacement at the top of the walls. The walls were pre-loaded vertically to produce an average stress of 0.12 MPa resembling moderate vertical loads; the top of the wall was allowed to rotate freely mimicking a 'cantilever' configuration. The in-plane displacement was applied cyclically with repeated steps of increasing amplitude; the surface of the walls was surveyed with Digital Image Correlation (Kumar et al. 2019;García-Macías and Ubertini 2019) to detect and monitor the formation and progression of cracks at a resolution of 20 µm. Two additional walls of identical material without a window opening ) were tested at a higher vertical overburden of 0.46 MPa and a 'double-clamped' configuration, meaning that the bottom and top edges of the walls were kept parallel during the tests. The first set of walls were seen to correspond to the flexural behaviour, while the latter represented the shear behaviour of masonry walls. Furthermore, two walls with a window opening and pre-damaged crack interfaces (Korswagen et al. 2019b) were tested in a manner identical to that of the five walls with opening for a window. The pre-damaged interfaces were achieved by incorporating a PVC sheet on some mortar joints strategically located at the corners of the window opening. The built-in pre-damaged pattern was designed to mimic cracking caused by a differential settlement profile and was determined with numerical models and assessment of real case studies (Staalduinen et al. 2018).
The tests imposed in-plane drifts to the walls of up to 0.1% with some cracks reaching a width of up to 2 mm. The cracks were observed to propagate when identical drifts were enforced repeatedly, though greater crack growth occurred at cycles of amplitude larger than that of the preceding cycle. Pre-existing cracks modified the initial stiffness of the wall and the ultimate crack pattern slightly (Korswagen et al. 2019b). Cracks in the flexural walls were mostly horizontal and stair-case-diagonal and propagated off from the corners of the window opening, while cracks in the windowless shear walls were mostly diagonal but steeper than the stair-case motif and started off from the centre of the walls. Additionally, ten spandrel-type wallets 1.5 m in width and 0.5 m in height were tested in a four-point bending setup such that a flexural vertical crack would cross the specimens in the centre. The cracks were observed to zig-zag around the mortar joints of the brick masonry resulting in a toothed crack pattern.

Background on calibration models (experiment-based)
The various experiments described, together with standard characterisation tests for masonry such as bond-wrench tests, compression tests on wallets, in-plane bending tests, and flexural and compressive tests on units and mortar specimens, were used to develop a modelling strategy capable of satisfactorily replicating the behaviour of the specimens in terms of stiffness, strength, hysteresis and crack patterns. These models (Korswagen et al. 2019a(Korswagen et al. , b, 2020a are succinctly discussed here. Modelling was done following the macro approach for masonry (Lourenço et al. 1995) and using the finite element method (FEM) with the DIANA-FEA program and two-dimensional, plane-stress, quadrilateral, 8-node, quadratic elements (CQ16M) with a mesh size of 50 mm considering a full 3 × 3 Gauss integration scheme. A non-linear material model was used to properly simulate the (cyclic) cracking behaviour of masonry; the orthotropic total-strain-based masonry model (also called Engineering Masonry Model-EMM in DIANA) was selected as constitutive model (Schreppers et al. 2016;Rots et al. 1985Rots et al. , 2016Lourenço et al. 1995), with different values of inelastic and elastic properties for the two main directions: perpendicular to and aligned with the bed-joints. A single set of parameters for the material model, collected in Table 1, was used for all calibration models since the material of the masonry in the experiments was identical within the expected variability; this is an important aspect of the calibration. The top and bottom steel beams were included in the model as compatible elastic beam elements and the masonry was tied to them; the aim was to accurately include the experimental boundaries. The models considered physical and geometrical non-linearities and were solved with a parallel direct sparse method, with a tolerance of 10 -8 , in a quasi-Newton (secant) incremental-iterative approach with either a force norm (tolerance 10 -2 ) or an energy norm (tol. 10 -4 ). An overview of the four experiment types used to calibrate the modelling strategy is presented in Fig. 2 Fig. 2 Graphic overview of experimental tests (left column), comparable crack patterns in finite element models (middle column) and comparison between experimental and computational force-displacement curves (right column). a Vertical displacement field captured with DIC from a laboratory wall with a window. b FEM model of the same wall. c Comparison of the force-displacement curve between the experimental (a) and the model wall (b). d Automatic crack identification from the DIC data for a laboratory wall without a window. e Crack strain plot of a FE model of the same windowless wall. f force displacement curve of (d) and (e). g Crack locations for a wall with pre-damage. h FE model of such a wall, for tensile (left) and shear cracking (right). i Force-displacement graph for (g) and (h). j Spandrel-type wallet (window bank) surveyed with DIC displaying crack strains. k FEM micro-model of said wallet. l Forcedisplacement graph of the wallet for several experiments with monotonic loading (j) and two FEM models differing in the Young's modulus of the masonry material (k) displacement fields and crack patterns and the middle column shows the equivalent view from a calibration FE model; the right column compares the force-displacement curve between experiment and model. Similarly, the first row depicts one of the tests on a virgin wall with an opening, while the second row outlines the case of a windowless wall. In both cases the loading protocol, inferred from the graphs, is asymmetric, consisting of a one-way cyclic displacement and a two-way cyclic portion with larger positive amplitude, respectively; and accordingly, the crack patterns displayed are not symmetric. The resolution of the imaging and processing systems for the laboratory cracks produced clear cracks, and the smeared models presented sufficient crack-localisation that the experimental patterns could be recognised; the crack locations, and their properties such as crack width, length and overall damage measured by Ψ (Korswagen et al. 2019a;, contribute to the calibration process. The following two rows showcase the experiments with walls with two pre-cracked mortar joints and the wallets or spandrels; for the former, the full two-way cyclic loading protocol is presented, while for the latter, some of the monotonic tests are compared to two model possibilities. The figure presents, at a glance, the various calibration campaigns undergone.

Setup of extrapolation models (based on field cases)
The above models, that were calibrated and validated against the laboratory tests, are subsequently extended into the nonlinear time-history domain so as to assess the effect of seismic vibrations and replace the laboratory boundaries with the constraints provided to a masonry wall by a real structure; in this light, these models become extrapolations from the validated laboratory scenario to real field cases while maintaining the geometry of the wall. The adaptations applied to the models are presented in Fig. 3   • Overburden a vertical pre-compression of 0.12 MPa is applied at the top edge after the self-weight gravity load. This stress corresponds to the weight of one concrete floor or two timber storeys. • Top boundary no steel beam is used at the top of the wall as in the laboratory; instead, a stiffness-free line mass element is incorporated to mimic the dynamic effect of a floor. The mass is set at 10 ton to simulate the presence of a fictitious rigid floor, an upper masonry wall, an upper floor and a roof. While the wall doesn't carry the vertical weight of these components, their mass does affect the horizontal in-plane behaviour of the wall; as such, the line mass does not influence the static condition. The mass is also used to calibrate the natural periods of the walls to match those of low-rise masonry buildings from which the walls are extracted. The flexible wall has a fundamental natural frequency of about 9 Hz and the stiff wall of about 13 Hz; depending on material, soil, and initial masonry condition, these periods can vary by up to ± 2 Hz. • Lateral boundaries two vertical, linear-elastic beam elements, with a fictitious crosssection of one brick, are added to the lateral edges to simulate the constraints of transversal walls. These elements have an elastic modulus equal to 1786 MPa and a Poisson's ratio of 0.16. The density is set to account for the out-of-plane effect of a wall of 3 m in length along the transversal direction. • Rayleigh damping a 2% damping based on the first two modes, as judged by order of participating mass, is included in the model. • Soil-structure interaction the bottom steel beam used during the calibration phase is replaced with a combination of springs and dashpots. These are placed underneath the wall in order to simulate the effect of the soil on the structure. Stiffness and damping are provided by two separate interface elements: (a) A non linear 2D line boundary interface (CL6TM-NL), which accounts for the stiffness and for the soil's behaviour, is characterised with a Coulomb friction law. The cohesion is set to 0.1 MPa while the friction and dilatancy angles to 0.61 rad. Additionally, a no-tension criterion is enforced for the springs representing the fact that no tension can develop between soil and foundation. (b) A similar, yet linear-elastic interface (CL6TM-LE), which includes the damping coefficients and offers no stiffness, is also included. The values of stiffness and damping are varied for two types of soil profiles (and two types of façades) as presented in Table 2. A stiff soil, with higher stiffness and better damping, is contrasted to a poor or soft soil. The methodology employed to determine these coefficients is adapted from NEHRP (2012) and described extensively in Longo et al. (2020), while the representative soil profiles were obtained from Deltares (2015) using insight discussed by Kruiver et al. (2017).
• Input motion the load introduced as an in-plane, quasi-static displacement at the top of the wall in the experiment is replaced by horizontal and vertical accelerations applied at the base of the wall. This extrapolation is twofold: first, the accelerations become inertia loads, focused at the top but also distributed throughout the wall, this modifies the validated behaviour slightly; and second, the appearance of cracks could be different when tested dynamically, however, the damage behaviour with dynamic tests has not yet been studied; hence, this becomes a reservation in the present study. Nevertheless, the assumption of transitory or transient damage as the residual damage is conservative .
The accelerations are exerted at the fixed end of the non-linear interface discussed in the previous point, and correspond to four different records: two belong to the earthquake occurred in Zeerijp on the 8 of January 2018 (M = 3.4) and measured at the stations of Garsthuizen and Appingedam, located 2.5 km and 7.5 km from the epicentre, respectively (Bezemer and Elk 2018; KNMI 2020), and two to the earthquake of Westerwijtwerd of May 22nd, 2019 (M = 3.4) recorded at the Stedum (3.2 km) and Hoeksmeer (9.7 km) stations. These natural motions are depicted in Fig. 4 and are categorised as 'near' and 'far' with the amplitude of the motions normalised to 1. This grouping, based on the distance of the station to the epicentre, also corresponds to similar properties between records even if they belong to different earthquakes and are measured at different stations. In fact, Fig. 4 shows that the horizontal component of the seismic motion arrives noticeably earlier at the near locations and contains a fewer number of effective cycles, while the far records display a larger number of cycles of similar amplitude. This is also reflected with a Fast-Fourier Transform analysis which reveals that the horizontal components of the far records are heavily loaded with low vibration frequencies in the order of 2 to 4 Hz, while for the vertical components the frequency content is more alike between the records. The Arias plots also show how the energy of the near records is focused in a small part of the record, while the far records spread their energy over a longer significant duration. In comparison to their horizontal components, the far vertical records exhibit a PGV amplitude of about one third of its horizontal component, while the near vertical records present a PGV amplitude of about one fifth of its horizontal component. Finally, in terms of spectral amplification, the horizontal components of the records reach similar amplification but the far records are active over a much longer interval of structural periods. In comparison, the vertical components, especially of the far records, display a higher dynamic amplification but are limited to periods close to 0.1 seconds. These variations are present when comparing both earthquake types as long as both are scaled to the same value of peak ground velocity. Note that the variations in the records, and in particular the spectral response, do not only originate from their distinct epicentral distances, but are also the result of the different site conditions, such as the local soil characteristics, the composition of the soil between epicentre and station, the orientation of the fault causing the earthquake in respect to the station, and many other subtle factors which do not play a role in this study (Green et al. 2020). Nonetheless, the four records selected, recorded at different stations and distinct epicentral distances, are used to provide the 'record to record' variability in the models. While structural analyses are usually conducted using a larger pool of records (of seven  Table 5. Records are normalised to the horizontal peak ground velocity for comparison purposes or more), the demanding computational expense of the NLTH-Analyses and their variations (see further on) necessitated this compromise.
In this manner, the light-damage impact of earthquake vibrations can be evaluated on a single masonry wall; the complexity of a three-dimensional structure and the effect of the soil are included in a model that can be reliably assessed for a multitude of variations while observing the initiation and propagation of damage in the form of cracks. In this context, the variability of light damage is characterised through seven key parameters listed in Table 3. The parameters are categorised into 'binary parameters', those that only have two or four distinct possibilities and for which no trend can thus be observed, and into 'continuous parameters', which may assume any value within a given range. All these parameters are discussed next: • Firstly, two wall geometries are modelled, with and without a window opening, to consider a more flexible wall with predominantly flexural behaviour, and a more rigid wall with prevalent shear behaviour; these two typologies correspond to the two types of wall geometry tested experimentally and represent the first binary parameter. • Secondly, the soil can be varied; the two soil variations used in this study have been discussed in Table 2. • Thirdly, the record-to-record variability is provided by four earthquake records; the type of earthquake motion is characterised by the distance to the epicentre with 'near' and 'far' scenarios and belong to two earthquake events; the difference in frequency content, number of effective cycles, and relationship between horizontal and vertical motions for example, are inherent to these records and serve to depict the variability in earthquake motions. • Fourthly, the variability in the intensity of the motion is considered by scaling the PGV of the horizontal component of the records (and proportionally, the vertical component) to eight different amplitudes: 2, 4, 8, 16, 32, 64, 96, and 128 mm/s. It's important to note that in this context, the distance to the epicentre is decoupled from the amplitude of the record. Moreover, this is the first continuous parameter since the PGV can adopt any value within this range. • Fifthly, the quantity of earthquake events is also varied; applying the same record for a second time leads to greater accumulated damage in the non-linear model. Thereby, the effect of repeated events on masonry structures can be quantified. • Sixthly, the masonry material must be varied; three scenarios are considered for the fired-clay brick masonry corresponding to a weak, standard, and strong masonry. Their differences are gathered in Table 4 and are based on an analysis of masonry characterisation tests (Jafari et al. 2017). • Finally, the initial condition of the masonry wall is also varied; this requires some discussion. The intensity of the (pre-)damage is measured using the Ψ damage parameter Material strength (relative) 0.5 to 1.5 Zeerijp (Z) or Westerwijtwerd (W) Initial damage (Ψ 0 ) 0 to 1.5 (Korswagen et al. 2019a;; this parameter summarises the number, width, and length of cracks into a single scalar such that the progression of damage can be quantitatively observed. Moreover, the value of Ψ = 1 corresponds to the threshold of visible damage when cracks have exceeded a width of 0.1 mm; see also Tables 5 and 6. A virgin (Ψ 0 = 0) and three pre-damaged initial conditions (Ψ 0 > 0) are considered; the three scenarios where the wall is already damaged when the earthquake motion is applied, correspond to distinct intensities of differential settlement damage. A settlement profile is enforced as a displacement of the fixed side of the non-linear boundary interface; the normalised profile is shown in Fig. 5 and was fabricated based on the study of real-case deformations via the measurement of masonry bed joint levelling (Staalduinen et al. 2018). The amplitude of the profile is adapted such that three intensities of pre-damage are achieved for the wall: low (0.25 < Ψ 0 < 0.75), moderate (0.75 < Ψ 0 < 1.25), and high pre-damage (1.25 < Ψ 0 < 1.75). It should be noted that  sometimes, models of different typology with the same pre-damage value may display a different failure mechanism and, that these initial values are used in the context of DS1 and so the terms 'moderate' or 'high' are relative to these values as pre-damage to assess light damage.
Each analysis is divided into three phases: self-weight and gravity loads, settlement (if present), and earthquake loading. All the models are investigated with the same analysis procedure, keeping load steps, iterative method, norms, and tolerances identical. For the settlement phase the load step changes based on the different properties and the intensity of pre-damage required. After the self-weight and settlement case, the non-linear analysis is continued in a time-history setting with the dynamic load being applied to the deformed model with time-steps of 2 ms. As with the calibration models, physical and geometrical nonlinear effects are considered; yet, in addition, transient effects are activated while considering the dynamic effects using a consistent mass and damping matrix. The Newmark-Beta method (β = 0.25, γ = 0.50) is used as time integration method for these non-linear time-history analyses. Moreover, the Rayleigh damping is kept using the then current stiffness matrix.

Results from extrapolation models
The permutation of the seven parameters to vary in the models results in more than six thousand different scenarios; however, only 3840 cases were run by neglecting cases were key, previously-run models indicated that no damage was expected. For example, for the cases where the earthquake intensity of 4 and 8 mm/s did not lead to any damage (increase) in the model, the case of 2 mm/s was not run. The result of each model was evaluated by automatically computing its damage value Ψ and contrasted to the initial damage Ψ 0 . Figure 6 presents a coloured overview of the results; note the eight zones delimited by dark lines which correspond to the 2 × 2 × 4 permutations of the three binary variables (wall type, soil, and earthquake type and event) and that each square is divided into two triangles for the two events of Zeerijp and Westerwijtwerd. The effect of each variable is discussed next while highlighting some examples from the large pool of analyses.
First, the wall type or façade type revealed that the flexible wall (A) showed a higher increment of damage (∆Ψ) in respect to façade B, consistently, at least for the neartype earthquake records. An example is provided in Fig. 7a in which two damage patterns of different façades are presented. The window opening of façade A contributes to an easier initiation of cracks due to stress concentration at the window corners. The  Fig. 7a, when the two façades are placed on soil A, having weak material, no pre-damage, and subjected to a single, near Zeerijp event of 32 mm/s, façade A experiences 25% more damage than façade B. In other words, for most of the cases, the wall with the window appears more vulnerable to damage than the rigid wall. It should be noted that these results display the maximum of crack strain in all Gauss points that occurred over the entire NLTH analysis. This explains the rather wide zones or 'clouds' of smeared cracks. At certain discrete times, the crack strain plots show a more localised pattern as some Gauss points soften while others unload. As an alternative, discrete crack/slip micro-models in predefined interface elements were used to study these localisations ); yet, these more detailed models are severely more computationally expensive and did not lead to significantly better calibration results; hence, they were not considered for this large number of analyses. Secondly, the soil type: the softer soil (B) requires a higher settlement amplitude to produce the same damage level obtained in a façade modelled with soil A. This is because, for the same enforced soil deformation, the more flexible soil better accommodates the desired deformation of the structure, while the stiffer soil leads to sharper curvatures and hence higher stresses and damage in the walls. Figure 7b shows the virgin, rigid wall with weak material and subjected to a single, near earthquake (ZN) of 96 mm/s as PGV, analysed with two different soils. A similar crack pattern is detected for the two models, which depict diagonal shear cracking at the middle of the wall. Cracks in models with soil A appear more concentrated and wider than the smeared shear cracking of the one on the weaker soil B. In this case the difference in damage at the end of the protocol reaches about 20%.
For the earthquake type, near or far from the epicentre using the two records of the Zeerijp event, Fig. 7c presents the effect on façade B on soil B with standard material and moderate pre-damage. Both earthquake types are scaled to a PGV of 32 mm/s and a single motion is run. The results show a different crack propagation; while for the near-field motion the existing crack pattern grows mainly in length with some additional parallel lateral cracks, the far motion depicts mainly an increase in width without affecting much its length. This behaviour is frequently seen in façade B when the far record often produces higher damage. On the other hand, the flexible façade (A) appears more sensitive to the near motion. This observation is likely related to the different frequencies present in the two signals which interact in distinct manners with the two types of walls.
Next, Fig. 8a. depicts the effect of different pre-damage situations on the final damage condition at the end of one motion; only the cases subjected to settlement are shown in the   graph. The analyses concern façade A with standard material on soil A and subjected to one near earthquake (ZN) of 16 mm/s. The pre-damage values (Ψ 0 ) are about 0.5, 1.0 and 1.5 respectively for low, moderate and high pre-damage cases. As can be seen, increasing the damage at the initial condition, produces higher final damage. The detected values are 1.2, 1.4 and 1.6 Ψ with increasing level of pre-damage. The model of the same category without initial settlement loading, generates a final damage of 0.8 Ψ. It must be noted that, the increment of damage ΔΨ diminishes when a higher pre-damage is adopted. In fact, considering the initial state, the difference between final and initial damage (ΔΨ) is equal to 0.8, 0.7, 0.4, and 0.1 Ψ, respectively for no, low, moderate and high pre-damage. This observation can be generalised for many other cases.
Regarding the quantity or number of earthquake records, the cumulation of damage gained from multiple earthquakes in a row is summarised in Fig. 8b. The presented façade B on soil A with weak material and no pre-damage is subjected to three different number (repetitions) of near-type earthquakes (ZN), all with PGV of 16 mm/s; one, four and eight motions are applied. As expected, a higher number of earthquakes leads to a higher value of final damage. After each earthquake motion, the cracks increase in width and in length. The increase is minor, about 2 to 5%, but consistent throughout the entire table of results; lower PGVs are also associated with a smaller increment rate.  Similarly, the variation in material presents a congruous observation: the weaker masonry material leads to a more flexible model, but also to lower values of resistance to tensile and shear failure, consequently, a larger increment of damage can be observed when the material properties are weak. An example of this trend is shown in Fig. 8c where the flexible façade A is subjected to a 96 mm/s far-field earthquake (ZF). The model is placed on soil A and has no pre-damage, and the results at the end of the motion of the three different material configurations are shown in the figure. The façade with weak material develops two cracks at the top window corners, together with damage in the top spandrel. The two stronger material variations depict additional tiny cracks located at the other window corners and a more distributed damage to spandrel and lateral sides; nevertheless, the width of these cracks and their length result in a final damage value lower than that of the case with the weak material.
Finally, varying the intensity of the earthquake motion, where the intensity is given by the scaled value of peak ground velocity, offers the expected relationship that, at higher intensity, the increase in damage is also higher. The effect of different PGVs on façade A is shown as an example in Fig. 8d; the model is placed on soil B, with standard material and virgin initial condition. The wall is subjected to one earthquake of the near type (ZN) and three PGVs are selected: 8, 16 and 64 mm/s. The final damage is remarkably different between the three models with an almost 100% increase between them. This underlines the large effect of this parameter on the results.

Regression model
The trends observed in the multitude of analyses run to produce Fig. 6 can be captured in a regression model that incorporates the effect of every parameter varied and the relationship between the most relevant of these parameters. To avoid a 'black box' effect where an arbitrary model is fit precisely to the existing data points but may produce unexpected results outside of the domain of the data, a 'reasoned model' is employed where the relationship between the input parameters and damage can be understood. First, clear trends are gathered in the following order: • The higher the PGV, the more damage is expected; • The higher the pre-damage value, the higher the final damage. However, the difference between initial damage and final damage diminishes when the initial damage is higher; • The weaker the material, the higher the final damage; • The more flexible soil (B) led to lower damage, albeit not significantly so. Note that the amplification that soft soils produce is not considered here since the PGV is being used directly; • The more flexible façade (A) showed higher damage for the near records, while the rigid façade showed more damage for the far records; • The larger the number of earthquakes, the higher the final damage. This increase was marginal but consistent.
To gain additional insight into these relationships, the three binary variables (soil, façade, and earthquake record) were considered independently while the four continuous variables (PGV, material, number of events and pre-damage) were included in a multivariate function, and a linear Pearson correlation between the increase in damage (ΔΨ) and the various parameters was observed; see Table 7. The Pearson correlation values range from − 1 to + 1, indicating perfect inverse correlation or perfect positive correlation, respectively. A perfect correlation is achieved when changes in one parameter are precisely reflected in the other. First, it can be observed that the correlation of the material is slightly negative; this means that as the material strength increases, the damage decreases. Similarly, the initial damage has also a negative correlation, meaning that, when the initial damage is high, the increase in damage reduces; this influence is more important than that of the material since the values are closer to − 1. Next, the number of earthquakes has a small positive correlation indicating that a larger number of earthquakes will also cause more damage. Comparably, the PGV has a strong positive correlation; this is expected as the PGV is the hazard variable causing damage. Lastly, the final value of Psi has been included in the table to illustrate that the increase in damage is not perfectly correlated with the final damage; this hints at the importance of the initial damage condition.
Moreover, differences between the soil, façade and earthquake type can be also recognised: First, the initial damage and material strength seem to be more relevant for earthquakes occurring nearby, while earthquakes farther away are more damaging regardless. This is also reflected in the higher correlation with the PGV. Conversely, the soft soil (B) increases the importance of the material and the initial damage, suggesting that the damage of structures on soft soil is more dependent on their condition, where also the repetition of earthquakes has a larger effect as seen from the higher values associated with the correlation of the number of earthquakes. Interestingly, the far records show a slight inverse of the correlation between the initial damage and the increase in damage for the stiff facade, suggesting that in these cases, the more damaged and thus more flexible state of the wall in combination with the specific spectral content of the far records leads to an increase in damage. Yet, a distinct behaviour difference cannot be established between the stiffer façade without a window opening (B), and the more flexible façade (A). Accordingly, the model function was shaped so as to fulfil the observations presented and the trends observed in Table 7 while simultaneously employing the fewest possible The blueish and greenish tints highlight negative and positive values respectively number of regression coefficients. Various regression models were tested and evaluated based on the mean value of the residuals, the number of coefficients, and their fit to the observations. The best model is herein presented; it constitutes a scaled sum of a logistic component and a linear component with the PGV as the argument and the increase in damage (ΔΨ) as the ordinate value. Both components of the function depend on the inverse of the material and the initial damage, while the number of earthquakes and initial damage also participate in the scaling factor.
where β 1 is the scaling factor in the function, β 2 is the exponent of the logistic component and β 3 is the linear component; these are in turn defined as: where N is the number of consecutive, identical events, m is the normalised material strength given by the ratio of tensile strength over the tensile strength of the standard material, and α are the eleven regression coefficients. Their values are gathered in Table 8 for each of the sixteen combinations of binary variables. The fit of the model is evaluated using the mean of the absolute values of the residuals, which is in average 0.15 Ψ points with a deviation of 0.16 points over the sixteen sets of coefficients. This is deemed a reasonable fit due to the scatter in the analyses' results; the logistic component of the model allows for a good fit to low values of damage for low values of PGV and introduces a steepening of the relationship at a certain PGV value (around 10 mm/s) which follows the data points adequately; this can be observed in Fig. 9. This figure also illustrates how damage remains close to zero for very low values of PGV but then starts to increase rapidly between values of 5 to 20 mm/s; however, for stronger materials (m = 1.5), this transition occurs more gradually into well around 50 mm/s. Therefore, the value of Ψ ≥ 1 corresponding to the start of visible damage, appears significantly earlier for the weaker masonry. The regression coefficients are restricted to an interval that leads to a smooth surface over the material parameter where only three data points are available (weak, standard, and strong). This follows the assumption that the behaviour of real world structures is expected to be a smooth curve. The values of the regression coefficients, like the Pearson correlation values, hint at the relationship between the parameters and the light damage of a masonry wall. The regression coefficient α 2 for example, indicates the influence of the number of earthquakes; a positive number slightly larger than zero suggests a proportional relationship between the number of earthquakes and the increase in damage. Similarly, α 3 (1) 1 3 establishes a link between the initial condition and the increase in damage; its values in Table 8 suggest that this link is stronger for the flexible façade but weaker for the stiffer façade. Next, α4 and α5 determine the effect of the initial damage on the logistic part of the function (principally for lower PGV values), while α6 and α7, indicate its effect on the linear component of the function. Likewise, α8 and α9 link the material to the logistic part of the model and α10 and α11 to the linear part.  Table 8) and its material strength where the regression coefficients (α8 to α11) show comparatively high values.
To improve the behaviour of the regression for the sensitive and more commonlyoccurring low values of PGV, the data points below or equal to 16 mm/s are considered twice as important as the higher values by using a weight of 2 in the least-squares model fit. The threshold of 16 mm/s was chosen, first, because extrapolation models were available for all four earthquake repetitions while the following model set, run at 32 mm/s, had twice the PGV value and was present in only three of the sets. This means that the confidence to the model fit up to 16 mm/s is higher than for points later on, and setting the threshold at any value between 16 and 31 would yield the same result. Secondly, in comparison to the PGV values that have been recorded in the field and to the four earthquake motions, which reached PGV values below 32 mm/s, a threshold above this interval would be too high. Since the present study focuses on damage at low PGV values, these should be preferred in the fit. Moreover, as seen in Table 5, the original records achieved PGV values between 1 and 32 mm/s; scaling these records too far beyond their original PGVs may introduce a bias further decreasing the confidence in the data points above 32 mm/s. On the other hand, the value of 2 for the weight of the fit is rather arbitrary; for a higher weight, the points above 16 mm/s become irrelevant and are ignored, so the largest reasonable value was selected.

Model uncertainty
The imperfect fit of the regression model, due partly to the non-smooth trends of the data points provided by the extrapolation FE models, can be captured by a model uncertainty parameter, ε. The parameter is considered by modifying Eq. 1 such that the increase in damage, ΔΨ, is affected by ε, which, using a Generalised Extreme Value, is randomly distributed, satisfying the condition that the final damage can never be lower than the initial damage; this is expressed in the following three equations.
The value and distribution of ε can be approximated to the average mean of the absolute value of the residuals of the model regression and distributed normally. However, the variability in the fit appears to be dependent on the PGV as is visualised in Fig. 10 and, moreover, the residuals are better represented with an extreme value distribution where the actual value of Ψ may be much larger than what is predicted by the regression model. A threeparameter Weibull distribution, or Generalised Extreme Value distribution, accurately follows the shape of the residuals and can be adjusted such that its parameters (shape, scale and location, k gev , σ gev and µ gev , see Table 9) become functions of PGV as is specified in Eq. 9 and are thus defined for any value of PGV (PGV ≥ 1 mm/s). The use of a Weibull distribution is slightly conservative because it allows for relatively large values of ε due to its long tail to the right, yet limits the underestimation of Ψ because of the truncated tail at its left. In this sense, the uncertainty parameter can be assumed to include not only a degree of aleatoric uncertainty but also the epistemic uncertainty in the FE modelling approach and crack characterisation processes, where unforeseen conditions may trigger much-largerthan-anticipated damage.

Probabilistic distributions and parameters
The regression model can be used within a Monte-Carlo simulation to determine the probability of damage. This type of simulation takes one sample for each of the partaking variables and determines the outcome with this combination of values; sampling is repeated until sufficient outcomes have been evaluated to establish a reliable result. Therefore, it is first necessary to assign a distribution from which to sample each of the relevant parameters. First, the material strength is present in the regression model as the ratio between the direct tensile strength and the tensile strength of the standard or mean material; this conveys a set of material properties from which the tensile strength was deemed most influential in the FE models. The coefficient of variation measured for the tensile strength in material tests is about 30% and was normally distributed; hence, the variable m can be characterised with a normal distribution with a mean of 1 and a deviation of 0.3. To avoid sampling material values too far off from the calibrated values (run for 0.7, 1.0 and 1.3), the distribution is truncated at two times the standard deviation.
Next, the initial condition of the masonry plays a large role in the final damage outcome; however, the condition of the structures in the region is still a topic of ongoing study and, the large variability in damaging actions like settlements and (autogenous) stresses from restrained shrinkage or thermal movements, make it difficult and also unreasonable to consider any distribution for this variable. Instead, greater insight can be obtained by Similarly, the PGV can be used as the hazard variable, that, against which the probability of light damage can be judged. Consequently, the simulations are run for various values of PGV in the interval from 1 to 100 mm/s to characterise the relation between PGV and damage. Moreover, the number of identical, repeated events is excluded from this simulation so as to consider the vulnerability of the masonry to one seismic event. Multiple seismic events are unlikely to be of identical amplitude so this aspect requires further study.
Next, two instances of a masonry wall have been modelled, a rigid, windowless shear wall and a flexible wall with an opening for a window. This was included in the model to provide some variability to the geometry and the dynamic behaviour of the structure; it is not feasible to attribute one of these two states to existing structures and then determine whether most structures correspond to the flexible or the rigid type; rather, it is expected that most buildings will have some masonry walls that behave more alike the former and some, more alike the latter. And, by including this parameter, the model is expected to be more applicable to a greater number of buildings. That being so, the parameter is characterised with a uniform distribution where each wall type is arbitrarily assigned a 50% probability.
Furthermore, characterising the type of earthquake motion, near or far, or the type of soil, stiff or soft, requires some additional discussion. First, intuitively, there is a correlation  Table 9 Coefficients to determine k, σ, and µ of the generalised extreme value distribution for the PGVdependent uncertainty of the regression model. See Eq. 9 between the PGV and whether an earthquake is near or far, with high values of PGV more likely to be caused by an event nearby. Secondly, the distribution of the epicentres in relation to the structures will also play a role in dictating whether earthquakes are distant or near. Therefore, a brief statistical study on historical data was conducted to characterise the type of earthquake and establish its correlation to the PGV. The events in the period between January of 1993 and December of 2019 with a magnitude of at least Mw ≥ 1.5 were observed and the PGV was computed throughout the region using ground motion prediction equations (GMPEs) as developed by Bommer et al. (2017); the PGV was sampled probabilistically. A simplified database with the number and location of masonry structures Arup (2018) was used to determine the distance and intensity of the seismic events in relation to the buildings in the region. A ratio between the near and far events, as categorised by a threshold varied between 2 and 5 km, was established for various PGV intervals; this is depicted in Fig. 11. This graphs shows how events with a local PGV greater than 40 mm/s for example, corresponded most of the times to events within a radius of 2 km; while events with a PGV smaller than 5 mm/s were much more likely to come from an event with a distant epicentre. As the threshold to distinguish between near and far events increases from 2 to 5 km, the likelihood of events being nearby reasonably increases. For a random event, this simulation employs a correlation function which establishes a 10% probability that the event is from a nearby source if the PGV is 1 mm/s and a 90% probability if the PGV is 35 mm/s. This function has been arbitrarily defined for any value of PGV based on the various thresholds plotted in Fig. 11. Finally, from the same analyses, the likelihood of stiff or soft soil was statistically determined for buildings exposed to at least 1 mm/s. The vertical shear wave velocity for the first 10 m of soil (V s10 ) was computed throughout the region based on a micro-zonation study ; locations with V s10 greater than 100 m/s were deemed to be stiff soil, while buildings on locations with a lower V s10 were considered to be on soft soil. This threshold value also matches the soil parameters presented in Table 2. The resulting binary distribution assigns a 95% to soil A (stiff) and 5% to soil B (soft).

Fig. 11
Correlation between PGV and near-type records from statistical, historical data. Low PGV values sometimes correspond to events with an epicentre located farther away, while PGV values above 10 mm/s always correspond to events originating less than 3 km away

Results of MonteCarlo simulation
The simulation samples a value of material strength, soil type, earthquake type, model uncertainty, and wall type from the aforementioned distributions; and, using Eq. 6, it determines the final damage for discrete values of initial condition (Ψ 0 = 0, 0.5, 1, and 1.5) and PGV. The final values are binned into categories for Ψ between zero and three in half point intervals, as presented in Table 6. The simulation was run to exceed sufficiency with at least ten million points for every discrete combination of initial condition and PGV. The results are illustrated in Fig. 12 as fragility curves. Here, the probability of exceeding various light damage levels is depicted based on the initial condition of the masonry structure and the value of PGV; the discrete values obtained from the simulation are marked and a simple spline is fitted. At 10 mm/s, the curves show only a 5% probability of reaching or exceeding visible damage (Ψ ≥ 1) when no initial damage is present (Fig. 12a), but this increases drastically to about 20% when some initial damage, yet invisible to the naked eye (Ψ 0 = 0.5, Fig. 12b), is present. The curves are relatively flat, indicating a large spread in the variability of the input parameters; this is also alike to field-cases where the multitude of structural configurations generate wide probabilistic ranges. However, for low pre-damage (Ψ 0 ≤ 0.5), the curves for visible damage are initially steeper suggesting that masonry is vulnerable to displaying light cracks for low values of PGV; this observation is not unexpected but the curves determined herein define it quantitatively. The graphs in Fig. 12c and d, which correspond to already visible initial damage (Ψ 0 ≥ 1), show that light damage in unlikely to aggravate severely unless the intensity of the earthquake vibrations becomes relatively large with PGV values greater than 40 mm/s. This value corresponds roughly to a PGA of 0.12 g, which in many codes of seismically-active regions is labeled as 'light' or 'very light' vibration intensity (FEMA 2013). These observations and some other key values are extracted from the curves and collected in Tables 10 and 11.
The simulations also revealed that the records of the near type were more damaging to the masonry walls (see also Fig. 6) and thus increase the probability of damage. In general, low-rise masonry buildings are more conservatively assessed using 'near-type' records with fewer cycles and a short significant duration. Moreover, this study employs two types of masonry walls to represent the in-plane behaviour of building components for masonry structures in general; it is clear that other types of façade, more complex, with more openings or different width-to-height ratios, will display a diverging picture of damage. Nevertheless, the variability provided by the two façades included should capture the probabilistic spread of the population and is reflected in the curves. Yet, it must be emphasised that this study focuses on the damage that appears in the bare masonry or at continuous corners of two masonry walls and not at the connections between walls and floors or between two adjacent walls. In this context, it is expected that joints between adjacent walls are more likely to develop cracks just as dilation joints are designed to gather strains and prevent cracks in other parts of the structure. This type of light damage is not contemplated by the curves and would also be more strenuous to quantify. Finally, the use of a continuous parameter to describe damage coupled to a description of the type of damage will result in some misleading scenarios. For example, the probability of a pre-damaged structure reaching visible damage (Ψ ≥ 1) approaches 100% as the value of Ψ 0 gets closer to 1, regardless of the PGV. For this reason, pre-damage or undetectable damage has been consistently assigned a value of Ψ 0 = 0.5 and the value showed by the curve of Ψ f ≥ 1 effectively represents the probability of doubling the initial damage value or that of a damage increase (∆Ψ) of 0.5. This means that cracks, initially ten times narrower than the threshold of 0.1 mm, reach or exceed this threshold and become detectable by the naked eye. The influence in the structural behaviour of these micro-cracks was deemed large enough to play a role in the development of further damage and thus impact the probability of displaying light damage. When observing the graphs, one peculiarity stands out: the beginning of the curves, between 1 and 10 mm/s, displays a plateau before sloping upwards. This means that, for low PGV values, a somewhat constant probability of damage is present. This result can be traced to the effect of the model uncertainty parameter ε which is based on the variability observed in the limited number of FE models. These suggest that low values of PGV can trigger visible damage in some cases where material, soil type, and earthquake type are unfavourably combined; especially when some initial damage is present, residual, latent stresses can be liberated by the seismic vibrations regardless of their intensity. This phenomenon requires further study with more analyses of very low vibrations in combination with (multiple) initial damage causes; for the Groningen region, where most buildings are in the areas of these low levels of vibrations, these analyses will be especially relevant.

Comparisons to other studies
It is key to compare these curves, which offer a range for DS1 for Psi between one and two (1 ≤ Ψ ≤ 2), against curves from other authors which are usually defined as a single line indicating the exceedance probability of DS1 or DS2, in many cases, without a clear distinction of whether the threshold corresponds to the incursion or the exceedance of the damage state. Before performing any comparisons with other studies however, several aspects need to be highlighted; chiefly among them are the metrics for damage and the representation of the hazard. The latter can somewhat be translated between studies. This paper has focused on presenting the PGV as a representation of the hazard, that is, the earthquakes. Guidelines and studies on low structural vibrations unrelated to earthquakes employ PGV because it better relates to the energy of the vibrations (SBR, 2017;Rutqvist et al., 2014;Elenas, 2000). Since this study is focused on low earthquake vibrations, it was deemed more appropriate to also make use of PGV as a measure of the hazard so as to make these comparable to other vibration sources such as trains, pile driving and construction activities, and normal usage of the structures. Traditionally, however, the intensity of the earthquake hazard is measured in Peak Ground Acceleration (PGA) as this relates to the forces acting on the structures which is of interest when analysing their strength capacity, something that is not as relevant when considering the light damage of the structure. Furthermore, the pseudo acceleration (Sa) is also used to describe the earthquake hazard since it describes the dynamically-amplified forces on the structure, but because the amplification depends on the natural vibration period of the structures, it becomes inconvenient for producing a single measure of the intensity of the hazard. These three measures can be proportionally related to one another for every earthquake record, but differ between records (and structures). For the comparisons ahead, the average of the transfer factors of the four signals of this study is used to obtain an approximate relationship between PGA and PGV: 26.4 s/s 2 ; for Sa, the upper envelope of the spectra in Fig. 4 is used to establish the ratio between Sa and PGA at different structural periods. Next, the differences in the measurement of damage must be brought forward. This study is unique in that it focuses on quantifying light damage as an explicit, direct function of visible cracks on a masonry wall; this is mostly an aesthetic and durability measure (Korswagen et al., 2019). Most, if not all studies consulted, use a qualitative description of light damage or infer light damage indirectly based on a parameter related to the ultimate (near-collapse) strength of the structures. In this sense, comparing light damage, or DS1, to the outcomes of other studies is difficult. Furthermore, the geometry and material, the structural configuration and connections between structural elements, the type of earthquakes and considerations for the soil and foundations, are all variables that make such a comparison only possible within an order of magnitude. In this context, Fig. 13 attempts this comparison in what is a busy and detailed picture and which is dissected next. The range for DS1, defined herein as an interval between values of Ψ, appears as a shaded area in the figure.
First, in Fig. 13a, the FEMA (FEMA, 2000) guideline specifies drift limits for the DS1 grade; since the limits are suggestions to help the design of a structure prevent the incursion of said structure into the damage grade, for these two curves, assumptions have been made. The limits are presumed to correspond to the fifth percentile, that is, that at the specified value, there is a 5% probability of the pool of structures exceeding DS1, and, that the probability of 50% is reached at twice the drift limit. So, the light damage of infill masonry walls (those with very low overburden and corresponding in the Dutch case to the nonload-bearing façades) seem to be allowed to display light damage with probabilities well within the highlighted area for PGV values up to 60 mm/s. The non-infill structural walls, on the other hand, sport a much lower probability of damage and their curve appears well to the right. Okada et al. (2000) conducted a study of general timber-frame buildings for typical earthquakes; this curve is also closer to the end of DS1 (Ψ≈2.5). This is not unexpected since the unreinforced Groningen masonry has never been designed to resist seismic action and is therefore comparably more vulnerable to earthquakes than the more flexible and purposely-design timber-frame structures. In contrast, Gehl et al. (2013), who studied the effect of different PGA/PGV relations on a masonry façades using the software 3Muri, produces a much steeper curve when employing the ratio of 30 s/s 2 . Similarly, Abo-El-Ezz et al. (2013) studied monumental rubble stone masonry in Canada and found the probabilities of light damage to be just in the middle of the range provided by Ψ = 1 and Ψ = 2.
Of greater interest however, are the empirical DS1 curves of Crowley et al. (2019) which are based on damage data collected in the region up to 2015 and associated to the PGV determined with GMPEs. Two curves are included from this study; one for unreinforced masonry post 1940 and another for the more vulnerable masonry farmhouses. The latter follows well the upper threshold defined by Psi=1 while the former fits precisely in the middle of the DS1 range. Nonetheless, these sections of the curves are log-normally extrapolated from the data points, available only for values of PGV below 20mm/s; for these lower values, the fragility curves determined herein appear to slightly underestimate light damage as compared to the empirical curves. The comparison study looked at plastered walls, which alter the perception of damage to the worse ); yet it is unlikely that this the only reason for the greater damage. Instead, it is probable that the existing old structures already had some initial damage and that the seismic action aggravated it sufficiently to become visible. The corresponding curve for this case is included in Fig. 13c (Ψ0 = 0.5), where the comparison against the curve for low-rise URM post-1940 is much more favourable at the very low values of PGV. Empirical data from farmhouses however, shows a higher probability of damage at these low PGVs. For values above 10 mm/s, the curves diverge and the fragility curves developed herein show a much larger probability of damage, but the range for DS1 still encloses the comparison curves well. The damage scale employed to characterise DS1 and the empirical nature of the data, which is often skewed to overrepresent damaged cases, are also factors to consider in this comparison.
Additional literature for DS2, which should correspond to the end of DS1 given herein, is compared in Fig. 13b. These comparisons are relevant since a larger portion of the set also evaluates Groningen masonry. First, Kallioras et al. (2018) based their results on dynamic, large-scale tests representing calcium-silicate masonry terraced houses. Here, two differences are noted: firstly, this treats terraced houses which usually have a weak structural direction not represented by the two types of façade of this study. Secondly, the calcium-silicate masonry seems to be more vulnerable to light damage than the clay-brick masonry considered (Korswagen et al., 2020). In this context, it is understandable that the curve of Kallioras et al. falls inside and not to the right of the range shaded as DS1 in Fig. 13b. Similarly, Crowley et al. (2019b) produced a DS2 fragility curve for the masonry typology of Kallioras et al. which matches the latter fittingly. However, the shape of the curve for terraced calcium-silicate houses of Crowley et al. is almost perfectly captured by the shape of the shaded area. This gives confidence in the dispersion of the variables used in this study and suggests that the curve of Kallioras et al. should have included more variability. Furthermore, Crowley et al. (2019b) also produced curves for other masonry typologies; an unidentified one has been selected to show that some are much less vulnerable than the terraced houses, which then also fall to the right of the DS1 area.
Moreover, the HAZUS guideline (FEMA, 2013) for analysing (earthquake) risk indicates that slight damage of URM can be expected to develop also within the fragility curves drawn, suggesting that indeed the Ψ≥1 curve indicates the exceedance probability of the start of DS1. Additionally, Van Elk et al. (2019) show an example of a building fragility curve for an undefined masonry typology in terms of the pseudo acceleration (Sa at T = 0.5s). How this compares to PGV is uncertain, since the earthquake signals or spectra of the comparison study are not known; nevertheless, using the envelope of the spectra of the four records (see Figure 4), the relationship between PGV and Sa can be estimated. Then, the masonry system exemplified by Van Elk et al. seems slightly more vulnerable but well inside the DS1 range. Further, Milosević et al. (2018), also explored DS2 fragility for unreinforced masonry with rubble stone and hollow brick. This material is expected to be more vulnerable than solid clay-brick masonry and this is confirmed by the placement of their curve on the left side of the DS1 range. Finally, Waarts (1997)  exceedance of the flexural strength in a probabilistic fashion; their curve also goes through the middle of the DS1 interval but indicates significantly higher probability of light damage at very low PGVs; in Fig. 13c, this difference is smaller. Since Waarts assumed that horizontal strains were present in the structures, the comparison against the curve for Ψ0 = 0.5 seems more appropriate.
In sum, the range of DS1 enclosed by the fragility curves showing the probability of exceeding Ψ ≥ 1 (start of visible damage) and Ψ≥2 (end of DS1, see Table 6) is, albeit the differences in structural typology and masonry material represented by the studies, comparable to other fragility studies looking at the low damage state of masonry. To aid future comparisons and expand the usability of the curves, a fit to a log-normal distribution has been performed and is displayed in Fig. 14; the parameters of the log-normal fit are gathered in Table 10.
While in contrast to other studies, this paper's strength is the more accurate damage definition and realistic monitoring of crack-based damage in complex FEM models that allow crack initiation, propagation and damage aggravation, these computationally-expensive models also limit the number of earthquake motions that can feasibly be analysed and thus, the reduced record-to-record variability included becomes the paper's weakness. Future improvements in computational solving times and modelling strategies will allow the analysis of additional earthquake records and so a better representation of the strong record-to-record uncertainty in the fragility curves.

Conclusions
This study has utilised a five step process to quantify the probability of light damage on bare, fired-clay-brick masonry walls. First, experiments on walls were used to calibrate continuum-based, nonlinear FE models also in terms of crack behaviour (second), which were further extrapolated to analyse the dynamic effect of seismic vibrations regarding the initiation and propagation of cracks (third). The models also contemplated pre-damage situations were the walls had been subjected to differential settlements that led to microcracks or invisible damage. The fourth step regressed relationships between the intensity of light damage and a multitude of parameters varied in the FE models such as the material strength, the initial (pre-damaged) condition, the type of supporting soil, and the Peak Ground Velocity and reiteration of seismic vibrations. The model variations established for instance, that weaker masonry experienced a larger increase in damage and that higher values of PGV led to even larger increases. While these relationships are not unexpected, they could be well quantified. For example, models of walls on peat (soft soil) accumulated 33% more damage than walls on the comparatively stiff sandy soil; while, models of walls with a window attained a 25% higher increase in damage than walls without a window. Similarly, models of walls subjected to two consecutive identical seismic events displayed on average 10% higher damage increase than those subjected to only one event.
As a fifth step, the regressed damage function was used to compute the probability of reaching or exceeding light damage for increasing intensity of PGV in a 'fragility or vulnerability curve', where light damage was quantified as visible cracks at least 0.1 mm in width. The curve shows that at 10 mm/s, there is a 5% probability of light damage in a masonry wall, but, if the wall has sustained prior, yet visibly undetectable damage, in this study induced by differential settlements, this probability rises to 20%. This is a quantitative conclusion about the importance of considering the existing condition of masonry structures when analysing seismic effects, especially for low PGV vibrations. As such, damage claim processes should consider the fact that seismic events not only may induce new damage but also aggravate existing, perhaps previously undetectable damage. In comparison to attempts from other authors to determine similar masonry fragility curves for 'damage state 1', the interval herein defined to correspond to DS1 captures well many of the comparable curves; in particular, empirical curves drawn from claims data match well with the fragility curves established for initial undetectable damage.
While the fragility curves determined herein are specific to the situation in the north of the Netherlands and are constrained by certain limitations such as the type of masonry and structural geometry, the reduced number of earthquake motions, and the type of predamage selected; the approach employed comprising FE calibration, extrapolation, and probabilistic extension, can be used in other situations to determine the probability of light damage.