Stochastic simulation of earthquake ground motions for the seismic assessment of monumental masonry structures: Source‐based vs site‐based approaches

Earthquakes are among the most destructive natural disasters and have resulted in a massive number of fatalities and economic losses all over the world. Simulated ground motion records are valuable, particularly for regions lacking seismic stations or with a limited history of large‐magnitude earthquakes. Notably, a significant percentage of monumental masonry buildings are located in regions with limited access to real records; hence, simulated records play a paramount role in their seismic protection. However, few studies have investigated the structural response of heritage buildings via response history analyses to assess the performance of simulated earthquakes against real ones. To accomplish this, this study simulates the recorded time‐series of the 9th of July 1998 Faial earthquake in the Azores (Mw = 6.2) at four available stations, using two different simulation approaches, that is, a source‐based stochastic finite‐fault method and a site‐based broadband stochastic method. First, two masonry facades with sidewalls characterized by different slenderness levels are adopted to conduct this research. Moreover, the proposed approach is also applied to an existing monumental structure, that is, São Francisco Church, located at Horta, which was affected by damage during the Faial earthquake. Results demonstrate that both simulation approaches provide similar results in terms of structural response prediction. The proposed framework also demonstrates that a small mismatch in terms of predicted damage patterns can result in a significant relative error in terms of displacement predictions.


INTRODUCTION
Recent events, especially the 2023 earthquakes in Turkey, have demonstrated that earthquakes are among the most destructive types of natural disasters, resulting in fatalities and significant structural damage and failures.2][3] When one refers to the preservation of heritage constructions, the probability of reaching a near-collapse state is unpredictable in time because of the intrinsic nature of the geological phenomenon.In some specific cases, authorities and governments are making enormous investments to protect monumental iconic buildings, even though such an approach is not sustainable on a large scale since the massive presence of heritage buildings populating cities and villages. 4In order to define a strategy for architectural heritage conservation, one can refer to studies from the U.S. National Institute of Building Sciences (NIBS) that show how the investment in mitigation saves six times the amount for damage repair ("prevention pays"). 5In this context, today's computational methods offer possibilities for assessing vulnerability and open new perspectives and challenges.3][14][15] Such approaches model the masonry material using different representation scales, that is, equivalent continuum, macro-blocks, or discrete representations and are suitable to perform the seismic assessment of masonry structures.FEM allows a more versatile application, as masonry can be represented either through a continuous equivalent media (designated macro-modeling) or by a discrete representation of units and joints (designated micro-modeling).The macro-modeling strategy is a well-accepted solution since simplifying the problems allows engineers to perform non-linear response history analysis (NLRHA) on large-scale structures.To accomplish this, the demand for a full time-series of ground motion records is an issue, specifically if one refers to a region lacking seismic stations or a limited history of large-and moderate-magnitude and potentially destructive earthquakes.
The paucity of ground motions in many regions, mainly for large-magnitude events at close distances, has motivated the utilization of ground motion simulations.Various simulation methodologies can be categorized as deterministic sourcebased, stochastic sourced-based, stochastic site-based, and hybrid models, each of which has its own limits due to modeling challenges. 16In deterministic methods, the rupture of the fault and the propagation of seismic waves along their path are physically described using a source model (e.g., kinematic or dynamic rupture model) and a material model (e.g., seismic velocity model), respectively.These complex methods can realistically reproduce recorded ground motion waveforms at low frequencies. 17However, stochastic methods are more accurate for high frequencies and represent ground-shaking time series at a specific site.9][20][21][22][23] Site-based stochastic methods implicitly account for these effects by fitting a random process to a recorded time-series of known earthquakes and site characteristics. 246][27] One challenge in the hybrid scheme is merging the two methodologies characterized by different parameter definitions. 28tochastic methods are more prevalent in engineering applications due to their general simplicity in light of the absence of a full wave propagation requiring precise source and velocity models.After calibration of their input model parameters for the specific region of interest, these methods can be powerful tools in assessing the possible effects of future earthquake hazards and disasters. 29,30In this regard, stochastic ground-motion simulation algorithms generate full time series of earthquake motions as alternative seismic inputs.The simulation algorithm requires reliable input-model parameters to generate final shaking, which matches specific features of the region-specific real records.2][33][34][35][36] Before ground motion simulations can be applied confidently in engineering domains, they must undergo thorough and rigorous analysis and validation.8][39][40][41][42][43] In a recent study, the authors used stochastically simulated ground motion records to evaluate the seismic damage to the historic masonry building of Arge-Tabriz (Iran). 44his study aims to compare simulated ground motions from a finite-fault stochastic source-based model 23 and a stochastic site-based model 24 and shed light on their applicability in structural engineering practice in order to perform seismic assessments of historic masonry structures.For this purpose, ground motion stations that recorded on the 9th of July 1998, during the Faial earthquake event with a moment magnitude (M w ) of 6.2, are used.In particular, records at four stations with epicentral distances of less than 150 km are considered.Furthermore, the input-model parameters for both F I G U R E 1 Overall representation of the proposed framework (EQ means earthquake).
approaches are calibrated and verified to minimize the misfit between the real and simulated ground motions.Hence, the goodness of fit (GOF) scores proposed in 45 are evaluated to provide seismologically acceptable ground motion records.Furthermore, masonry facades with two sidewalls (U-shape) characterized by various slenderness levels are adopted as masonry prototypes to conduct this study and reveal the performances of the two simulated approaches with respect to the seismic performance of monumental masonry structures.Finally, the proposed approach is applied to a real monumental structure, the so-called São Francisco Church, located at Horta, which was affected by damage during the 1998 Faial earthquake.
The paper is organized as follows.Section 2 describes the objective and the theoretical framework of the proposed work.The 1998 Faial earthquake and the investigated structures are presented in Section 3. Section 4 presents the results in terms of ground motion simulations performed with both approaches, that is, source-based and site-based.Section 5 discusses the obtained results by comparing the two ground motion simulation approaches with respect to the real records.Finally, relevant conclusions are drawn in Section 6.

