A new method for earthquake-induced damage identification in historic masonry towers combining OMA and IDA

This paper presents a novel method for rapidly addressing the earthquake-induced damage identification task in historic masonry towers. The proposed method, termed DORI, combines operational modal analysis (OMA), FE modeling, rapid surrogate modeling (SM) and non-linear Incremental dynamic analysis (IDA). While OMA-based Structural Health Monitoring methods using statistical pattern recognition are known to allow the detection of small structural damages due to earthquakes, even far-field ones of moderate intensity, the combination of SM and IDA-based methods for damage localization and quantification is here proposed. The monumental bell tower of the Basilica of San Pietro located in Perugia, Italy, is considered for the validation of the method. While being continuously monitored since 2014, the bell tower experienced the main shocks of the 2016 Central Italy seismic sequence and the on-site vibration-based monitoring system detected changes in global dynamic behavior after the earthquakes. In the paper, experimental vibration data (continuous and seismic records), FE models and surrogate models of the structure are used for post-earthquake damage localization and quantification exploiting an ideal subdivision of the structure into meaningful macroelements. Results of linear and non-linear numerical modeling (SM and IDA, respectively) are successfully combined to this aim and the continuous exchange of information between the physical reality (monitoring data) and the virtual models (FE models and surrogate models) effectively enforces the Digital Twin paradigm. The earthquake-induced damage identified by both data-driven and model-based strategies is finally confirmed by in-situ visual inspections.


Introduction
Several methods have been validated in the literature for earthquake-induced damage identification, aimed at the preventive conservation of historic masonry buildings. Among them, vibration-based structural health monitoring (SHM) methods resorting to modal parameters as damage-sensitive features and based on unsupervised learning techniques are eminently efficient for early-stage damage detection (Bennati et al. 2005;Ivorra and Pallarés 2006;Alvandi and Cremona 2006;Brencich and Sabia 2008;Ramos et al. 2010;Spina et al. 2011;Ditommaso et al. 2012;Brando et al. 2015;Gentile et al. 2016;Dolce et al. 2017;Masciotta et al. 2017;Clementi et al. 2017;Saisi et al. 2018;Kita et al. 2019;Formisano et al. 2020;García-Macías et al. 2020a;Giordano et al. 2020a, b;Bianconi et al. 2020;Standoli et al. 2021). Unfortunately, their main drawback is that they are limited to detection only, even though (Saisi et al. 2015;Ubertini et al. 2018;Venanzi et al. 2020) providing a first clue on damage quantification, e.g. in the form of permanent variations in natural frequencies.
Methods addressing further levels of damage identification need to be developed and validated in the literature for cultural heritage (CH) preservation. The introduction of the concept of Digital Twin, whose first applications concerned mechanical engineering for life prediction (Grieves and Vickers 2017), has recently gained increasing scientific interest in civil engineering in the case of historic masonry constructions. According to (Angjeliu et al. 2020), the Digital Twin is composed of three main items: (i) physical model and experimental reality, (ii) collected data characterizing the experimental reality (e.g. vibration-based monitoring data) and (iii) the virtual reality or virtual simulation model. In other words, it can be defined as a model or a suite of models of the structure that continuously exchanges information with the physical reality through experimental/field data, thus providing evidence of novelties or changes, for instance, related to a faulty or damage condition.
While damage detection can be considered typically a data-driven process, damage localization may require the inverse calibration of a Finite Element (FE) model of the structure. The main goal is to minimize the mismatch between the numerical and experimental responses (typically natural frequencies and mode shapes) by the calibration of uncertain damage-related parameters (e.g. material properties, or boundary conditions) of the model. In this context, changes in the modal features can be directly related to damage-induced variations in the mechanical parameters of the structure (Atamturktur and Laman 2012; Sehgal and Kumar 2016;Pellegrini et al. 2018). The complex geometry and spatial distribution of different masonry materials of most historic buildings coupled with a large number of simulations needed for solving the minimization problem determine an unsustainable computational burden posing major limitations in the real-time SHM practice. Surrogate models (SM) have been recently approached as rapid and computationally efficient tools to reproduce the behavior of a time-demanding numerical model, enabling damage identification through continuous model updating compatible with continuous SHM systems (Moustapha and Sudret 2019;Liu et al. 2020;Pepi et al. 2020;Micheli et al. 2020a, b). Some experiences in the literature can be mentioned about the application of surrogate-based model updating of historic structures (Cabboi et al. 2017;Torres et al. 2017;García-Macías et al. 2021Venanzi et al. 2020). The damage localization task can be performed by solving an inverse calibration problem on a surrogate representation of a linear FE model and by using longterm vibration data. In formulating the inverse problem, equivalent elastic properties of macrostructural elements can be considered as uncertain parameters to be identified by continuously minimizing an objective function. The latter can exploit experimentally identified and numerically predicted (damage-induced) decays in natural frequencies as well as variations in mode shape components. Noteworthy to say, the abovementioned works on SM report applications to masonry towers considering artificially simulated damages, whereby validations with real earthquake-induced damage are still missing in the literature.
On the other hand, the recent growth of computer processing power and structural analysis software packages has enabled a constant evolution towards increasingly accurate but at the same time more complex analysis methods for the seismic assessment of masonry constructions, from (elastic) linear static (LSA) to linear dynamic analysis (LDA), and from non-linear static (NLSA) to non-linear dynamic analysis (NLDA). The latter can provide a more accurate and reliable assessment of the structural response of masonry buildings (Zampieri et al. 2017;Manzini et al. 2018;Abbati et al. 2019). It consists of running several non-linear time-history analyses considering earthquake records as inputs (each record produces a ''single-point'' analysis) and can be performed either by using a high number of seismic records (Cloud Method and/or multiple stripe analysis (MSA)) (Lagomarsino and Cattari 2015) or considering a proper (still significant) selection of time histories to be scaled for performing incremental dynamic analysis (IDA) (Vamvatsikos and Cornell 2002). The proper definition of the non-linear behavior of the masonry material is required for NLDA, allowing a rigorous comprehension of the seismic behavior of the structure along with its possible failure mechanisms. To the actual state-of-the-art, the application of IDA to CH masonry structures is comparatively less common with respect to other types of structures, such as RC frames, mainly due to the higher computational costs (Basone et al. 2017;Masaeli et al. 2018). Nevertheless, a few applications are worthful to be mentioned in (Pena et al. 2010;Karanikoloudis and Lourenço, 2018;Kita et al. 2020). The research carried out in  proposed the use of an IDA-based approach for damage localization and quantification in the case of a laboratory low-rise masonry structure. Although effective on such a scaled-down model, the application of the method in full-scale masonry structures that have experienced real earthquake-induced damage has not yet been explored. This paper proposes a new method for damage detection, localization and quantification in long-term monitored historic masonry towers subjected to earthquakes. It firstly relies on the data-driven OMA-based approach for a rapid and automated damage detection. Subsequently, the extension of the OMA procedure is proposed by further implementing the paradigm of Digital Twin, with the main goal to increase the reliability of the identified damage. In the context of the present paper, the Digital Twins are models of the considered structure that are updated over time based on on-site continuous monitoring data and models that yield a prediction when an earthquake is recorded by the permanent vibrationbased monitoring system. More specifically, the combination of two independent innovative methods is proposed for damage localization and quantification: (i) SM-based damage identification and (ii) IDA-based damage identification. With the former relying on linear numerical modeling and long-term vibration monitoring data, and the latter on a non-linear modeling and seismic on-site response to earthquakes, their integration and cross-validation represent the pivot feature of this work.
The proposed method is applied and validated on the San Pietro bell tower. It is an iconic slender masonry construction located in Perugia, Italy, that has been continuously monitored by the authors since December 2014. The results demonstrate that the proposed method enables a rapid post-earthquake damage identification, allowing to immediately detect the presence of damage and to subsequently localize and quantify it with an acceptable level of confidence in different macrostructural elements.
The paper is organized as follows. Section 2 illustrates the proposed method for rapid earthquake-induced damage identification in long-term monitored historic masonry towers. Section 3 presents the case study, while Sect. 4 reports the results and validation of the proposed method. Finally, Sect. 5 summarizes the main conclusions of the work.

