Development of parametric seismic fragility curves for historical churches

For both spiritual and cultural reasons, churches are an essential part of the historical heritage of several countries worldwide, including Europe, Americas and Australasia. The extreme damage that occurred during the 2016–2017 Central Italy seismic swarm highlighted once again the noteworthy seismic vulnerability of unreinforced masonry churches, which exhibited several collapses and caused uncountable losses to the Italian artistic heritage. The seismic performance of 158 affected buildings was analyzed in the aftermath of the main shocks. The failure modes activated by the earthquakes were identified making reference to the local mechanisms currently considered in Italy for post-seismic assessment of churches. The structural damage of the investigated buildings, related to 21 mechanisms rather than to an overall global response, was explained resorting to empirical statistical procedures taking into account ground motion intensity and structural details that can worsen or improve the seismic performance. Finally, parametric fragility curves were derived selecting those structural details that mostly influence the damage by means of the likelihood-ratio test. Developed models can be used in future territorial-scale scenario or risk analyses.


Introduction
A strong sequence of earthquakes struck Central Italy during 2016-2017, starting on August 24, 2016 causing severe damage and hundreds of casualties over a wide area (Fig. 1), within the boundaries of Marches, Umbria, Latium and Abruzzi regions (Mazzoni et al. 2018). The affected area experienced several earthquakes such as that of 1859 in Norcia (Marotta et al. 2019), that of 1979 in Valnerina and that of 1997 in Umbria-Marches , which required restorations and interventions. Besides masonry and reinforced concrete ordinary buildings, Central Italy is characterized by a 1 3 building portfolio rich of historical constructions, in fact as much as 4325 unreinforced masonry (URM) churches were inspected during the 2016-2017 seismic sequence (Penna et al. 2018). A large proportion of them delivered a poor performance, with about 45% considered unsafe and restricted for use. It is well known that churches exhibit a seismic vulnerability higher than ordinary buildings, because of their architectural and structural characteristics (D'Ayala 2000;Lagomarsino 2006). Historical URM churches respond to earthquakes with local mechanisms rather than with a global behavior. For this reason, the vulnerability assessment of URM churches is commonly based on the identification of the damage associated to possible local mechanisms. This approach was originally developed after observations following the 1976 Friuli and 1980 Irpinia earthquakes (Giuffré 1988;Doglioni et al. 1994), while after the 1997 Umbria-Marches and 2002 Molise earthquakes systematic sets of mechanisms were proposed by Lagomarsino and Podestà (2004), wherein a wider description of each of them is available. Currently, a list of 28 mechanisms is adopted in the Italian form for post-earthquake survey of churches (PCM-DPC Mibac 2006).
The large number of churches present in Italy frequently requires territorial-scale risk analyses, resorting to fragility assessment delivered from basic information collected after a rapid visual inspection, rather than performing structural analyses based on accurate geometry and details survey, as well as on material parameters investigation (Pirchio et al. 2021). Such fine verifications can be performed only at a second stage on individual buildings, shortlisted because of their higher risk.
In order to define parametric fragility curves associating the damage related to each collapse mechanism with ground motion intensity and vulnerability indicators of churches, the observed behavior of a sample of 158 Central Italy URM churches is herein analyzed. Information about these buildings is the result of inspections performed by the authors or by photo and video documentation made available by the Corps of Firefighters. Since the seismic vulnerability of URM churches is strongly modified by structural details that can worsen the seismic performance or improve it through the presence of earthquake-resistant elements (Borri et al. 2019), information about different vulnerability modifiers were carefully collected. Among these modifiers also fall strengthening interventions such as reinforced-concrete ring beams and overlays to brace the roof pitch, tie rods and so on (Sorrentino and Tocci 2008;Giuriani et al. 2009;AlShawa et al. 2019;Sferrazza Papa et al. 2021) and their role will be accounted for in the procedures described hereinafter.
The influence of each modifier on each mechanism was firstly addressed in a disaggregated fashion, through use of multiple-linear regressions and following the approach originally proposed for unreinforced masonry churches damaged by the 2010-2011 Canterbury, New Zealand, seismic sequence (Marotta et al. 2017(Marotta et al. , 2018. Whereas previous models widespread in literature explained damage by means of ground motion intensity only, the proposed regressions account also for vulnerability modifiers, delivering different damage levels to different churches subjected to the same level of shaking. Several studies were conducted over the years for the definition of churches empirical fragility curves, based on post-earthquake observations. In these works the probability of exceeding a certain level of the global damage of a single class of churches is represented as a function of macroseismic intensity (Lagomarsino and Podestà 2004;Díaz Fuentes et al. 2019) or of peak ground acceleration (PGA) (Hofer et al. 2018;Cescatti et al. 2020;Morici et al. 2020;Palazzi et al. 2020). Typological fragility curves for global damage estimation were proposed by Lourenço et al. (2013), who considered brickwork, stonework or mixed brickwork-stonework churches in New Zealand, and by de Matteis et al. (2019), who considered three-naves churches in Abruzzi region. In order to have more refined estimation of damage and fragility, the models presented hereinafter will not just consider one or few typological classes but rather several vulnerability modifiers, thus deriving tailor-made predictions for each church. Moreover, models will be derived for 21 local mechanisms, rather than for a global estimate of the church performance. The adopted approach does not require a preventive subdivision of the surveyed sample in smaller samples homogenous from a typological point of view, thus allowing to work on larger numbers and delivering comparatively more robust models. Regressions providing expected damage and parametric fragility curves delivering a probability distribution of damage will be presented.
Regressions and fragility curves will be derived in terms of PGA, because the latter is probably the intensity measure most familiar to practitioners and adopted by Italian Civil Protection Department for scenario and risk analyses (Mollaioli et al. 2018;Borzi et al. 2020;Zucconi et al. 2020). Local values of recorded PGA are attributed to each church location by means of triangulation-based linear 2-D interpolations based on shakemaps (INGV 2016a, b). Source shakemaps were computed by the Italian Institute of Geophysics and Volcanology and account for local stratigraphic conditions by means of a site classification into five classes based on geological maps of Italy (Michelini et al. 2008). This procedure is able to account for site conditions only approximately, because site stratigraphy can only crudely be estimated from small-scale maps and because topographic effects are neglected. Nonetheless, current shakemaps are a substantial progress compared to bedrock shakemaps and represent the best information currently available. Either one of the two main earthquakes of the sequence, the August 24, 2016 (M W 6.0) or the October 30, 2016 (M W 6.5) shock, is used as the reference event depending on the location of the church and the consequent date of survey. Subsequently, the damage of 49 churches out of 158 is interpreted making reference to the first event, and that of the remaining 109 churches to the end of October shock. The latter group of churches did not suffer any significant damage from the prior earthquake, because their inspection was not requested by the stakeholders until after the second main earthquake.

Characteristics of the churches
Churches have great importance for the Roman Catholic community and exhibit characteristics of unquestionable aesthetic as well as symbolic values. Unfortunately, they are also significantly vulnerable to seismic actions, since characterized by slender walls and large-span horizontal structures, often presenting thrusting elements. For the purpose of understanding the scale and nature of the seismic vulnerability of existing URM churches in Central Italy, it is useful to investigate their characteristics, both in terms of architectural traits and structural features, since damage is related to architectural types and construction details that may vary from zone to zone.
Among the sample of 158 URM churches assessed in the struck regions, 73 were located in Marches, 41 in Umbria, 32 in Latium, and 12 in Abruzzi regions. In Fig. 2 the investigated churches are sorted according to the century during which construction started, an information known for a sub-inventory of 73% of them. The majority of this sub-inventory was established in the Middle Ages (5th -14th centuries, 48%) and in the Renaissance (15th-16th centuries, 27%). The age statistic confirms the heterogeneity of the sample and that Central Italy ecclesiastic masonry-construction heritage was built over a wide time span, compared to other countries. Several churches were enlarged and/or refurbished during the centuries following their foundation, further contributing to the complexity of the investigated sample.
Following the approach previously developed for New Zealand churches (Marotta et al. 2015), a typological classification based on plan and spatial features is developed, and churches are classified into seven types, as outlined in Fig. 3a, whereas the frequency of the types among the entire stock is shown in Fig. 3b. The relative majority of churches (41%) belongs to the A type, while the B type, that includes the presence of the apse, accounts for 36% of the portfolio, highlighting that the majority of churches has a simple plan. The remaining 23% is distributed among churches with transept and apse having one (C) or Fig. 2 Occurrence of investigated URM churches according to century during which construction started. Sub-inventory of buildings for which this information is available 1 3 three naves (D), central plan churches (E), basilica churches similar to D but much larger and with three or more naves (F) and churches with two naves or with an asymmetric nave arrangement (G). Similar results can be found in Canuti et al. (2019), over a sample of 541 URM churches belonging to the Marches region alone.
Symmetry and regularity are recorded for each church. It was observed that nearly 12% of churches are asymmetrical and not regular in plan, and cases of asymmetry are often due to extensions in plan that occurred during the life of the building, or the presence of adjacent buildings. In fact, only 40% of churches are isolated and not connected to adjacent buildings. Chapels were found in 15% of the inventory, sometimes according to an asymmetric layout. A porch/nartex at the church entrance is not widespread (10%) and it is usually located in front of the main façade and opposite the church altar (75%), but sometimes is positioned on a side of the building (25%).
Domes are present only in sixteen churches of the inventory, belonging to central plan or basilica types. Bell-towers were observed in 42% of the inventory, always connected to the nave and thus possibly cause of vulnerability due to their different dynamic properties, while soaring elements, such as pinnacles and bell-gables, were identified in 46% of the sample.
Construction practice along the Apennines privileges the use of easily-available natural stone for the units of masonry (Rovero et al. 2016;AlShawa et al. 2021), although fired clay bricks are used in some areas of the Marches region . These traditions are reflected in the investigated inventory, wherein masonry units are made of stone and bricks in 82 and 17% of the churches, respectively, while in two cases stonework is used as facing for brickwork. As already widely known, the quality of masonry plays a key role in the response of URM buildings (Brando et al. 2017), especially when natural stone units have irregular shape, are small and arranged in separate leaves that do not guarantee a monolithic behavior. In the ordinary buildings of Central Italy, in the aftermath of the earthquakes, different levels of unit and mortar quality were detected and it was confirmed that the use of undressed stone units, in conjunction with low-strength mortar, often led to poor earthquake response (Liberatore et al. 2016). Dramatically low levels of lime were measured in mortars collected in the Amatrice area from damaged buildings, with the same mortar type being used both in the external leaves joints and in the internal core of the walls (Mirabile Gattia et al. 2019;Roselli et al. 2019). Churches usually present a higher masonry quality in comparison to ordinary buildings, but differences need to be recognized between main historical centers and small settlements.
The nave is covered by URM vaults in 39% of the cases, while roof battens are visible from the nave and supported by timber trusses in 20% of the inventory or by masonry structures (transversal arches or façade and end wall in small buildings) in 16% of the churches. In the remainder 25% of cases the roof was concealed by a ceiling. It is important to highlight a spread presence of "camorcanna" vault (16% of the vaulted structures), i.e. a reed mat vault made by plaster and reed laths hanging from timber centerings (Fig. 4), and thus characterized by high flexibility and low mass. Occurrences are different from those encountered in De Matteis and Zizi (2019), in which however the analyzed sample consists of 68 single-nave URM churches.
The abovementioned structural characteristics, and many others presented in the following, are relevant because they strongly modify seismic vulnerability of URM churches (Marotta et al. 2017). Accordingly, these information are used in the following for the calibration of the vulnerability models.

Damage mechanisms
The earthquake response of historical URM churches can be described by identifying separate macro-elements or macro-blocks, which are specific architectural portions (e.g., façade, transept, apse, bell tower) whose seismic behavior is only faintly concurrent with flanking parts (D'Ayala and Speranza 2003;da Porto et al. 2012;Brandonisio et al. 2013;Milani and Valente 2015;Giresini 2016). Accordingly, the damage survey and interpretation of Central Italy URM churches were based on the identification of 28 possible local mechanisms, as currently adopted in Italy for post-earthquake assessment of churches (PCM-DPC Mibac 2006). The mechanisms are listed in Table 1, graphically schematized in Fig. 5 and a selection of the most frequent of them is photographed in Fig. 6. These mechanisms are related to 10 macro-elements, namely: façade, porch, nave, transept, triumphal arch, apse, chapels, adjacent buildings, soaring elements, bell tower. For several of these macro-elements the mechanisms describe a transversal (or out-of-plane, overturning) failure mode (Andreotti et al. 2015;AlShawa et al. 2017) and a longitudinal (or inplane, shear) failure mode (Lagomarsino et al. 2019). Performance of the corresponding horizontal structures, vaults and roof, is described by additional mechanisms (Sorrentino and Tocci 2008). Damage in special building elements, such as porch, soaring elements, adjacent constructions, bell tower, is more difficult to categorize and is covered generically, without a further description of the structural response (Gizzi et al. 2014;Sorrentino et al. 2014;Degli Abbati et al. 2018). Prominent structural elements are granted supplementary mechanisms to locate more precisely the damage: gable of the façade, colonnade between the naves, lantern of the dome above the transept, belfry of the bell tower. Six levels of damage are assigned to each mechanism relying on expert judgment, ranging between 0 (no damage) and 5 (total collapse), according to the European Macroseismic Scale (Grünthal 1998).
Not all churches present all macro-elements, for instance only 16 churches have a transept with a dome and thus mechanism #14 is possible only in about 10% of the sample. Additionally, not all domes have been actually damaged, because mechanism #14 was activated (i.e., occurred) only in 11 of them, that is in about 69% of the domes. The percentages of mechanisms whose activation is possible due to the presence of the relevant macroelement are shown in Fig. 7, along with the ratio of activated-over-possible mechanisms. Some mechanisms (#9, 11) showed systematic activation (above 80%) but their macro-elements (vaults in the aisles, transept) are present in few buildings. Because of their rather poor sample size, these mechanisms, together with #10, 12, 15, 20, 24, were not evaluated further. Within the 21 remaining mechanisms, the most vulnerable is the one referred to the interactions between the nave and its roof (#19) followed by the shear response of the longitudinal walls (#6), activated respectively in 81% and 80% of possible cases. Damage in the triumphal arch (#13) and overturning interactions between the apse and its roof (#21) present the same activation rate (79%), while damage in the nave vaults (#9) and shear in the apse (#17) show a slightly lower occurrence (78%). Interactions next to irregularities (#25) and shear in the façade (#3) were observed in 77 and 76% of cases, respectively. Furthermore, all the 21 mechanisms, except the overturning of the chapels (#22), were activated in more than half of the cases, highlighting that the sample includes a large proportion of damaged buildings.
As anticipated, the seismic damage is not attributable to the intensity measure alone and vulnerability modifiers need to be accounted for (Parisi and Augenti 2013;de Matteis et al. 2016). Interactions between the nave and its roof 6 Shear in longitudinal walls 20 Interactions between the transept and its roof 7 Longitudinal response of the columns 21 Interactions between the apse and its roof 8 Damage in the vaults of the main nave 22 Overturning of the chapels 9 Damage in the vaults of the aisles 23 Shear in the chapels 10 Overturning

Regression analysis
The damage of each of the 21 mechanisms can be expressed according to a linear formulation (Marotta et al. 2017(Marotta et al. , 2018: where d is the observed damage, I the ground motion severity, x 1 , x 2 , … x v the vulnerability modifiers (Table 2), m 0 , m 1 , … m v the regression coefficients, v the number of considered vulnerability modifiers, b the intercept and ε the error term. For each mechanism, a set of fifteen modifiers is considered, related to earthquake-resistant details or vulnerability increasers. Each modifier scores between 0 and 1: an effective earthquake-resistant detail  Two statistical procedures, namely the stepwise and best subsets (Draper and Smith 1998), were used to determine the variables of the predictive regression: the stepwise selection method allows the determination of the variables that generate the most efficient predictive model, inserting variables in turn until the regression equation and each variable involve a p-value below the selected threshold of 0.05; the best subsets procedure selects the subset of parameters that optimize an objective criterion, such as having the largest adjusted coefficient of determination: where R 2 is the coefficient of determination and n is the sample size. The coefficients of the multiple-linear regressions of Eq. (1) resulting from both procedures for the considered 21 mechanisms, together with the corresponding values of R 2 adj , are reported in Table 3. The statistical checking test for excluding multicollinearity was performed for each regression model, but results are not reported for the sake of brevity. No indicator coupling was found among selected modifiers.
(2)  Table 1 Table 2 List of the vulnerability modifiers used in the multi-linear regression models

Ref. # Description
x 1 Tie rods Thrusting elements x 6 Large openings x 7 Ring beam x 8 Heterogeneous materials x 9 Connections x 10 Braced roof pitch x 11 Slenderness x 12 Asymmetry conditions x 13 Poor masonry quality x 14 Vertical-stacked-bond vaults x 15 Lunettes Table 3 Computed parameters of the selected regression models according to Eqs.  The first line for each mechanism represents the outcomes of the stepwise regression; the second line shows the outcomes of the best subsets regression Table 3 (continued)  Although not included in the current Italian form, poor masonry quality was found to be crucial for 20 mechanisms. It is also important to highlight that the presence of poor masonry quality can lead to wall disintegration (examples in Fig. 8), before a rigid-body mechanism can be activated. In the case at hand, this phenomenon was observed in 21% of the activated mechanisms.
Other very relevant modifiers are connections between intersecting walls or between walls and horizontal structures. Although it is widely known that tie rods help to reduce the overturning of the walls (Giresini et al. 2018), they seem to play a negligible role, probably due to the predominance of other modifiers such as connections. The relevance of buttresses has been investigated as well, considering them present only if made of masonry and not just of a few tens of mm of plaster, and giving them a score close to 0 if widely spaced or of small depth compared to wall thickness. Buttresses only slightly influenced the predicted damage, possibly because they are present only in about 15% of the investigated churches. Other parameters, such as large openings (whose combined length exceeds 1/3 of the wall length), heterogeneous materials (assigned in case of "camorcanna" vaults for mechanism #8 and when two adjacent structural elements are made of different masonry types), asymmetry conditions (e.g., due to eccentricity of a soaring element with respect to the underlying masonry, or due to juxtaposition of a new extension) are relevant for specific mechanisms. The removal of PGA from the set of parameters for some of the mechanisms is related to the statistical procedures that, aiming at identifying the most efficient predictive regression, treat all variables equally.
As previously stated, the two statistical procedures use a different criterion for selecting predictor variables, and thus result in different models both in terms of number of vulnerability modifiers and their coefficients. This number is smaller in the stepwise regression, despite it delivers an adjusted coefficient of determination usually close to that of the best subsets regression. Additionally, in order to maximize the adjusted coefficient of determination, the best subsets procedure may introduce variables involving a p-value above the selected threshold. Hence, the stepwise procedure is suggested for a more robust and faster territorial-scale vulnerability assessment.
The use of statistical models based on the collection of typological data implies uncertainties in the predictions. In order to assess the reliability of the proposed procedure, an investigation on the confidence intervals of the observed damage was carried out. A confidence interval is an interval estimate of the mean value computed from the statistics of the observed data, and its width provides an idea of the uncertainty about its estimation. It has an associated confidence level that quantifies the level of confidence that the mean value will lie in such interval. The confidence level is usually chosen equal to 0.05, meaning that there is a 95% probability that the linear regression line of the population will lie within the confidence interval computed from the sample data (Ross 2014). Accordingly, the lower the confidence level specified, the larger the estimated range that is likely to contain the line. Regression lines of the observed damage for two of the 21 considered mechanisms are shown in Fig. 9, together with confidence intervals and, for comparison reasons, prediction bounds. Prediction bounds provide information on individual predictions of the accounted dependent variable, giving a range of values where an additional observation of the variable can be expected to be located. In fact, prediction intervals provide a range of values where one can expect future observations to fall and are useful when the aim is using the model to predict individual values of the response. As shown in Fig. 9, the confidence interval is much smaller than the prediction interval, because the former is an interval estimate for an average value while the latter is an interval estimate for a single observation. This result confirms the recommendation of using territorial scale analyses prominently for average estimation rather than individual forecasts.

Parametric fragility curves
Availability of site-specific hazard curves has made risk analysis an attractive complement to standard structural analysis related to just one or two damage states. Similar considerations apply also to territorial-scale scenarios based on statistical approach and, thus, it is of great interest to compute not just damage regression, as in the previous section, but a complete set of fragility functions.
The fragility function F i (PGA/g) for a certain damage level i (i = 1,…,5) is defined as the probability to reach or exceed damage level i, given the PGA (Porter et al. 2007): The probabilities to reach a certain damage level i, given PGA/g, are: The typical expression adopted for fragility functions is the lognormal distribution: In order to derive empirical fragility curves, starting from the observation of earthquake-induced damage and based on the vulnerability modifiers of Table 2, a linear law providing µ i for each damage level was used: where the term m 0i depends on damage level, and the coefficients m 1 , m 2 , …, m ν account for the effect of the modifiers. It is worth to notice that the latter are assumed to be independent from damage level. In other words, the modifiers operate a simple "shift" of the set of fragility curves, as functions of ln(PGA/g), or a scaling of the curves, as functions of PGA/g. An example is shown in Fig. 14 in the Appendix. Although in principle this hypothesis can be easily removed introducing coefficients depending on damage level, the sample size does not allow for a robust estimation of such coefficients. Standard deviation, besides independent from damage level, is in addition assumed to be independent from the modifiers.
The parameter set θ contains the coefficients m 0i (i = 1,…, 5), the coefficients m 1 , m 2 ,…, m ν and the standard deviation β. To appraise the parameters, the maximum likelihood estimation (Benjamin and Cornell 1970) was used, whereby the parameters of the model are obtained by maximizing the likelihood function L(θ), so that under the assumed statistical model the observed data is most probable. Fig. 9 Regression lines of the observed damage with scatter plot of the data with confidence (in red) and prediction (in dashed blue) intervals for: a mechanism #1; b mechanism #17 In order to identify the set of vulnerability modifiers that mostly influence the damage, as performed with the regression models of Sect. 4, the likelihood-ratio (LR) test was carried out, by removing from the initial statistical model one at a time the modifiers representing the structural characteristics of the churches. The null hypothesis, i.e. the modifier can be removed, is then checked according to the LR test (Neyman and Pearson 1933). In fact, the LR test assesses the goodness of fit of two competing statistical models based on the ratio of their likelihoods, specifically one found by maximization over the entire parameter space Θ and another constraining the parameters in a subspace Θ 0 of Θ. In the case at hand, Θ contains all the coefficients m 1 , m 2 , …, m ν related to the vulnerability modifiers, whereas Θ 0 contains a reduced set. Both Θ and Θ 0 enclose the coefficients m 0i (i = 1,…, 5) and β. If the constraint (i.e., the null hypothesis) is supported by the observed data, the two likelihoods should not differ by more than a threshold. The LR test statistic can be expressed as: where, in the case at hand, sup ∈Θ 0 L( ) represents the likelihood maximized over the subspace Θ 0 with some parameters removed, and sup ∈Θ L( ) the likelihood maximized over the complete initial space Θ. The test statistic can also be expressed as a difference between the log-likelihoods: where If the null hypothesis is true, Wilks' theorem ensures that λ converges asymptotically, for large samples, to a chi-square (χ²) distribution with k 0 degrees of freedom, equal to the difference in dimensionality between Θ and Θ 0 (Wilks 1938).
When a small difference occurs between the likelihood maximized over the complete space Θ and the likelihood maximized over the subspace Θ 0 , the value assumed by λ is small, its right-tail probability α is large, and the null hypothesis is accepted, i.e. the parameters can be removed from the model. An α-value greater than or equal to 0.05 is commonly interpreted as justification for accepting the null hypothesis, i.e. there is no significant difference between the complete and the reduced model, thus allowing the removal of the parameters from the model. For each mechanism, the set of significant parameters was obtained by removing one at a time each parameter, computing the corresponding value of α, and selecting the parameter to be removed as the one associated to the reduced set with the largest value of α. The procedure was repeated until a set of parameters delivering a value of α < 0.05 was reached (Table 4). The mean likelihoods, i.e. the geometric means of the likelihoods over the sample, are reported in Table 5 for the reduced model L m0 (with significant modifiers) and for the complete model L m (with all the modifiers): where n is the sample size. The test statistic λ, the degrees of freedom k 0 and the α-values are reported as well. In addition, the corresponding quantities for the model without modifiers, denoted as • , are reported. Appreciable values of likelihood are obtained for most mechanisms, with lower values for those less populated. For almost all the mechanisms, at least two modifiers are significant. Moreover, values of the mean likelihood of the reduced model, L m0 , are in most cases close to the mean likelihood of the complete model, L m , whereas the mean likelihoods of the model without modifiers, − L m0 , are sensibly lower. Thus, the introduction of modifiers entails a significant improvement of the models and it is confirmed that only a reduced set of modifiers is required. Exceptions occur for mechanisms #4, #7 and #23, where the procedure delivers a final model depending only on the ground motion intensity measure and disregarding any other vulnerability modifier.
Once again, poor masonry quality was found to be crucial for most mechanisms, followed by the presence of connections. It is also worth noting, for mechanism #8, that the high value associated with the presence of heterogeneous materials assigned in case of "camorcanna" vaults is due to the fact that this type of vault was identified only in case of damage or collapse, consequently the presence of the modifiers was always associated with high levels of damage. As for mechanism #1, the negative value associated to the presence of buttresses is counterintuitive and, as already mentioned in Sect. 4, can be explained by the significant spacing usually observed between buttresses in façades and their small depth compared to the wall thickness, leading to a 0.33 score in all façades where they were present. Such façades were just 11 out of 132, so additional data may lead to a future revision of this result.
The proposed procedure delivers tailor-made fragility curves for each church, based on the vulnerability modifiers actually present and on their score. Therefore, the fivedamage-levels fragility curves of each possible mechanism need to be computed on a case-by-case approach, as shown in a worked-out example in the Appendix. Sample fragility curves related to eight mechanisms are presented in Fig. 10. These mechanisms are among the most populated (Table 5) and related to several macro-elements, namely: façade, nave (including longitudinal walls and colonnade), triumphal arch, soaring elements, bell tower. The curves of another two mechanisms are shown in the Appendix. The scores assumed for the relevant modifiers are the average values of the churches in the investigated portfolio and are reported in Table 6 for all 21 mechanisms, not just the eight presented in Fig. 10. In the same table the standard deviation of the scores and the frequency of each modifier, scoring more than zero, for each mechanism are reported as well. Therefore, combining Tables 4 and 6, it is possible to compute the fragility curves for each mechanism, assuming an "average" church. The curves in Fig. 10 differ from mechanism to mechanism under several aspects. First of all, some mechanisms (e.g., #13, 21) are more fragile than others (e.g., #2), encouraging a mechanism rather than a global approach as privileged so far in the literature. Additionally, mechanisms #2, 28 present curves closer one to the other compared to #6, 13, emphasizing that for the former once light damage is activated, heavy damage follows more rapidly. Despite all curves being log-normal by definition, some mechanisms (e.g., #7, 21, 28) present an Table 4 Computed coefficients of the parametric fragility curves  evident "S" shape, while others present a very short initial horizontal tail, noticeable if a logarithmic scale is adopted for the horizontal axis. Such shape is customary for large β values compared to µ values. It is possible that a future, more populated database will deliver updated curves with a more pronounced "S" shape, without affecting the feasibility of the method proposed here. Additionally, for each mechanism, the damage observed for each church in the aftermath of the earthquakes was compared with the predicted damage obtained with the stepwise and best subsets regression models according to Eq. (1), and with the mean damage, D m , derived from the fragility curves: In all mechanisms a good correlation and a coefficient of determination R 2 equal to or greater than 0.5 are found. On the left of Fig. 11 comparisons related to four mechanisms, among the most populated and corresponding to different macro-elements (namely: façade, nave, apse, bell tower) are shown. Performance of regression models and fragility curves in reproducing observed damage is similar. In fact, the comparisons between predicted damage of regression models and mean damage of fragility curves show a noteworthy correlation, with values of R 2 in all cases larger than 0.8 and in most cases larger than 0.9 (right of Fig. 11).
Finally, it is worth highlighting that regressions provide a predicted damage, whereas the proposed parametric fragility curves deliver a distribution of probability for each damage level, resulting in an enriched statistical model.

Conclusions
The 2016-2017 Central Italy earthquakes strongly affected the architectural and historical national heritage, among which unreinforced masonry churches suffered heavy damage. The extensive data collected in the aftermath of the two main shocks for a sample of 158 religious buildings highlighted the intrinsic structural fragility of this architectural type. In literature, the fragility of churches generally was analyzed evaluating global damage, as a function of ground motion intensity, for just a class of buildings or for very few types. In this study, damage was interpreted mechanism by mechanism, accounting for differences in vulnerability. The damage suffered by the churches was correlated to shakemap PGAs by means of 21 local mechanisms and was initially expressed as a multiple-linear regression of intensity measure and modifiers worsening or improving the seismic response. The regression models were evaluated according to the stepwise and best subsets procedures, and the corresponding coefficients were provided. In addition to the statistical regression models, fragility curves were developed, whose parameters are identified by maximum likelihood estimation. The introduction of modifiers resulted in an increase of the likelihood compared to models relying on ground motion intensity alone. The main benefit of the proposed model was that each mechanism of each church had its own fragility curves, "tailored" on the basis of the vulnerability modifiers. The most influential modifiers, among which were masonry quality and connections, were identified according to the likelihood-ratio test. A comparison of the used procedures was conducted, showing good correlation with the observed damage and between the different models. Developed models can be used in future territorial-scale scenario or risk analyses, resorting to basic information acquired after a rapid visual inspection. The procedure described in previous sections can be used again to update the fragility curves if refined ground motion shakemaps and updated observation databases become available.

Fig. 11
Left column: comparison of observed damage with predicted damage from regression models and mean damage from fragility curves; Right column: comparison between predicted damage from regression models and mean damage from fragility curves. SW = stepwise regression; BS = best subsets regression; FC = fragility curves Table 5 Computed statistics of the fragility curves for complete models, reduced models and models without modifiers: n     Table 6 (continued) x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9 x 10 x 11 x 12 x 13 x 14  Table 6 (continued) x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9 x 10 x 11 x 12 x 13 x 14

Appendix: calculation section
To better explain the procedure presented above, a worked-out example on one of the considered Central Italy URM church is reported hereinafter. The building under consideration is a stone masonry church (Fig. 12), whose occurred PGA during the August 24, 2016 shock was 0.30 g. Only eight mechanisms are analyzed, because some macro-elements, specifically those related to mechanisms #4, 7-8, 14, 16-18, 21-23, 25, are not present. Mechanisms #27-28 are not considered because at the time of the survey the bell tower had suffered a complete collapse and it was not possible to retrieve information about its characteristics. In Table 7 an overview of the modifier attribution is reported, varying between 0 and 1 depending on the presence and effectiveness of vulnerability increasers and earthquakeresistant elements. For the sake of conciseness, only some modifiers are discussed in detail. Since the façade does not present buttresses whatsoever, a zero score is attributed to the corresponding modifiers of mechanism #1. Likewise, the absence of tie rods implies a zero score for the same mechanism, whereas the presence of an undersized ring beam leads to a 0.67 score assigned to mechanism #3. A unity score is introduced to account for the poor masonry quality, because the church is made of undressed natural stone units with an external regular facing and an inner core of rubble masonry.
The second part of Table 7 shows the parameters of the fragility curves: the values of μ i associated with the different damage levels computed according to Eq. (6), the probabilities to reach each damage level i given PGA/g according to Eq. (4a-c), the mean damage D m derived from fragility curves according to Eq. (11) and the damage d observed after the August 24, 2016 shock. The good agreement between D m and d can be noticed. In Fig. 13, fragility curves of two mechanisms of the sample church are shown. An example of the effect produced on the set of fragility curves by the modifier x 13 , related to the poor masonry quality, is shown in Fig. 14.  Table 7 Modifier attribution and computed parameters of the fragility curves for each mechanism of the sample church   licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.