FRAMEWORK AND OBJECTIVES
This study's goal is twofold: simulate the 1998 Faial earthquake with alternative stochastic ground motion simulation approaches and determine the impact of these approaches on the structural assessment of historic masonry structures.To accomplish these objectives, the framework described below has been followed: • Collection and processing of recorded accelerograms at available stations.
• Calibration of input-model parameters corresponding to stochastic ground motion simulation approaches.
• Validation of simulations is twofold: by comparing simulated and real earthquake records in time and frequency domains, and by evaluating GOF scores.• Implementation of finite element (FE) models of monumental masonry buildings to perform reliable non-linear response history analysis (NLRHA).• Definition of a framework for the comparison of the structural response for masonry structures applicable for both farand near-field records.
Figure 1 is a holistic representation of the proposed framework.Key aspects and theoretical background are detailed and described in the following subsections.

Ground motion simulations and validations framework
The following subsections describe the theoretical frameworks of two simulation approaches to generate artificial records.
In addition, the mathematical formulation to derive the GOF scores is presented.

Stochastic source-based approach
The stochastic source-based approaches encompass both point-source and finite-fault methods.The initial version, which was a point-source approach, was initially introduced by. 180][21][22] This study employs the most recent version of the stochastic finite-fault ground motion simulation methodology within the EXSIM12 platform. 46The method is modified from the algorithm suggested in, 22 which was developed based on the original FINSIM code in, 21 integrating the improvements proposed in. 20It is noted that the modified version of the stochastic finite-fault approach improves the low-frequency portion of simulations.This algorithm identifies the fault rupture by parameters such as earthquake magnitude, fault geometry, strike, dip, and slip distribution.The source contribution is incorporated with the path propagation and site effects to reach the seismic signal in the time domain at any observation site.
In the stochastic finite-fault ground motion simulation approach, the ruptured fault plane is expressed as the grid of smaller sub-sources by assuming a point-source for each sub-source with an  −2 source spectrum. 18,19The rupture of each sub-source occurs with an appropriate time delay contingent on the distance of the sub-source from the hypocenter.Summation of the delayed sub-source contributions is accomplished in the time domain as follows: where () is the total seismic signal at a time ,  corresponds to the total number of sub-sources,   represents the seismic signal of i th sub-source, which is the inverse Fourier transform of its spectrum, 19,20 Δ  is the sum of the fracture initiation and time delay due to the distance of the i th hypocentral distance, the term   corresponds to the fraction of rise time considered for additional randomization, and finally, the term   resembles the normalization factor of the i th sub-source introduced for the conservation of energy.The normalization factor is calculated as follows: where  0 stands for the corner frequency of the entire fault plane,   is the j th frequency ordinate,  0 is the total seismic moment, and the terms  0 and  0 are, respectively, the seismic moment and corner frequency of the i th sub-source formulated as follows: 0 = 4.9 ⋅ 10 6 ⋅   ⋅ where   represents the slip of the i th sub-source in Equation (3).In Equation (4), the term   represents the total number of sub-sources that are activated when the i th sub-source triggers and Δ is the stress drop in bars.The term PP is the pulsing percentage.The algorithm is based on a dynamic corner frequency approach where the corner frequencies of the activated sub-sources descend with rupture progress until reaching a specified level, which is PP for the rest of the sub-sources, and the corresponding corner frequency remains constant.

Stochastic site-based approach
The stochastic site-based model in 24 modulates and filters a stochastic white-noise process using a linear time-variant system to describe a ground motion record at a site.This procedure represents certain recorded ground shaking time-series characteristics instead of modeling the earthquake source and path effects.The stochasticity comes from the white-noise process, while the temporal and spectral non-stationarities are constructed by modulation and filtering of the input process.The form of the model is expressed as follows: where () is the acceleration time series, (, ) is a deterministic time modulation function that shapes the temporal evolution, ℎ( − , ()) is the impulse response function of the filter, () is the white-noise process, and  ℎ () is the standard deviation of the integral process.
In accordance with, 24 a piecewise function that describes the temporal variation of the waveform is employed for the time modulation as follows: where the process starts at  0 and experiences strong shaking phase between times  1 and  2 .The intensity of the strong shaking phase is denoted by the coefficient  1 , while the shape of the decaying end of the function is controlled by the coefficients  2 and  3 .
The pseudo acceleration response of a single-degree-of-freedom linear oscillator is selected for the filter function, which contains time-varying parameters.The filter function is given as: where  and  are the frequency and damping of the filter.A linear function is considered for the variation of the frequency parameter as follows: Additionally, the following piecewise constant function is selected for the damping ratio that accounts for the observed variation of bandwidth of the real records: where  1 ,  2 , and  3 are constant damping ratios,   1 and   2 are the times at which the damping ratios change, and   is the total duration of motion.
To identify the model parameters for a given accelerogram, the evolving intensity and frequency contents of the stochastic process described in the Equation ( 5) is fitted to those of the recorded ground motion in accordance with the objective functions reported in. 24This involves matching the cumulative expected energy of the stochastic process with that of the recorded ground motion, which yields the six parameters of the modulation function.Similarly, the cumulative expected number of zero-level up-crossings of the stochastic process and the average cumulative number of positive minima and negative maxima of 10 simulated processes are fitted to those of the recorded ground motion, resulting in identification of the seven parameters of the filter function.