Description
The proposed method consists of a data-driven vibration-based approach (using OMA) combined with two independent innovative model-based methods, aimed at the detection, localization and quantification of earthquake-induced damages in historic masonry towers under continuous long-term monitoring. The method is named DORI and it addresses the Damage identification based on the use of Operational modal analysis, Rapid surrogate modeling and Incremental dynamic analysis. The description is provided in the next lines (steps I, II, III, IV) and the general framework is summarized in the flow chart illustrated in Fig. 1.
Preliminary activities for deep knowledge acquisition are: (i) geometrical survey, (ii) material characterization, (iii) on-site inspection, (iv) structural condition assessment, (v) damage survey related to possible crack pattern analysis and material degradation mapping (Lourenço, 2006;Ubertini et al. 2016;Caprili et al. 2017;DeMatteis et al. 2019). They also provide the basis for ambient vibration tests (AVT) that allow evaluating the dynamic characteristics of the building. The overall gained information is essential for the implementation of the permanent vibration-based SHM system on the one hand, whose continuous data are necessary for long-term OMA, and the construction and calibration of Digital Twins on the other.
Step I: Data-driven OMA-based damage detection It comprises four consecutive stages: (i) continuous measurement of the dynamic response of the tower through a permanent monitoring system, (ii) continuous automated modal identification and frequency tracking, (iii) multiple data regression for the removal of changing environmental and operational conditions from identified modal frequencies and (iv) novelty analysis for damage detection, i.e. statistical analysis of damage-sensitive features to detect any novelty or anomaly in the investigated system (Farrar and Worden 2012). In order to define damage-sensitive features that are minimally affected by changing environmental conditions, tracked modal frequencies contained in an observation matrix, ∈ ℝ n×N , and the same modal frequencies independently estimated using a statistical model and stored in a matrix ̂ are used to construct the so-called residual error matrix, ∈ ℝ n×N , which is defined as follows: Methods for constructing ̂ in Eq. (1) are Multiple Linear Regression (MLR) and Principal Components Analysis (PCA), among others. Both methods work well when there is a linear correlation between natural frequencies and environmental parameters (e.g. temperature), which is often an acceptable simplification. While MLR statistically constructs a linear model relating a set of independent variables (predictors) to (1) = −̂ the natural frequencies (estimators) (Dervilis et al. 2015), PCA is a model dimension reduction technique that works by (i) projecting the estimators into the space of their principal components, (ii) retaining only a few principal components explaining most of the variance in the estimators and (iii) going back to the initial space of the estimators through the use of such a reduced number of principal components (Yan et al. 2005). Under the hypothesis that the retained principal components represent the effects of the unobserved predictors, PCA can yield similar results when compared to MLR. Hybrid combination of MLR and PCA may also be proposed when only some of the most significant predictors are measured .
The damage condition is identified as an anomaly in the residual errors contained in , e.g. a change or deviation in their distribution. Statistical process control tools (novelty analysis in stage iv) are adopted for this purpose, using control charts constructed by properly defined statistical distances, such as Q 2 and T 2 statistics (Worden et al. 2000Sohn et al. 2002;Fuller 2009;Mosavi et al. 2012;Farrar and Worden  where r is the group averaging size parameter, ̄ contains the mean values of the residuals in the last r observations subgroup, while ̄ and are the mean values and the covariance matrix of the residual errors, respectively, statistically estimated in the training period. Subsequently, in the observation period, an anomaly in the data is identified in the form of outliers (values of T 2 greater than the UCL) steadily observed over time after a seismic event.
Step II: SM-based damage identification The primary scope of a SM is to bypass in a cost-effective way the input/output relationship of a computationally demanding model. Its construction can be described by the following stages: (i) selection of design variables, (ii) definition of the design space, (iii) generation of the training population and (iv) construction of the SM (García-Macías et al. 2020b). The design variables, denoted as x i ∈ ℝ, i = 1, … , m (x = [x 1 , …, x m ] T ), define the design space and allow the parametrization of the FE model and investigation of potential damage scenarios. They represent damage-sensitive parameters and are defined in a way that meaningfully vary within a certain range. Thus, the vector of m design variables defines the m-dimensional design space ℕ = x ∈ ℝ ∶ a i ≤ x i ≤ b i , with a i and b i representing the range limits. Afterwards, the generation of a training population of N individuals requires drawing input samples uniformly over ℕ , considering the m x N matrix of design sites X = [x 1 , …, x N ] T ∈ ℝ m×N . In this way, Z = [z 1 , …, z N ] T is defined as an observation vector, where z i ∈ ℝ represent the system's response to the input x i . In this context, a SM can be regarded as a model-free mathematical representation of the response z(x) of a computationally demanding FEM, with x being the vector of design parameters and modal properties assumed as outputs. This representation indicates that, as long as the SM is optimally calibrated at first, it can directly relate damage-sensitive mechanical properties (e.g. stiffness constants) to modal properties (e.g. natural frequencies and mode shape components) without the need for running any FEM simulation.
The SM-based method proposed in this paper, previously developed and validated in the case of similar masonry structures (Venanzi et al. 2020;García-Macías et al. 2020b), relies on the use of long-term monitoring vibration data and numerical model updating. The realtime damage detection and localization are performed by solving an inverse FE model calibration problem using a SM. To this aim, the calibrated linear FE model of the structure is ideally subdivided or partitioned into meaningful macrostructural elements (referred to as macroelements hereafter) on the basis of expected damage and engineering judgement. In this work, the scalar multipliers (k i , i = 1, …, m) of Young's moduli (E i ) and shear moduli (G i ) of all materials contained in m macroelements ( Fig. 1) are defined as damage-sensitive design variables. It is noteworthy to stress that each material contained in the i-th macroelement is associated to the same i-th stiffness multiplier, so that the procedure aims to identify a sort of average damage within a macroelement without attempting at its precise definition. According to the adopted partition, E i and G i of all the elements contained in the single macroelement are thus defined as random variables as follows: (2) where E i 0 and G i 0 are the undamaged elastic moduli of a material within the i-th macroelement, respectively, while parameters k i represent the design variables x = [k 1 , …, k m ] T in the SM-based damage identification approach. The inverse calibration problem continuously identifies these parameters by minimizing a certain objective function, defined in the following lines. Distinct permanent reductions (e.g. induced by a seismic event) in the components of x would reveal the presence of damage condition in the corresponding macroelements.
Among several models available in the literature, the response surface method (RSM) showed to be simple, yet effective for SM representation (García-Macías et al. 2020b). The second-order quadratic formulation of the RSM is adopted in this work, written as follows: where 0 , j , jj and ji are coefficients that require an initial calibration, while the term θ refers to the error between the original FEM and the SM.
The objective function is defined as follows: where relative differences between experimental and numerical modal parameters are represented through the terms ε i and δ i for natural frequencies and mode shapes, respectively, while α and β are the corresponding weighting coefficients. The term η has a regularization purpose against possible ill-conditioning problems (García-Macías et al. 2020b). The minimization of the objective function yields the time histories of identified damage-dependent parameters (from hereafter stiffness multipliers k i ) for each macroelement. The SM-based damage localization and quantification procedure proposed in the present work consist of: (i) online SM calibration, i.e. continuous calibration of SM frequencies on the basis of experimentally identified natural frequencies, highlighting variations in stiffness multipliers k i due to random errors in modal identification, errors in SM calibration, residual effects of changing temperature and operational conditions that the statistical model does not effectively depurate and variations due to damage; (ii) automated detection of damageinduced shifts in the time histories of k i through control charts by applying Eqs. (1) and (2) to k i within each macroelement; (iii) quantification of earthquake-induced decays for each macroelement through a comparison between post-earthquake and pre-earthquake values of k i .
Step III: IDA-based damage identification The IDA-based method was previously developed and validated by some of the authors in the case of an international benchmark laboratory masonry structure, called Brick House . It consists of multidimensional seismic IDA simulations carried out using a non-linear FE model and seismic data recorded via an on-site vibration-based SHM system. After numerical NLDA simulations, sets of IDA curves are constructed with reference to distinct structural macroelements, allowing multidimensional relations between a set of meaningful local damage-dependent parameters or Damage Measures (DMs) and properly selected earthquake ground motion Intensity Measures (IMs). Unlike the Brick House, where only seismic input IMs have been investigated, the IDA-based procedure is here proposed for a long-term monitored bell tower, whereby the use of seismic response IMs is also considered, using recorded seismic events.
Local multidimensional IDA curve sets are proposed in this work using a reasonable DM (e.g. tensile damage parameter, d t ) and a suite of uncorrelated IMs according to the following equation: where m again indicates the number of macroelements in which the structure is ideally subdivided for damage identification purposes and n * the total number of IMs. It is noteworthy to stress that the IDA-based procedure uses the same ideal partitioning/subdivision into macroelements defined in the SM-based approach for damage identification. Given the high computational cost required by IDA, the curve sets can be constructed during the training period of monitoring, during which there are presumably no seismic events. At the occurrence of an earthquake, the seismic/response intensity parameters are measured or directly computed from the measurements of the monitoring system of the structure, and local damage conditions in each i-th macroelements ( DM * i,j , with i = 1, …, m and j = 1, …, n * ) are immediately estimated using the previously constructed multidimensional IDAs. Local damage conditions are finally averaged in terms of weighted mean ranges and weighted mean values as follows: where RMSE i,j are root mean square error coefficients computed on the basis of the IDA curve sets dispersion, and their inverse values are considered as weights of the n * most efficient IMs.
Step IV: Cross-validation of results from previous step This is the pivot step of the DORI method where the cross-validation of results obtained from steps I, II and III is applied and, in particular, multiple consecutive agreement checks are investigated for consistency in the different stages of earthquakeinduced damage identification. First, data-driven damage detection is double-checked through damage detection achieved with SM and IDA-based approaches. When the agreement is ascertained at this stage, it is then progressively investigated for damage localization and quantification tasks. More specifically, the OMA-based and the SMbased approaches assess a potential increase of the relative frequency of the outliers to be observed in the control chart after a seismic event: (i) the former indicating the detection of permanent changes in global structural behavior, while (ii) the latter localizing damage in the macroelement which exhibits the highest values of such outliers. Also, in the short-term approach, right after the occurrence of a seismic event, an immediate earthquake-induced damage identification consists of: (i) evaluating instant decays in the time histories of stiffness multipliers obtained by the SM-based approach and (ii) estimating damage by the IDA-based approach (through IMs calculated from on-site recorded seismic responses) in terms of reduction of the adopted DM. Finally, in the long-term approach, the relative frequency of the outliers from the SM-based procedure computed over a significant amount of time after the seismic event can be reasonably related to the earthquake-induced damages quantified by means of the IDA-based method. However, in case of no (or weak) agreement between the results at each stage, additional preliminary investigations would be necessary for more refined knowledge acquisition and for improving the SM and the non-linear FEM model.
Beyond the concept of the cross-validation of results, it is crucial to understand that techniques of steps I, II and III, are complementary in terms of modeling strategies and achieved damage identification levels. Therefore, while their separate implementation is certainly useful, their combination is the key to enrich the accuracy of the proposed DORI method. It is also worth noting that, once the statistical models are constructed (OMA), once linear (SM) and non-linear (IDA) FEM models are built and calibrated, once the SM is validated and once the IDA curves are generated, which is all done during the training period (typically during the first year of monitoring), the DORI method is fully automated and operates in real time with a very limited number of sensors. Given the minimal hardware cost required, a similar approach is believed to be highly convenient when dealing with towers of high historical value which justifies the initial modeling effort.

Innovative aspects
The DORI method integrates several innovative aspects. The main one is the introduction and combination of Digital Twins in the context of long-term OMA-based monitoring for localization and quantification of earthquake-induced damage. Two independent innovative model-based methods are jointly investigated, relying on the use of SM and IDA, respectively, in order to reduce uncertainties in damage identification. Overall, the integration of different and complementary methods aims at the most accurate earthquake-induced damage identification (detection, localization and quantification). OMA alone has been successfully applied in the literature for damage detection. The SM is already implemented in a few examples in the literature but needs to be further investigated in the case of historic masonry bell towers, while IDA is newly proposed, not being observed so far in the literature for damage localization and quantification in historic masonry structures. On the one hand, the proposed SM-based method uses a linear model and combines long-term vibration monitoring data and numerical modeling, as a rapid and computationally efficient tool for damage identification. On the other hand, unlike the SM, IDA requires a non-linear model and applies the on-site recorded seismic response to earthquakes. OMA, SM and IDA-based methods are all calibrated during the training period of the monitoring system, while subsequently applied practically in real time after the occurrence of a seismic event. Several benefits can be exploited by comparatively investigating the outputs of the three approaches (cross-validation): (i) first of all, damage detected by continuous OMA can be double-checked with damage identified with SM and IDA-based approaches, allowing to highly mitigate the risk of false-positive results; (ii) there is an integration between linear (SM) and non-linear (IDA) modeling, as well as between long-term vibration monitoring data (OMA and SM) and the seismic response to earthquakes (IDA). This integration can be exploited for comparative and/or complementary damage identification purposes, potentially reducing the uncertainties related to IDA and enriching the accuracy of damage localization and quantification; (iii) the SM-based approach can quantify damage in terms of earthquake-induced decays of damage-dependent parameters, while the IDA-based method provides a quantitative damage assessment through DMs. To the best of the Authors' knowledge, no other method is available in the literature to identify earthquake-induced damage in masonry towers with such a limited number of vibration sensors as proposed in this paper.

Description of the structure
The monumental bell tower of the Basilica of San Pietro is considered to be a major symbol of the city of Perugia due to its historical and architectural value. Dating back to the thirteenth century, it was subjected to several structural and architectural interventions throughout history, where the last one was carried out in 2002 after damages caused by the 1997 Umbria-Marche earthquake (Cattaneo et al. 2000). Damages were concentrated predominantly in the belfry (severe cracks at the base of the columns and diffused cracking pattern on the external walls) and the cusp. Today, it is in a good state of preservation and no significant and visible damages were observed after the major earthquakes occurred in Italy in recent years: L'Aquila earthquake in 2009, Emilia earthquake in 2012 and Central Italy seismic sequence in 2016 .
With a total height of 61.4 m and restrained up to the first 17 m by the surrounding Basilica, main cloister and abbey, the bell tower can be considered composed of three main structural parts, as illustrated in Fig. 2: (i) the shaft with a dodecagonal cross-section, (ii) the belfry with a hexagonal cross-section and (iii) the pyramid-shape cusp completing the tower on the top. A low-cost permanent vibration-based SHM system was installed on top of the tower in December 2014 with the purpose of monitoring its structural integrity. It currently comprises three high-sensitivity uni-axial accelerometers (model PCB 393B12 with 10 V/g sensitivity) placed at the base of the cusp (see Fig. 2a), two temperature sensors (K-type thermocouples) and a data acquisition system. The use of three uni-axial accelerometers at the top floor of the structure (A1, A2 and A3) demonstrated over time to be sufficient for monitoring the global dynamic behavior of the building, and more specifically, for the tracking of the flexural modes and the first order torsional mode, while minimizing the cost of measurement hardware, which is crucial for a broader application of the proposed monitoring solutions. The monitoring system is remotely connected to a dedicated server in the Laboratory of Structural Dynamics of the Department of Civil and Environmental Engineering, University of Perugia. Continuous monitoring data have been processed through a MatLab code based on a fully automated SSI output-only modal identification technique (Ubertini et al. 2013). The adopted OMA technique includes the following main steps: (i) continuous vibration data acquisition, (ii) run of SSI (CVA formulation) data analysis, considering the model's order varying from 40 to 60 (with step increments of 2) and the number of output block rows of Hankel matrix ranging from 140 to 200 (with step increments of 10), (iii) noise modes elimination through reasonable relative rejection thresholds between modal parameter estimates (0.01, 0.03 and 0.01 for relative changes in frequencies, in modal damping ratios and differences between MAC values, respectively), (iv) clustering analysis of remaining modes and, finally, (v) extraction of mean values of modal parameter estimates with 95% confidence intervals. Investigation of seasonal and daily fluctuations of identified modal parameters, associated with changes in environmental conditions (primarily in ambient temperature), can be found in (Ubertini et al. 2017;Giordano et al. 2020b).

The FE model
A 3D numerical model of the San Pietro bell tower has been built within the framework of the Finite Element Method (FEM) (see Fig. 2c), in the ABAQUS platform (Simulia 2010). The structure is discretized with hexahedral C3D8R elements in the shaft, while tetrahedral C3D4 elements have been adopted in the belfry and the cusp. The FEM model is composed of 300,722 elements, 134,868 nodes and 340,989 degrees of freedom. At the constitutive level, stones and brick masonry have been modeled as a single homogeneous material, by using both isotropic and orthotropic formulations for different purposes, as specified in the following sections. The interested reader can find the tuned mechanical parameters of the constituent materials in , together with more details on the FE model, boundary conditions and the calibration procedure.
The non-linear FE model has been constructed wherein the mechanical behavior of the masonry material has been reproduced using the concrete damage plasticity (CDP) constitutive model (Lubliner et al. 1989;Lee and Fenves 1989), which is broadly acknowledged in the literature for seismic vulnerability and NLDA (Formisano et al. 2018;Valente et al. 2019;Milani and Clementi 2019;Acito et al. 2021). As available in ABAQUS, the use of CDP requires the adoption of an isotropic constitutive behavior for the FEM material, because, from a theoretical point of view, it is based on isotropic elasticity in combination with isotropic plasticity. In the present case, the behavior of the masonry material has been characterized considering damage in tension only. The strength parameters assigned to the FE model have been earlier estimated in . They have been defined in 1 3 terms of tension stiffening and tension damage, whereby the failure condition and the postpeak behavior is described by quantities such the tensile stresses, σ t , the cracking strains, ∼ ck t , and the tensile damage variable, d t .

Seismic input for IDA
Seismic input choice/selection is a key issue in NLDA of structures (Vamvatsikos and Cornell 2002;Kita et al. 2020). Among different types of accelerograms, natural ones are the best to be used. Seven groups of natural seismic records (two horizontal components each) have been obtained by the software REXEL (Iervolino et al. 2009) to be used as input for IDA of the San Pietro bell tower. Considering the local seismic hazard of the structure, the seismic input has been defined in terms of acceleration time series which are spectrumcompatible, i.e. whose response spectra result compatible with the site-specific response spectrum of the building. The unscaled spectrum-compatible accelerograms are plotted in Fig. 3, where the main original characteristics of the seven selected ground motions are also reported. They have been progressively scaled at increasing levels using appropriate scaling factors (SF). The seismic loading is bidirectional, whereby components have been applied in the two horizontal directions to the FE model of the tower: x component in E-W and y in N-S direction, respectively.

OMA-based earthquake-induced damage detection
The OMA-based procedure (step I of the proposed method) has been applied before and after the San Pietro bell tower was hit by the three main shocks of the 2016 Central Italy seismic sequence: (i) Accumoli Mw6.0 earthquake of August 24th and the following shocks, (ii) Ussita Mw5.9 earthquake of October 26th and (iii) Norcia Mw6.5 of October 30th. The bell tower is located at a distance of about 85 km in the NW direction from the epicenter of the first earthquake and approximately 70 km from the epicenters of the other two. Five years of monitoring data have been used to this aim (December 9th 2014-December 9th 2019), while one year of training period has been considered. Time series of identified natural frequencies (collected in ) are illustrated in Fig. 4a.
Fx1 indicates the first flexural mode in the x-direction, while Fy1, Fy2 and Fy3 are the flexural modes in the y-direction (x and y directions are illustrated in Fig. 2a, referring to East-West and North-South, respectively). Finally, T1 refers to the torsional mode. The effects of changes in environmental and operational conditions have been removed by applying a statistical model based on a combination of MLR and PCA to the monitoring data (estimated modal frequencies are collected in ̂ ). Afterwards, the statistical model for novelty analysis is built according to Eqs. (1) and (2). Despite no relevant structural damages were observed in any part of the structure (above all in the belfry), the longterm vibration-based procedure detected earthquake-induced damage: (i) small consistent decays in natural frequencies have been identified (see zoom plots in Fig. 4a) and subsequently (ii) an anomaly has been automatically revealed as a result of the control chart. The latter is illustrated in Fig. 4b, where small permanent changes of the structural behavior of the bell tower (anomalous deviations from normal conditions) can be highlighted after the Accumoli earthquake (sharp increase in the relative frequency of outliers), and also later after Ussita and Norcia shocks. Another indication of anomalies is obtained by further analyzing the control chart. Probability density functions (PDFs) of the T 2 statistical distance computed before and after the seismic sequence are displayed in Fig. 4c. Higher average value and increased dispersion of T 2 can be observed after the seismic sequence (the histogram is relocated towards the right as flatter and longer). With an UCL equal to 95%, the relative frequency of outliers changes from 5.44% before the earthquake sequence to 43.69% after the sequence.

SM construction and validation
For the construction of the SM, the orthotropic calibrated linear FEM has been ideally subdivided into three distinctive macroelements M i , i = 1, …, 3, and then parametrized through a set of damage-sensitive features according to Eq. (3). Macroelements M 1 , M 2 and M 3 refer to the shaft, the belfry and the cusp, respectively, as shown in Fig. 2, whereby the stiffness multipliers k i constitute the design variables x = [k 1 , k 2 , k 3 ] T of the San Pietro bell tower in the SM-based damage identification approach. Thus, the SM has been constructed on the basis of a training population, whose design space has been formed by k 1 , k 2 and k 3 , whereby permanent reductions in one of these components would conceivably indicate the presence of damage in the corresponding macroelement. Random variables have been uniformly generated by using a Latin Hypercube Sampling (LHS) method for a homogeneous representation. LHS is an effective statistical method representing the generalization of the concept of the Latin square (only one sample) to an arbitrary number of dimensions (multidimensional distribution). In the present surrogate modeling approach, LHS has allowed the generation of the training population in terms of random uniform variables-hereby stiffness multipliers k i -to accurately calibrate the surrogate model on the basis of the linear FEM model. The coefficients k i have been defined with upper and lower bounds such that the design space assumes the form of ℕ = x ∈ ℝ 3 ∶ −0.25 ≤ x i ≤ 0.15 . Generating a training population of N = 2048 individuals, Fig. 5 depicts the PDFs of the natural frequencies obtained from the 2048 analyses on the FEM, while Fig. 6 illustrates them plotted versus those predicted by the SM. The very low scatter around the diagonal suggests that the SM is constructed with high accuracy. The frequencies of the original FE model are also reported in both figures to highlight the distance between the undamaged and the damaged FEM.

Damage identification
Subsequently to OMA, earthquake-induced damage has been identified by applying the SM-based approach before and after the 2016 Central Italy seismic sequence (step II of the proposed method). Time series of cleansed natural frequencies and experimental eigenvector components have been used as input for solving the inverse FE model calibration problem on the RSM-type SM defined by Eqs. (4) and (5). Cleansed natural frequencies have been obtained by summing the vector of average values of natural frequencies tracked during the training period to the vector containing statistical model residuals E (see Eq. 1).
The damage localization and quantification consist of investigating the output of the SM-based optimization procedure, i.e. tracking the variation in time of the coefficients k i that multiply the elastic and shear moduli of materials belonging to M 1 , M 2 and M 3 , respectively (Eq. 3). The optimization procedure has allowed investigating different combinations of weighting coefficients α and β, as reported in Table 1, with the purpose of weighting distinctly the two terms of the objective function in Eq. (5), search for stable solutions avoiding unfeasible results and aim for accurate damage quantification. Figure 7 depicts earthquake-induced relative decays of k i obtained by applying the investigated combinations of α, β and η into the SM-based procedure. Relative decays have been quantitatively evaluated for each macroelement by comparing values identified in the last two days before each shock with those identified in the first two days after. The k 2 multiplier experiences the highest variation (of the order of 10%), denoting that damage has occurred in the belfry macroelement, while the shaft and the cusp have significantly lower values. Also, the cusp presents slightly higher decays if compared with the shaft (of the order of 3-5% and 1-3%, respectively). It is worth noting that damage localization in the belfry is consistent with the damage scenario predicted in ) through non-linear FE simulations, as well as with the damage pattern that was observed on the structure after stronger earthquakes occurred in the past (cfr. Sect. 3.1). Another observation is that Accumoli and Ussita earthquakes have induced similar decays in k i , a consistent result if compared  . 7 Earthquake-induced relative decays of k i for each seismic event obtained by applying different combinations of α, β and η into Eq. (5) of the SM-based procedure with OMA-based earthquake-induced relative decays of natural frequencies. In combinations C1, C2 and C3 (cases when α < β) some solutions that deviate from the others are observed. According to the adopted objective function in Eq. (5), the inclusion of the term δ i (related to mode shapes) in addition to the term ε i (related to resonant frequencies) has allowed investigating the robustness of the solution with respect to the weighting coefficients variation. In this context, Fig. 7 illustrates that the contribution of the mode shapes' term to the solution accuracy is quite limited because only three eigenvector components are available in the present case. As better discussed later on in the manuscript, the best combination in Table 1 seems to be A1 (α = 1 and β = 1). Figure 8a depicts the time histories of the stiffness multipliers k i obtained for such a combination. Zoom plots illustrate small consistent earthquake-induced decays that have been identified for the belfry (k 2 ) in correspondence to the seismic events, instant decays that are almost absent for the shaft and the cusp. The three k i vectors (shaft, belfry and cusp) have been then used to build the statistical models for novelty analysis. In particular, distinct control charts have been built by separately using the k i time series into Eq. (2). They are illustrated in Fig. 8b, where a clear deviation from normal conditions in the case of the belfry macroelement can be highlighted after the main earthquakes. In the case of the shaft and the cusp, the T 2 statistical distance keeps the vast majority of values below the UCL (95%).
The localization of earthquake-induced damage in the belfry is confirmed by further analyzing the distribution of control charts. PDFs of T 2 before and after the seismic sequence are displayed in Fig. 8c for each macroelement. It is noted that the histogram of the belfry is more shifted in comparison with that of the shaft and the cusp. The relative frequency of outliers in the case of the belfry changes from 2.21% before the earthquake sequence to 36.41% after the sequence, while those of the shaft and the cusp remain under 2.5% (lower than UCL). Finally, to emphasize the fact that combination A1 is the best for damage detection and its localization and quantification in the belfry macroelement, control charts obtained from combinations A2 (α = 5 and β = 1), B2 (α = 5 and β = 0) and C2 (α = 1 and β = 5) are illustrated in Fig. 9. It can be observed that all combinations consistently localize the damage in the belfry (larger frequency of outliers), but damage detection is generally less effective in A2, B2 and C2 compared to A1.