Validation approach
Simulated motions need to be validated qualitatively and quantitatively against real records to be seismologically valid.On the other hand, validation is application specific, and depending on the application of records in engineering, synthetics may accurately represent observed ground motions only in terms of a few characteristics of the waveform relevant to the application.To this end, the time and frequency content of recordings can initially be analyzed qualitatively.For the statistical evaluation of simulations as the next stage, GOF measures are relevant.One of the well-known GOF-based validation approaches is proposed in 45 and is based on the complementary error function of the seismological parameters using the following equation: where R sim and R real are simulated and real ground motion parameters, respectively.A number of seismological characteristics, such as Fourier amplitude spectra (FAS), pseudo-spectral acceleration (PSA) for 5% damping ratio, and other amplitude, energy, or duration-related intensity measures, are candidates of metrics to evaluate the mismatch between the real and simulated time series.The classification for the GOF score is as follows: In this study, to validate the simulated records of the 1998 Faial earthquake, a combination of a set of seismological measures representing peak ground motion responses, Fourier and response spectra, total energy, and shaking duration are considered.These metrics include peak ground acceleration (PGA), peak ground velocity (PGV), peak ground displacement (PGD), PGV/PGA, Arias intensity (I a ), cumulative absolute velocity (CAV), acceleration spectrum intensity (ASI), modified acceleration spectrum intensity for the period range of 0.1-2.5 s (ASI*), 47 velocity spectrum intensity (VSI), Housner intensity (HI), bracketed duration (tb), FAS within the frequency range of 0.1-25 Hz, and PSA within the period range of 0-4 s.It is important to acknowledge that the classification of GOF scores from excellent to bad fit is subjective, and the ranges can be subject to alteration based on the intensity measures employed.As mentioned earlier, this study incorporates various sets of ground motion parameters to address this concern.

Structural analysis framework
The following subsection presents the numerical framework adopted to perform the NLRHAs.Specifically, details concerning the FE implementation (i.e., material constitutive relationship, damping, mesh discretization, and solution algorithm) and the structural response measures (SRM) for comparing results are presented.

Numerical finite element modeling
Geometrical modeling of the investigated structures is developed using the software Rhino 3D, since it allows one to efficiently represent spatial geometries through a non-uniform rational B-spline (NURBS) approach. 48The FE model is generated in the advanced commercial software ABAQUS CAE. 49he so-called macro-modeling approach is followed, meaning the masonry arrangement is smeared over a homogeneous material.This is particularly convenient for the analysis of large-scale structures 10 since it offers the possibility to reproduce the macroscopic masonry mechanical behavior through several models, for example, the smeared crack concrete, the brittle crack concrete, and the concrete damage plasticity (CDP) models.Specifically, CDP couples plasticity with a scalar-based damage model, and, as it was originally developed for concrete, an isotropic elastic behavior is assumed. 50This is a limitation when adapting CDP to masonry, as anisotropy may have an important role, especially when periodic or non-periodic masonry arrangements are present.However, it is worth remembering those anisotropy directions are hard to find in large-scale structures as many masonry arrangements may exist concurrently. 51,52The importance of anisotropy is thus lessened, and isotropic behavior is an acceptable simplification.The CDP material model accounts for different strength values, post-peak behaviors, and damage descriptions for tensile or compressive regimes.It has been extensively used for the study of large masonry structures, 7 and the results indicate that it offers a good compromise between computational time and accuracy.
The quasi-brittle nature of masonry is represented by a linear type of softening in tension.In compression, a plateau exists after the compressive strength, followed by a linear type of softening.Damage variables are adopted when softening is active and aim at reducing the initial (undamaged) elastic modulus through the following equations: where  0 is the elastic modulus of the undamaged masonry,   is the effective stress value;   is the damage parameter relating the effective stress with the corresponding inelastic strain,   is the total strain value, and    is the inelastic (plastic) strain value.The subscript  reads as  or , if associated with the compressive or tensile regime, respectively.A scalar-based damage model describes the damage in tension   (cracking) and compression   (crushing), with a value between zero (no damage) and one (fully damaged).
When cyclic loading is applied, CDP assumes a non-associative flow rule given as a Drucker-Prager hyperbolic function and requires the definition of several physically based parameters.In ABAQUS CAE, the Drucker-Prager strength domain criterion is implemented, and its parameters are summarized in Section 3.2, which details the specific structures.
Damping is a further important factor that influences masonry structures' seismic performance.In historical masonry, some authors performed calibration of the damping ratio, which is attested to be approximately 3%. 53Other authors, (i) knowing that the damping ratio of masonry structures that experienced damage is greater than that in the undamaged state and (ii) assuming a contribution of hysteretic damping to the total damping of the structure, have adopted a damping ratio of 5%. 44In terms of damping modeling, in this work, the Rayleigh damping model has been adopted.The Rayleigh model approximates the damping coefficient as a linear combination of mass and stiffness: where , , and  represent the mass, stiffness, and damping matrixes, respectively.The coefficient  and  denote the mass-and damping-proportional coefficients, respectively.These coefficients are computed as follows: where  is the damping ratio, and   and   are the frequencies of two selected modes.According to, 44   is taken as the first natural frequency in the considered direction, whereas   is the frequency of the highest mode with a significant participation factor.The three-dimensional FE model enforces the use of solid elements; therefore, the mesh discretization is achieved using tetrahedron (Delaunay) FEs, due to their adaptability to more complex geometries.The so-called TETC3D4 FEs in ABAQUS CAE, 49 based on a tetrahedral geometry with linear interpolation, have been used.
Furthermore, appropriate boundary and loading conditions are implemented in order to perform the structural simulations.NLRHA is performed by applying a real or simulated record.Two phases are idealized in the loading process: initially, the gradual application of gravity loads, and subsequently, the record is applied to the base of the structure, separately along the two principal geometrical directions, namely X and Y.In order to integrate the equations of motion, an implicit time integration scheme has been adopted with non-linear geometric effects taken into consideration. 49 I G U R E 2 Macro element damage pattern comparison for a masonry portal with an opening: Proposed method (the shaded blue and red areas indicate, respectively, the damage pattern generated by two different records, e.g., real and artificial ones).