Non-linear seismic IDA curve sets
The choice of IMs and DMs for IDA purposes is a very important issue (Vamvatsikos and Cornell 2002). Nineteen weakly correlated and meaningful IMs have been investigated in the present paper, as obtained from a statistical correlation analysis regardless of the structure under consideration (Peak Ground Acceleration PGA, Root Mean Square Acceleration RMSA, Characteristic Intensity I C , Arias Intensity I A , Destructive Potential Factor P D , Cumulative Absolute Velocity CAV, Spectral acceleration S a (T1), Acceleration Spectrum Intensity ASI VT , Peak Ground Velocity PGV, Root Mean Square Velocity RMSV, Specific Energy Density SED, Cumulative Absolute Displacement CAD, Spectral velocity S v (T1), Housner Intensity I H , Peak Ground Displacement PGD, Root Mean Square Displacement RMSD, Spectral displacement S d (T1), Displacement Spectrum Intensity DSI NH and Input Energy E I ), while tensile damage d t has been used as DM .
The IDA curve sets of the San Pietro bell tower have related d t to the aforementioned seismic input IMs and allowed investigating local damage conditions for the three distinct macroelements M i , i = 1, …, 3: shaft, belfry and cusp (see Fig. 2). Despite the seismic loadings for IDA are bidirectional in the present case, with two horizontal components applied to the FE model (see Fig. 3), the IDA curve sets are graphically represented by means of IMs computed using their mean direction. The DM is a volumeaveraged damage parameter, i.e. it has been computed as an average value weighted over the volume of the finite elements of every single macroelement.
The curves presenting a relatively high dispersion have been firstly discarded, and only the less scattered have been considered for IDA-based damage identification. By operating in this way, the most efficient seismic input IMs have been identified as those contained in vector IM* = (PGA, RMSA, I C , I A , S a (T1), ASI VT , PGV, RMSV, S v (T1), I H , S d (T1)). Figure 10 shows the plots of curves obtained with some significant IMs: acceleration-related IMs (PGA, I A and ASI) and velocity-related parameters (PGV and I H ). Sets of seven curves and their mean curve are graphically illustrated in each graph for the shaft, the belfry and the cusp. By a rapid visual investigation, it is possible to immediately and directly identify the belfry as the most vulnerable macroelement, whereby a clear increase in the slopes of the curves can be highlighted and slight-to-moderate seismic events can conceivably produce higher damage compared with the two other macroelements. The maximum numerically obtained tensile damage d t is equal to about 0.85 in the belfry, while greater IMs' levels are needed to observe similar important damage in the cusp. The shaft seems to be the least damaged with a damage parameter value of about 0.2 (corresponding to a 20% reduction of Young's modulus). These first interpretations of the IDA curve sets are consistent with numerical damaging results presented in past works in the literature conducted by some of the authors Ubertini et al. 2018).
The dispersion of the IDA curve sets has been investigated using the root mean square errors RMSE i,j of single curves computed with respect to the mean curve (Eq. 7). Figure 10 reports the coefficients of dispersion, RMSE i,j , for the shaft, the belfry and the cusp. In general, a greater dispersion can be observed in the case of the cusp, meaning that the adopted DM is not robust enough in relation to IMs of this macroelement. However, the main focus is concentrated on the belfry, whose dispersion is quite limited: the most efficient IM in terms of consistently lower dispersion is I A .