Comparison framework
Beyond a reliable FE model to perform high-fidelity NLRHAs, an appropriate comparison framework that considers several structural response measures (SRM) must be defined.This subsection briefly describes the methodologies used to treat the SRMs obtained by simulating the seismic behavior of the masonry structures via both real and simulated (i.e., source-based and site-based) records.Relative errors in the prediction of SRMs are calculated with respect to the response obtained with the real records.The following response measures are considered: • Maximum base shear.
• Maximum horizontal displacement (a control-point network is considered).
Base shear and displacement as response measures arise because it is common engineering practice to assess seismic performance according to force or displacement capacities.On the other hand, as supported by the literature, accelerograms' characteristics may also affect the damage pattern of unreinforced masonry structures.To this end, it has been assumed as a third response measure.
Regarding the first two comparisons, the relative errors are computed with the following equations: where  and  refer to a specific SRM obtained using simulated and real records, respectively.One should note that maximum base shear and maximum control-points displacement comparisons might be performed even if the structure under investigation does not present damage, whereas the damage pattern becomes relevant only when non-linearities are experienced.Different approaches to quantify the damage on masonry structures have been proposed in the literature [54][55][56] ; however, no studies assess the similarities in terms of damage semiotics, giving them an integral score.To this end, a methodology involving the Hausdorff distance score method 57 has been integrated with concepts arising from (i) the structural damage occurring in masonry structures and (ii) engineering-based knowledge from several post-earthquake field observations of the authors of this paper.
The proposed approach consists of discretizing the structure in macro elements, for example, piers and spandrels, and analyzing them independently.The crack patterns of the macro element Ω  are represented in Figure 2, where the damaged area resulting from the reference set () and the simulated one () are represented.
The Hausedorff distance method [58][59][60][61] is applied to each point belonging to the damaged area of one set (or crack line) to find the minimum distance with respect to the damaged area of the other set.In contrast to the traditional use of the Hausdorff distance technique, which involves measuring the maximum distances; this approach calculates the average distance considering all the pixels belonging to the damaged area.This is done to avoid potential errors when interpreting damage near the edges of the macro elements.Additionally, only the damaged area that is not overlapped with respect to the second set is considered (see Figure 2).Such a step is performed in both directions, that is,  and .In particular,  indicates an array collecting the Hausdorff distances of all points belonging to the cracked area of the set  with respect to the set . Hence, for the sake of clarity,    indicates the shortest distance between the non-overlapping points affected by damage  belonging to the set  and the crack pattern  (see Figure 2).The average value can be computed as follows: At this stage, both   distances are normalized by the maximum value of the Hausedorff distance computed at each point belonging to the macro element domain Ω  with respect to the reference set, that is,  (see Figure 2).Such a distance is represented with the symbol  max  : Finally, the score for each macro element can be computed via the following equation: where, Ω  and Ω  indicate the non-overlapping damaged area (or crack line), whereas Ω   and Ω   are the total damaged areas.Both quantities are referred to as the sets  and , respectively.One should note how two perfect equal crack patterns will assume a value of 1 ( = ), whereas lower values indicate more pronounced dissimilarities (assuming at most a null value).Equation ( 17) provides a score referred to the single macro elements.However, such an approach can be further mathematically manipulated to provide a measure of the full structure or of specific groups of macroblocks, as: where  corresponds to the number of considered macro elements and Ω  is the area of the  − ℎ macro element.To study masonry structures that are likely to be affected by localized failure mechanisms, Equation (18) assumes particular relevance for defining a meaningful score when isolating that specific part.

DESCRIPTION OF THE CASE STUDY: EARTHQUAKE AND MASONRY PROTOTYPES
The next subsections aim to describe the earthquake event used as the case study.Hence, the adopted masonry prototypes' geometry and mechanical features are detailed and described.

9th July 1998 Faial earthquake (M w = 6.2), Azores, North Atlantic
The Azores plateau, with three major tectonic plates, is located at the North America, Eurasia, and Nubia triple junction 62 (Figure 3).The plateau includes Faial, Pico, São Jorge, São Miguel, Santa Maria, Graciosa, and Terceira islands, which experienced several earthquakes during their history: 1926 with a body-wave magnitude (M b ) ∼5.6, 1958 with the maximum modified Mercalli intensity (MMI max ) of X, and 1998 with M w of 6.2 that struck Faial Island.The 1998 Faial earthquake occurred between the islands of Faial and Pico in the Azores plateau.The earthquake was felt in almost all Azores islands and caused predominant structural damage in Faial and Pico with an MMI max of VIII. 63The number of fatalities, injured people, and homeless was reported to be 8, 150, and 1500, respectively. 64Additionally, 2100 buildings experienced partial or total collapse. 64The earthquake was recorded at five strong ground motion stations, as illustrated in Figure 4.A literature survey reveals the existence of different source models in the region (FL1 and FL2).In, 65 fault line 1 (FL1 in Figure 4) has been proposed as ruptured, whereas 66 investigated various fault plane solutions that could rupture during the 1998 Faial earthquake.The findings suggested that the ruptured fault plane is the dextral strike-slip east-north-east (ENE)-westsouth-west (WSW) solution with strike and dip angles of 264 • and 83 • , respectively (FL2 in Figure 4).Table 1 lists detailed information in terms of epicentral distance (R epi ), coordinates, PGA, and PGV of the stations that recorded the 1998 Faial event for the east-west (EW) and north-south (NS) horizontal components (esm-db.eu). 68Among the stations, the highest shaking level was recorded at the near-field station Horta (HOR) with an epicentral distance of 11 km, as recorded by the European strong ground motion database (ESD). 69At the other far-field stations, lower shaking  levels were observed.For ground motion simulation, this study considers the recorded motions at the four stations, Angra do Heroismo (GZC), Praia Vitoria (PVI), HOR, and Sao Sebastiao (SEB), with less than 150 km epicentral distances.

Masonry prototypes: Geometrical and mechanical features
As mentioned above, two masonry facades with sidewalls (U-shape) that differ in height are considered.Figure 5A represents the overall characteristics and geometry of both masonry prototypes.These structural benchmarks are an abstraction of single-nave churches, frequently incurred in European earthquake-prone regions and reported as one of the most damaged masonry typologies in post-earthquake surveys. 14,70Here, it is noted that size and slenderness play a major role in rocking dynamics and OOP failure.Finally, the São Francisco Church, located at Horta, is analyzed since the authors own a detailed geometrical survey performed post-earthquake (Faial 1998), and the same was damaged by the earthquake and subsequently strengthened (Figure 5B).As mentioned in Section 2.2.1, the CDP model is used to model the non-linear behavior of the masonry material.Elastic behavior is described by elastic modulus and the Poisson ratio, which are set equal to 1000 MPa and 0.2, respectively.The U-shape prototypes' mechanical properties are summarized in Tables 2 and 3. A Rayleigh damping factor equal to 3% has been considered.These values have also been utilized as a starting point for calibrating the procedure of the real case of study, that is, the São Francisco church, since they are similar to those adopted in the technical report developed for the post-earthquake assessment of the church. 71However, even if out of the scope of this paper, the calibration procedure has regarded the damping factor in order to reach damage of the numerical model reasonably close to those surveyed post-earthquake.The damping factor has been calibrated to a value of 5%, whereas compressive and tensile non-linear behavior is the same as those adopted for the U-shapes reported in Table 2.

TA B L E 2
Compressive and Tensile behavior of the masonry: U-shape configurations.

Compressive behavior Tensile behavior Stress [MPa]
Inelastic In this study, the recorded ground motion time-series of the 1998 Faial earthquake are simulated using two stochastic approaches: the source-based stochastic finite-fault ground motion simulation approach in, 46 which is a modification to the initial version of, 22 and the site-based stochastic approach proposed in. 24The real records are taken from the ESD website. 69The real and simulated accelerograms for the four stations of interest are baseline corrected and bandpass filtered between 0−25 Hz.The output of the source-based stochastic finite-fault method is a single representative (random) horizontal component at each station.On the contrary, independent simulations for the site-based stochastic method are conducted for each horizontal EW and NS component at the considered stations.There are a total of eight real timeseries at the four selected stations, in addition to the 12 simulated time-series from both methodologies.In the following subsections, the calibration of input-model parameters for simulations based on the approaches provided in Sections 2.1.1 and 2.1.2is discussed in depth.

Input-model parameters (source-based)
The 1998 Faial earthquake was previously simulated in 31 at only near-field station HOR on the bedrock.One should note how those simulations have limitations in reflecting site characteristics.In this work, simulations are conducted by incorporating site effects and region-specific input parameters in order to model and enhance the outcomes at all recorded stations realistically.In parallel, for the sake of comparison, records are also simulated based on the input parameters provided in, 31 updated to account for site effects.A review of the relevant literature in the region demonstrates the prevalence of several seismic sources and attenuation models.The ruptured fault plane employed in 31 (FL1 in Figure 4) has a length and width of 16.5 and 9.4 km, respectively, estimated according to. 72The strike and dip angles by that study are provided as 165 • and 85 • , respectively.Later, in 66 potential fault plane solutions that could rupture during the 1998 Faial earthquake were studied.The results of that study proposed the ruptured fault plane as the dextral strike-slip (SS) ENE-WSW solution with length and width of 12 and 5.5 km, respectively (FL2 in Figure 4).The strike and dip of the fault plane solution in 66 are given as 264 • and 83 • , respectively.Two alternative geometric spreading models are used and tested for the path effects.The first model was employed in. 31ater, a region-specific model was proposed in. 67In, 73 a model for South Iceland for the quality factor was proposed.Due to having a similar geodynamic setting to the study area, 31 employed this model for the simulation process.Later, in 67 a region-specific model for the Azores region was proposed.Here, two alternative sets for the input-model parameters are tested and calibrated to simulate the records at the stations, that is, Set 1 and Set 2. The records of Set 1 are simulated based on the model proposed in, 31 modified for the site effects.On the other hand, the records of Set 2 were generated based on the region-specific input-model parameters.For both sets, site responses were modeled through the generic soil amplification factors proposed in. 74Soil classes at the stations are taken from. 75To this end, at stations HOR and SEB soil class D is employed, at station PVI generic soil is employed, and at station GZC soil class C frequency-dependent amplification factors of 69 are employed.Finally, the kappa values are taken from the regional study of. 67Verified input parameters used for modeling the two alternative sets of simulated time-series at the selected stations are listed in Table 4.

Input-model parameters (site-based)
The parameters of the time-modulation function and the frequency and damping of the filter describing the impulse response function are derived by fitting the stochastic model to the recorded ground motion time-series at each station.The modulation function is characterized by matching the expected cumulative I a (Arias intensity) of the model to each record.For the frequency and damping, the cumulative number of zero-level up-crossings and the average cumulative number of positive minima and negative maxima of 10 simulated processes are fitted to those of the recorded ground motions in accordance with the objective functions described in. 24The fitted model parameters are given in Tables 5  and 6.