3
Seismic response IMs have been also investigated for IDA purposes. They have been computed from numerical seismic response acceleration, velocity and displacement time histories from a pivot node of the FEM during IDA. It is the node corresponding to accelerometers A1 and A2 of the monitoring system illustrated in Fig. 2. In this context, the IDA curve sets relate mean seismic response IMs (computed from East-West and North-South components) to local volume-averaged damage parameters. The most efficient response parameters are the following seven: Peak Response Acceleration (PRA), . PRA, PRV and PRD are peak parameters, while RMSRA, RI C , RI A and RMSRV represent integral measures. It is worth noting that the definition of RI C and RI A represents a novelty in the literature. Some IDA curves are depicted in Fig. 11 (PRA, RI A and PRV). A rapid investigation indicates that there is a clear consistency with the observations reported from the IDA curve sets obtained with seismic input IMs: the belfry is confirmed as the most vulnerable macroelement, while the cusp and the shaft present significantly smaller IDA-based damages. Also, the greatest dispersion is again observed in the case of the cusp. Figure 11 also reports the RMSEs coefficients for each IDA curve set, subsequently needed for obtaining weighted mean damage values through Eq. (7).

Damage identification
The same main earthquakes of the 2016 Central Italy seismic sequence have been used for IDA-based damage identification in the case of the San Pietro bell tower (step III of the proposed method). The ground motion records were appropriately downscaled from the nearest accelerometric station , while the availability of experimentally on-site recorded seismic responses is very precious. They have been used to calculate seismic input and seismic response IMs, respectively, which combined have allowed Fig. 11 The IDA curve sets and corresponding mean curves: plots of tensile damage, d t , versus seismic response IMs (PRA (a), RI A (b) and PRV (c)). Plots of the last column report the dispersion (RMSEs coefficients) for each IDA curve set further refinement of IDA-based damage quantification. Plots of the ground motion waveforms and the time histories of seismic accelerations recorded on top of the bell tower (East-West and North-South components) are illustrated in Fig. 12.
Once all the IDA curve sets have been constructed and their dispersion investigated (see Figs. 10 and 11), local damage conditions have been assessed by using the mathematical relations of Eq. (6) between DMs and IMs of the selected real seismic events. Eleven seismic input IMs have been computed from the three base seismic records and collected in vectors IM *A , IM *U and IM *N , while seven seismic response IMs have been calculated from the three top seismic records and collected in vectors RIM *A , RIM *U and RIM *N . Figure 12 illustrates damages estimated in the belfry macroelement by applying Eq. (6) Fig. 12 Belfry macroelement: IDA-based tensile damage estimated by using eighteen IMs (calculated from Accumoli, Ussita and Norcia earthquakes) into Eq. (6). FEM damage (d t ) is also reported for comparative purposes separately using all eighteen IMs as computed from Accumoli, Ussita and Norcia earthquakes, respectively. IDA-based damages are reported in terms of minimum, maximum and mean values (obtained from the mean curves), with a range estimation zone (between the minimum and maximum values) highlighted by a light grey color. Damage predicted through a direct non-linear FEM simulation ) is also reported for comparative purposes. It consists of numerically computed mean tensile damage values (for the shaft, belfry and cusp) obtained from the FEM non-linear dynamic analyses after applying accelerometric inputs as seismic loadings. A good consistency between the latter and IDAbased predictions can be observed, whereby d t values fall within the corresponding ranges. However, both aspects of consistency between FEM damages and mean IDA predictions and the desirable small ranges of IDA damage are of particular importance in IDA-based damage quantification, therefore, their compromise is reasonably needed. To this aim, the idea of not using all IMs, but carefully selecting only some of them (proper combinations), can improve the results obtained from the IDA procedure . The best combinations here exploited are reported in Table 2.
Earthquake-induced damages obtained from every single IM as illustrated in Fig. 12 have been averaged according to Eq. (7) through weight coefficients which depend on the IDA curve sets dispersion (computed as the inverse of RMSE i,j ). Figure 13 reports IDA-based tensile damages estimated in terms of weighted average ranges (delimited by  Table 2): Accumoli (a), Ussita (b) and Norcia (c). Direct comparison in the case of the belfry (d) maximum and minimum values) and weighted mean values for the shaft, belfry and cusp. A good agreement between FEM damages and IDA-based estimated damages can be highlighted. Also, a careful selection of only some IMs is possible due to the wide possibility given by eighteen IMs, aiming at improving the accuracy of IDA-based damage quantification. Three proper combinations have been investigated: combination 1 comprises all eighteen IMs (n * = 18), while the 2 nd and the 3 rd combine only some parameters. Figure 13 also illustrates colored contour plots of damage pattern that have been numerically obtained at the last time-step of IDA (contour range from 0-no damage to 1-full damage), where significant concentration of tensile damage can be highlighted on the belfry. Finally, Fig. 13d illustrates weighted mean IDA-based tensile damages plotted versus FEM damages, d t , in the case of the belfry. Switching from combination 1 to 3, there are better matching results with a closer allocation along the diagonal. It is worth noting that the best consistency is obtained when seismic input and seismic response IMs are combined.