Results and discussion
This section displays and compares graphically and statistically the real and simulated time-series from two methodologies, with reference to the GOF validation approach outlined in Section 0. The simulated motions are initially validated against recorded data statistically to be accepted by engineering practitioners.The simulations from the source-based stochastic method are compared against the geometric mean of the real horizontal components.A comparison of the GOF scores for the simulated records generated using Set 1 and Set 2 alternative input-model parameters introduced in Section 4.1 demonstrates that Set 2 yields better outcomes, resulting in higher GOF scores for all stations (Table 7).
Statistical results from the source-based approach based on input-model parameters of Set 2 demonstrate the validity of simulations in the category of very good fits compared to the real records from different seismological aspects mentioned in Section 2.1.3(Table 7).The only exception is the simulated time-series at station SEB, which has a fair fit.Consequently, the results of Set 2 are verified and analyzed further below.For the site-based stochastic simulations, a component-wise comparison is conducted (in each direction separately).The GOF scores for the simulations show an excellent fit against the observed records for all stations and components except for the EW component of station GZC and the NS component of station HOR.This is expected, given that the approach relies on a fitting procedure to simulate characteristics of recorded motions as accurately as feasible.Consequently, the input parameters employed for the simulation of the 1998 Faial event at the selected stations are validated.
Results for the real ground motion records besides simulations based on both stochastic approaches in terms of full time-series, FAS, and PSA at the selected four stations are shown in Figures 6-9.It is noted that the results of the stochastic source-based approach are only provided for Set 2 as the best model.Among the selected stations, station HOR has the closest source-to-site distance with the shortest recorded duration of ground shaking.At this near-field station, the maximum horizontal observed PGA is recorded as 435 cm/s 2 .PGA of the simulated records from source-based and site-based approaches are estimated as 417 and 488 cm/s 2 (in terms of the maximum of the EW and NS components), respectively.Overall, the spectral contents of the simulated motions are close to the observed ones for the frequency range between 0.1 and 25 Hz.Station GZC has a recorded PGA of 17 cm/s 2 , while the source-based and site-based methods simulate PGA of 16 and 15 cm/s 2 .The source-based simulation slightly overestimates FAS within a frequency band of 0.4-1.0Hz, which is also obvious in the spectral ordinates.Station SEB has an observed PGA of 21 cm/s 2 , while the corresponding source-based and site-based methods simulate PGA values of 11 and 22 cm/s 2 , respectively.At this far-field station, the source-based approach underestimates the PGA level.
Additionally, the simulated FAS from the source-based approach underestimates the observations for frequencies above 1.0 Hz.The underestimation of low-period content might be attributed to the insufficient modeling of soil effects at this station.In addition, the observed PSA at station SEB shows multiple dominant periods that are not captured by the one-mode filter of the site-based method.This is due to the limitation of using a single degree of freedom impulse response function as the filter.The recorded at station PVI is 10 cm/s 2 while both stochastic simulations result in a PGA of 8 cm/s 2 .At this station, the simulated FAS closely resembles the real FAS, except for a slight overestimation of the amplitude at frequencies below 1.0 Hz, which is visible from PSA.
At near-field station HOR, the maximum observed PGV is 37 cm/s, while lower values of 24 and 29 cm/s are simulated using source-based and site-based methods.In general, for far-field stations, the source-based method results in a slightly higher discrepancy in simulated PGV compared to the site-based method.The maximum observed PGVs in stations GZC and PVI are 1 and 0.8 cm/s, respectively.For these stations, the site-based method estimates PGVs of 0.9 and 1 cm/s, while the source-based method results in higher estimates of 1.6 and 1.3 cm/s, respectively.In contrast, at the farthest station SEB, the observed maximum PGV is 2.2 cm/s, but the source-based method underestimates the simulated PGV at 1.3 cm/s, while the site-based method estimates it to be 1.9 cm/s.In contrast, the site-based method yields larger discrepancies in terms of PGD compared to the source-based method.The maximum observed PGD in station HOR is 5.5 cm, whereas the source-and site-based methods yield simulated PGDs of 4.0 and 10.0 cm, respectively.At stations GZC and PVI, the maximum observed PGDs are 0.55 and 0.3 cm, respectively.The source-based method estimates PGDs of 0.5 and 0.4 cm, while the site-based method results in lower estimates of 0.25 and 0.35 cm, respectively.Finally, at station SEB, the observed maximum PGD is 0.45 cm, and the source-and site-based methods estimate simulated PGDs of 0.4 and 0.55 cm, respectively.
Overall, results at all stations reveal that two alternative sets of simulated ground motion records are seismologically acceptable since the records are in agreement with the real records in approximately the entire time, frequency, and period domains.In this study, the input-model parameters of the site-based approach are calibrated for the full ground motion waveforms, resulting in more realistic accelerograms, including duration and shape.However, the source-based approach can only simulate the shear-wave (S-wave) portion of the ground motion record, limiting its ability to capture the P and surface waves accurately.Consequently, the duration and shape of the simulated motion may not be consistent with the full waveforms of the real records.Nevertheless, incorporating region-specific source, path, and site parameters in the simulation process allows the source-based approach to realistically capture the S-wave part of the motion, which is more relevant to structural damage.The records are then further examined for engineering purposes in terms of estimating the seismic response of various masonry structures.

RESULTS OF THE NON-LINEAR RESPONSE HISTORY ANALYSES
In this section, the results of the NLRHAs are critically discussed.In particular, the framework introduced in Section 2.2.2 has been adopted to compare the SRMs obtained from the simulated site-and source-based accelerograms and assess their performance with respect to the structural responses using real records.At first, two masonry facades with sidewalls (U-shape) characterized by different slenderness levels are considered.Finally, São Francisco Church, located at Horta close to the HOR station, is analyzed.

U-Shapes with different slenderness
The overall characteristics and geometry of both masonry prototypes, that is, P1 and P2, and the numbering of the control points adopted are reported in Figure 5A.The same mechanical properties have been assumed for P1 and P2, whereas a variation in terms of facade height is adopted, resulting in different fundamental periods T 1 , that is, T 1 = 0.20 s for P1 and T 1 = 0.28 s for P2.
As mentioned above, the source-based approach generates only one component, whereas two components per record characterize the real and site-based datasets.Since the ground motion is applied in the out-of-plane direction of the facade at both polarities, that is, ±Y, real and simulated site-based record components are run separately.One can note that the far-field stations have relatively low PGA, ranging from 10 to 22 cm/s 2 , generating linear elastic responses of the masonry structures, which do not experience any damage.On the other hand, the near-field station is characterized by higher PGA, which is expected to generate damage (see Section 4).However, for the sake of comparison, the near-field records have also been used to run response history analyses with a linear elastic constitutive material relationship to quantify how the appearance of damage affects the accuracy of the two investigated ground motion simulation approaches, that is, site-and source-based.The total number of NLRHAs performed is equal to 100.
The results obtained by the source-based simulation approach have been used to compute two relative errors, one with respect to each component of the real records (NS and EW).On the other hand, the site-based relative errors have been computed with respect to the actual component simulated.
At first, the stockier U-shape benchmark is discussed, namely P1 (Figure 5A). Figure 10 presents the relative errors between the simulated and real records in terms of maximum displacement (  ) experienced during the time history for each control point (1 to 15). Figure 11 presents the errors computed analogously, though in terms of the base shear (  ).The P1 prototype presents an overall balance between underestimations and overestimations with relative errors that show good accuracy of the two simulation approaches.It is worth remarking that such a trend is recorded for both polarities, ±Y.
Concerning the near-field station, the simulated records again provide a consistent relative error magnitude and an overall balance between underestimations and overestimations.Observing the relative errors in terms of base shear prediction (Figure 11), one should note how underestimation and overestimation trends of the simulated records are reflected for the far-field stations.In contrast, some inversion of polarities is observed for the near-field station, creating doubts about the accuracy of such approaches to simulate near-field records.More significant errors in displacement prediction have been obtained when a non-linear constitutive behavior of the material has been adopted.However, such a result was expected since the material and geometrical non-linearities affect macro elements' local stiffness and amplify the predicted error, particularly when no perfect matching of the damage pattern is obtained between real and simulated records.In order to verify this, the damage crack pattern comparisons for the near-field simulations showing the greatest error are represented in Figure 12.The results also underline how a relatively small damage area can drastically drop the prediction accuracy of the site-based simulation approach.A similar discussion can be drawn for the slender prototype, namely P2 (see Figures 13 and 14).At the far-field station, the predictions of both site-based and source-based simulations are slightly more accurate because the structure presents a higher fundamental period.Furthermore, more consistency is visualized concerning the near-field simulations performed using linear and non-linear material constitutive relationships.Such an observation suggests more consistent results in damage pattern prediction, which is underscored afterward by applying the modified Hausdorff score (see Section 2.2.2).On the other hand, compared with P1, predictions show less accuracy in simulating displacement of the control points for near-field records, even in the case of linear elastic constitutive relationships.For clarity, a brief summary of the relative errors obtained for each part (i.e., sidewalls and facade) is reported in the Appendix for both P1 and P2 prototypes.
In this study, it is observed that near-field site-based simulations have larger errors in terms of maximum displacement (as shown in Figure 10).This can be attributed to several factors.Firstly, the correlations between the two components are not accounted for in the simulations.It is important to note that these correlations are more significant in nearfault motions.Additionally, directivity pulses are not accounted for, which may have contributed to the larger errors observed in the simulations.Interestingly, these shortcomings do not result in a larger error in the case of base shear.The same observation is valid for source-based simulations.This demonstrates that the validation of simulations is highly application-specific and that the accuracy of simulations may vary depending on the specific parameters being examined.
Finally, Figures 15 and 16 represent the  score, introduced in Section 2.2.2, which is evaluated for each macro element of P1 and P2 prototypes, where the nomenclature for each macro element is reported in Figure 17.A global score is also computed by considering the weighted value with respect to the area of the single macro elements (see Equation ( 18)).Overall, results show better agreement for the P2 case.Regarding the P1 +Y case, the damage patterns generated by the simulation approaches are in good agreement with those of the real records at the Facade and Sidewalls-2 level.On the other hand, lower accuracy has been obtained at the Sidewall-1.By inverting the polarities, the results slightly change.The damage pattern score at the facade level drops, whereas a better agreement has been obtained at the sidewalls.Such a result also agrees with the relative error displacement reported in Figure 10, where the pattern is not consistent when inverting the polarities.One should note that, overall, the Site 2 -NS comparison is the worst, as reflected in Figure 10.The P2 prototype overall has obtained better correlation scores with some exceptions for macro elements 2 and 5.

São Francisco Church
The FE model of São Francisco Church and the numbering of the control points monitored (1 to 9) are represented in Figure 5B.One-dimensional NLRHAs have been performed by applying the ground motion in turn along the X and Y directions at both polarities.The total number of simulations performed is equal to 100.Hence, the simulated records' performance is now tested considering a real church, to assess the consistency of the results obtained for the U-shapes and draw preliminary conclusions regarding the reliability of such approaches in engineering practice.Figures 18 and 19 report the relative errors (obtained with the two simulated approaches) of control point displacements and base shear with respect to the NLRHAs performed with the real ground motions.However, referring to the control point displacement predictions only, slightly smaller errors are registered, with some peaks rarely exceeding 50% and mainly in the near-field simulations considering a non-linear constitutive relationship.The comparison also confirms the same trend in base shear prediction, where errors appear within a reasonable range of accuracy, underlining the good performance of the simulated ground motions.However, referring to just the near-field simulations, what attracts attention is the trend to underestimate the prediction, particularly when non-linear constitutive relationships have been used.Such results open new perspectives and the need to perform further studies assessing the reliability of the simulated ground motion time-series as a function of the structural typology, both for progressive damage and limit condition states in international codes, such as damage limitation, significant damage, or near collapse.