Cross-validation of results: OMA + SM + IDA
This section is dedicated to the application of step IV of the proposed DORI method. It consists of critically and comparatively investigating the results hierarchically obtained from steps I, II and III. Figure 14 progressively summarizes each step of damage identification after the main events of the 2016 Central Italy sequence in the case of the San Pietro bell tower. The earthquake-induced damage has been detected by the OMA-based procedure first, and then ascertained from the SM and IDA-based approaches in terms of consistent decays of damage parameters and reduction of volume-averaged tensile damage Fig. 14 Step IV. Hierarchical illustration of results of earthquake-induced damage identification on the San Pietro bell tower: OMA, SM and IDA of the three macroelements, respectively. Afterwards, the agreement checks on earthquakeinduced damage localization have been carried out. The SM-based procedure has provided time histories of stiffness multipliers k i , where substantial decay in the 2nd damage parameter and the corresponding control chart have localized damage in the belfry macroelement. On the other side, the IDA-based approach has confirmed damage localization predominantly in the belfry by using all seismic responses recorded both at the base and the top of the tower. The last task concerns earthquake-induced damage quantification. The SM-based procedure has provided the entity of decays of stiffness multiplier k 2 of the belfry macroelement equal to about 10% for each seismic event, as output directly depending on the combination of α and β weights of the objective function. In order to overcome the uncertainties of the SM-based approach, also related to the use of a linear model, the damage has been complementarily estimated by the IDA-based method through IMs calculated from seismic responses recorded during the earthquakes. It has been quantitatively assessed equal to about a 30% tensile damage parameter for each shock. The IDA-based approach can be expectedly viewed as a more conservative method for damage quantification because it directly uses the recorded seismic responses, while the SM-based approach relies on longterm data that may recover gradual changes over time. Also, a brief discussion about the adopted DM for IDA purposes may reveal useful in the case of the San Pietro bell tower, whose non-linear mechanical behavior of the FEM masonry material is described by the CDP model that is based on isotropic damaged elasticity in combination with isotropic tensile plasticity. If compressive plasticity had also been considered, the DM of each macroelement would have been a sort of mean value between tensile damage (d t ) and compressive damage (d c ), because these parameters are directly computed from tensile and compressive equivalent plastic strains, respectively (Simulia 2010). In this way, the earthquake-induced reduction of volume-averaged DM could be reasonably considered lower, for instance about half of what reported in Fig. 13, i.e. equal to about 15%. In this context, it can be stated that a relatively good agreement on damage quantification is achieved between the SM and IDA-based approaches.
Overall, all agreement checks of step IV of the DORI method have been satisfied in the case of the San Pietro bell tower. As a final result's investigation, Fig. 15 illustrates another, more intelligible perspective on the cross-validation of the results obtained from steps I, II and III. Figure 15a compares the OMA-based with the SM-based results. The former Fig. 15 Step IV. Cross-validation of the results of earthquake-induced damage identification on the San Pietro bell tower: OMA vs. SM (a) and SM vs. IDA (b) highlighted a clear increase of the relative frequency of the outliers, from 5.44 to 43.69%, indicating the detection of changes in the global structural behavior after the earthquakes of the 2016 Central Italy seismic sequence. The latter localized damage in the belfry where the relative frequency of outliers goes from 2.21% before the earthquake sequence to notably 36.41% after the sequence (instead, in the case of the shaft and cusp very low values are observed, below 2.5%). Figure 15b depicts earthquake-induced damages quantified by the IDA-based method in terms of mean values from each seismic event that are plotted against the SM-based damages in terms of relative frequency of the outliers after the seismic sequence. Insignificant damages can be observed for the shaft and the cusp, while notable damage is observed for the belfry, as consistently quantified through both approaches.
To conclude, visual inspections have been carried out by the authors, revealing a limited, yet clearly visible damage state. Figure 16 illustrates in-situ photo evidence of the Fig. 16 Damage pattern developed in the belfry macroelement after the 2016 Central Italy seismic sequence: tensile damage's contour plot obtained with non-linear seismic IDA finite element simulation and corresponding in-situ photo evidence of the main cracks main cracks developed in the belfry macroelement after the 2016 Central Italy seismic sequence: (i) in correspondence of key sections of the pointed arches and (ii) at the base portions of the columns and in the bottom masonry walls (comprised between the columns) of the belfry, some of which are also characterized by openings. This damage pattern appears consistent with the non-linear seismic IDA finite element simulation.