CONCLUSIONS
This manuscript aims to simulate the 1998 Faial earthquake with alternative stochastic ground motion simulation approaches and determine the impact of these approaches on the structural assessment of historic masonry structures.Both seismological approaches are calibrated and tested against the recorded time-series of the 9th of July 1998 Faial earthquake (Mw = 6.2) at four available stations.Overall, the stochastic site-based approach yields superior results over the source-based approach in terms of higher goodness of fit (GOF) scores.The source-based stochastic method necessitates more specific information regarding the source, path, and site effects.On the other hand, the site-based stochastic method is more straightforward, practical, and computationally efficient, though it is record-dependent, and the model parameters are calibrated based on the actual motions.Findings reveal that, with the source-based technique, the results are closer to the real records in the time and frequency domains when regional input-model parameters are employed as opposed to global parameters.This observation is also valid in terms of the GOF values in which higher scores are obtained from the region-specific input-model parameters.
The site-based method yields comparable GOF scores for both near-field and far-field stations.The explanation could be linked to the methodology's inherent fitting technique.In contrast, the source-based method yields the highest GOF for the near-field station, indicating that the input-model parameters have been better calibrated.The worst match is obtained at the most distant station, possibly due to a poor attenuation or inappropriate soil model.The source-based method allows the remaining stations' GOF scores to fall within the same range.Due to the time domain simulation, the duration content of the records simulated using the stochastic site-based technique is more representative of actual motions than the source-based method.Furthermore, the source-based method only simulates the S-wave component of the ground motion records.
The results of the study suggest that the larger errors observed in the near-field site-based simulations in estimating maximum ground displacement can be attributed to the lack of accounting for correlations between components and the absence of directivity pulses in the model.The correlations between the two components are found to be significant in near-fault motions, and their neglect can lead to larger errors in predicting the maximum ground displacement.In addition, directivity pulses, which are caused by the rupture of the fault, are not accounted for in the simulations and could also have contributed to the observed errors.
Both stochastic simulation approaches are used at the structural level to perform non-linear response history analyses (NLRHAs) of adequately selected masonry prototypes.In order to assess their accuracy, a proper comparison framework accounting for structural response measures is introduced.Such a tool permits comparing the simulated damage patterns against the real ones (i.e., obtained from the real records) via the introduction of an objective measure (the modified Hausdorff distances score method).This measure involves concepts related to the structural damage occurring in masonry structures and engineering-based knowledge obtained from several post-earthquake field observations.NLRHAs underscore how the appearance of localized structural damage affects the prediction of the other response measures, that is, displacements and base shear, emphasizing the need to consider the damage pattern in the comparisons.
To conclude, the accuracy of simulations can vary depending on the application, and the study stresses the need for further research to explore the effects of correlations and directivity pulses in non-linear dynamic analysis results, particularly for near-field site-based simulations.Such research would help to refine simulation models and improve the accuracy of predictions in earthquake engineering.As a limitation of this study, it is important to note that the investigation focuses solely on the Faial 1998 events and utilizes a limited number of records.For the sake of the general applicability of the proposed methods, it is recommended that future investigations employ the framework proposed herein to analyze a broader range of events, specifically including a larger number of near-field stations.In addition, in an initial endeavor to validate the suitability of the stochastic site-based approach for predicting the response of masonry monuments, this study involves simulating various realizations of the same record corresponding to the event at the respective station.However, it is acknowledged that once the approach's effectiveness in modeling structural response is confirmed, future studies can expand to encompass the modeling of scenario earthquakes for seismic assessment of masonry structures.
To summarize, future developments should consider including: (i) comparing both stochastic simulation approaches for other seismic events, (ii) adding features such as component correlations and near-fault directivity pulses to improve the site-based simulations, (iii) investigating other structural typologies to assess the reliability score for each of them, iv) exploring stochastic approaches within other numerical structural simulation tools commonly used in engineering practice, and (v) conducting simulations of scenario earthquakes and their application in the field of earthquake engineering.

D I S C L A I M E R
ANY use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

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.

F I G U R E 3 F I G U R E 4
The epicentral locations and the focal mechanisms of 26 events with Mw greater than 5.0 that occurred in the Azores region between 1939 and 2007.The yellow curves depict the active spreading centers, while the black box corresponds to the 1998 Faial event (figure adapted from67 ).Recording stations of the 1998 Faial earthquake (as shown by GZC, PVI, HOR, MOS, and SEB) along with alternative ruptured fault lines (FL1 and FL2) and their corresponding epicenters, as shown by stars.

F I G U R E 5
Investigated masonry structures: (A) U-shaped benchmarks, (B) São Francisco Church (the geometrical models are available upon request to the corresponding author).

F I G U R E 6
Real and simulated time-series, FAS and PSA at station GZC.F I G U R E 7Real and simulated time-series, FAS and PSA at station PVI.

F I G U R E 8
Real and simulated time-series, FAS and PSA at station HOR.F I G U R E 9Real and simulated time-series, FAS and PSA at station SEB.

F
I G U R E 1 5 P1 prototype: Comparison in terms of modified Hausdorff distance score (Upper +Y, Lower -Y).

F I G U R E 1 6
P2 prototype: Comparison in terms of modified Hausdorff distance score (Upper +Y, Lower -Y).

F I G U R E 1 7 8
Macro elements discretization adopted for the U-shape prototypes.F I G U R E 1Relative error-maximum displacement for each control point: (A) +X, (B) -X, (C) +Y, (D) -Y.
Shaghayegh Karimzadeh: Methodology, Conceptualization, Implementation, Collection and processing of recorded motions, Source-based stochastic simulations; Misfit evaluation methodology; Analysis and interpretation of results; Supervision on site-based simulation approach; Writing-original draft.Marco F. Funari: Methodology, Conceptualization, Implementation, FE simulations, NLRHAs simulations, Analysis and interpretation of results, Supervision, Writing-original draft, Writing-review & editing.Simon Szabó: Implementation, FE simulations, Data Analysis, Writing-original draft.S. M. Sajad Hussaini: Methodology, Implementation, Site-based stochastic simulations, Writing, Data analysis-original draft.Sanaz Rezaeian: Supervision on site-based simulation approach, Writing-review & editing.Paulo B. Lourenço: Funding acquisition.Supervision, Writing-review & editing.A C K N O W L E D G M E N T SThis work was partly financed by FCT/MCTES through national funds (PIDDAC) under the R&D Unit ISISE under reference UIDB/04029/2020, and under the Associate Laboratory Advanced Production and Intelligent Systems ARISE under reference LA/P/0112/2020.This study has been partly funded by the STAND4HERITAGE project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 833123), as an Advanced Grant.This work is also partly financed by MPP2030-FCT PhD Grants under the R&D Unit Insti-tute for Sustainability and Innovation in Structural Engineering (ISISE), under reference PRT/BD/154348/2022.This work is partly financed by national funds through FCT-Foundation for Science and Technology, under grant agreement UI/BD/153379/2022 attributed to the 4th author.This study has been partly funded by Foundation of Science and Technology, under grant agreement PRT/BD/154348/2022.
Information on the available five stations that recorded the 1998 Faial earthquake.
TA B L E 1 Input-model parameters of the stochastic source-based ground motion simulation approach.Modulation parameters of the stochastic site-based ground motion simulation approach.
TA B L E 4 Filter parameters of the stochastic site-based ground motion simulation approach.GOF scores evaluated for the real and simulated ground motion records of the 1998 Faial earthquake.
TA B L E 6