Conclusions
The major objective of this paper has concerned the development and the validation of the proposed DORI method addressing earthquake-induced damage detection, localization and quantification in historic masonry towers. The extension of the OMA-based data-driven damage detection method through the introduction and implementation of Digital Twins, in the present case two independent computational models (SM and IDA) exchanging information with the physical reality, represents the most innovative aspect of the method. These models provide significant contributions to reduce uncertainties in damage identification, given their complementary nature in terms of input data (i.e. long-term stationary data in the case of SM and short-term seismic data in the case of IDA) and modeling approaches (i.e. linear and non-linear FEM for SM and IDA, respectively). Overall, the integration of the three methods is conceived to allow an accurate earthquake-induced damage identification (detection, localization and quantification) with a minimal hardware requirement.
The results of the proposed method have been presented in the case of the monumental San Pietro bell tower, an historic masonry structure located in Perugia, Italy, that has been continuously monitored for five years using only three accelerometers and two temperature sensors and that was hit by important seismic events during the 2016 Central Italy seismic sequence. The main results of such a field validation study are summarized below.
(i) While most of the published works on SM applied to masonry towers relied on artificially simulated damages, a meaningful field validation with real earthquake-induced damage has been achieved in this paper for the first time. The control charts of the stiffness multipliers k i of the three selected macroelements (shaft, belfry and cusp) outputted by real-time SM calibration allowed to clearly localize the damage induced by the 2016 Central Italy seismic sequence in the belfry, highlighting a consistent and permanent decay in k 2 which represents a first metric of damage severity. (ii) The use of IDA for SHM was early proposed by the authors in a recent paper, where the investigated structural typology was a small brick structure tested in the laboratory and validation of the proposed method was fully based on numerical simulations. In this work, the IDA-based method has been applied for the first time to a full-scale masonry tower undergoing real earthquake-induced damage. In this context, seismic response IMs (computed from the top measurements of the tower), have been shown equally effective with respect to more classical seismic input ones (computed from the base recordings), which is relevant for minimizing hardware cost. Results have confirmed damage localization in the belfry, while also allowing damage quantification in terms of tensile damage. (iii) The capability of SM alone to enable accurate damage predictions is questionable as one may argue that a linear model with local stiffness multipliers is not able to catch the complex damaging mechanisms of masonry and that the inverse SM calibration may be ill-conditioned or may get stuck in a local minimum. On the other hand, IDA analysis may also result in uncertain predictions due to epistemic uncertainties related to the difficulty in calibrating a non-linear mechanical model and due to the high variability in earthquake input records. The combination of SM and IDA, i.e. the DORI method, has been therefore newly proposed in this work, with the main goal to increase the reliability of the identified damage and reduce the probability of a missing or a false alarm. In the present case study, the combination of the two methods has allowed confirming damage localization in the belfry and, more importantly, to relate the shift and the relative frequency of outliers observed in stiffness multiplier k 2 by application of the SM to the weighted average values of tensile damage parameter, thus achieving a more accurate damage quantification. It is also worth noting that in-situ visual inspections have confirmed the existence of a slight damage in the belfry, consistent with the DORI method predictions.
Based on the results presented above, it is concluded that the DORI method is a valuable tool for earthquake-induced damage detection, localization and quantification in masonry towers, whose initial modeling effort, which should be carried out during the first year training period, is believed to be justified in application to towers of high historical value, given the fully automated and real-time implementation of the method afterwards (during the observation period) and its minimal hardware requirements.