The assessment and risk management of masonry churches and, more generally, of all cultural heritage buildings require a structured and ad hoc approach. National and international guidelines (ICOMOS 1964, 1999) pursue this intent by providing methodologies and best practices to safeguard cultural heritage while ensuring the safety of surrounding communities. Many of the lessons learned by experience over the years were conveyed, for example, into the Management guidelines for World Cultural Heritage sites − 1998 (Feilden and Jokilehto 1998), the ISCARSAH Principles from the International Scientific Committee for the Analysis and Restoration of Structures of Architectural Heritage – 2003 (ICOMOS/ISCARSAH Committee 2003) together with the Annex on Heritage Structures of ISO/FDIS 13822 (ISO 13822:2010), the Risk Preparedness: A Management Manual for World Cultural Heritage – 2007 (Stovel 2007), the Operational Guidelines for the Implementation of the World Heritage Convention – 2024 (UNESCO WHC 2024), and the ISCARSAH Guidelines for the Analysis, Conservation and Structural Restoration of Architectural Heritage − 2024 (ISCARSAH 2024). These documents have highlighted critical perspectives on the preservation and risk management of cultural heritage, as well as the profound and far-reaching impact of relatively rare catastrophic events on sites of cultural heritage. Moreover, they advocate a paradigm shift in conservation efforts, moving from reactive approaches to a more holistic, proactive framework centred on prevention. This preventive approach not only ensures better preparedness for potential disasters but also enhances the long-term resilience of cultural heritage sites by addressing vulnerabilities before they manifest into irreversible damage. However, limited funding, lack of technical expertise, and competing priorities often hinder the effective application of these frameworks. Acknowledging these challenges, this paper proposes a new seismic vulnerability index to quantify the physical vulnerability of historical masonry churches to seismic events. Despite the already existing seismic vulnerability indices, which analyse churches by decomposing them into macro-elements, each associated with one or more potential collapse mechanisms, the proposed index relates the building's seismic performance solely to geometrical characteristics. Similar proposals can be identified, for example, in (Lourenço and Roque 2006; Palazzi et al. 2019), both of which introduced simplified indices based on a few geometric and material indicators, that are intended for the immediate screening of a large number of buildings, followed by more detailed investigation. However, the latter, as with most vulnerability index proposals available, parameterises vulnerability solely through resistance indicators, without accounting for seismic event intensity.

On the contrary, in the present study, physical vulnerability is defined as the expected degree of loss to vulnerable elements resulting from interaction with the hazardous event for which the hazard is estimated (ISSMGE TC32 2004). The aim is to provide a simplified, geometry-based index suitable for comparative analysis across a large sample of similarly structured churches.

Quantitatively, vulnerability is calculated by the aggregate effect of the intrinsic physical capacity of vulnerable elements to withstand an event, called resistance (R) and predominantly determined by physical features, and the damaging potential of an event, parameterised by its intensity (I). This aggregate effect is represented through fragility curves derived from multiple Nonlinear Static Analyses (NSA) of representative archetypes (Del Carlo et al. 2025), where material variability is accounted for by treating masonry mechanical properties as random variables. Indeed, as these properties are difficult to ascertain in most real cases, incorporating them directly within the nonlinear analyses allows their influence to be captured directly within the fragility curves.

The several fragility curves, each obtained through the multiple NSA performed for each archetype, are then aggregated using the Generalised Linear Model (GLM) technique to derive a quantitative seismic vulnerability index. The index represents the expected tendency of vulnerability across a properly selected range of seismic intensities, scaled from 0 to 1, where 0 implies a zero probability of exceeding a given damage level and one the opposite.

Overall, the proposed framework provides a geometry-based, simplified yet effective vulnerability estimate, for preliminary regional-scale vulnerability analysis. This could be particularly useful in contexts where only limited information is available, enabling vulnerability classification based on a small set of basic and measurable geometric indicators. These indicators can often be acquired through tools such as GIS or existing databases, avoiding costly and time-consuming on-site survey campaigns. Moreover, using a 0–1 scaled index, conceptually grounded in the relationship between structural resistance and hazard intensity, allows comparison of the obtained seismic index with vulnerability indices for other types of events, thereby ensuring broader applicability across diverse scenarios.

Historical churches are especially vulnerable to seismic events due to their typical masonry construction, lack of seismic-resistant features, and peculiar architectural configurations. The structural response of these buildings under seismic loading is further complicated by their architectural and historical significance, which restricts invasive retrofitting. Because of these factors, the study of historical churches requires dedicated, simplified methodologies to analyse their physical vulnerability at a large scale. The methodology here introduced combines simplified structural modelling, probabilistic capacity and demand analysis, with data-driven vulnerability modelling. The proposed workflow consists of the following three main steps (Fig. 1): capacity analysis, derivation of fragility points – defined as single probability values of exceeding a selected damage state for specific PGA levels – and regression-based vulnerability modelling. Each step is detailed in the following sections, together with the symbology adopted.

Capacity analysis

Capacity analysis aims to investigate the structural behaviour of buildings under external loads. The collapse or near-collapse response of masonry structures is usually analysed through incremental-iterative analysis procedures. These analyses apply external loads in successive increments, with each step producing a corresponding displacement increment. Incremental-iterative analyses can be categorised into non-linear static (pushover) and non-linear dynamic approaches. Although non-linear dynamic analysis is usually regarded as the most accurate method for capturing the seismic response of masonry structures, the application of time-history analysis is limited in practice by its high computational cost, the need for precise constitutive models to describe the non-linear hysteretic behaviour of masonry, and the challenges in selecting and scaling suitable ground motion records (Lagomarsino and Cattari 2015). Pushover analysis, despite being based on the static approach, remains a widely used and valuable tool for analysing historical masonry buildings. Its relative simplicity, reduced computational demand, and ability to distinguish sequential activation of failure mechanisms make it particularly suitable for preliminary evaluations, and it is therefore considered the standard tool for the seismic Performance-Based Assessment (PBA) of existing buildings (D’Ayala 2014; D’Ayala et al. 2014; Guardiola-Víllora, Molina, and D’Ayala 2023; Kocaman and Kazaz 2023).

Pushover analysis is performed on a set () of representative archetype models. Church archetypes refers to representative structures that reflect the variability in geometric and material characteristics within a church typology (Mata, Kalagasidis, and Johnsson 2014). Following the approach proposed by Rota, Penna, and Magenes (2010), the geometric configuration of each archetype is treated deterministically, with fixed values assigned to the plan and elevation layouts and the wall thickness. In this way, the overall geometric variability of the church typology investigated is preserved, even though each archetype is defined deterministically. In contrast, material variability is explicitly accounted for within each archetype by treating masonry's mechanical properties as random variables. The set of mechanical properties realisations () can be sampled from predefined probability distributions derived from national code references, literature, or experimental data. For each archetype, a set of (with ) pushover analysis is conducted, each corresponding to a different realisation of the material properties, Fig. 2a. The resulting pushover curves describe the relationship between base shear (or equivalent spectral acceleration) and lateral displacement at a selected control point on the structure, allowing tracking of the structural performance up to collapse or significant damage (Fig. 2b). Displacement capacity points can be identified along the capacity curve based on selected performance parameters (e.g., base shear or displacement). Since curves are generated for each archetype, they collectively define a distribution of displacement capacity points associated with the selected damage state. This distribution can be further processed, for example, using a Kernel Density Estimation (KDE) (Węglarczyk 2018) method to obtain smooth empirical Probability Density Functions (PDFs), as shown in Fig. 2c. This process must be repeated for each archetype to describe the overall geometric variability of the church typology.

Fragility points derivation

Analytical fragility curves enable the establishment of a direct mechanical relationship between the building’s construction features, its structural response to seismic loads, and the probability of failure or the probability of reaching a given damage state. In this framework, this probability is computed using Monte Carlo Simulation (MCS). First, the displacement demand imposed on the archetype by different levels of seismic intensity is obtained using the N2 method (Fajfar 2000; FEMA 2003). Assuming each archetype behaves as a Single-Degree-Of-Freedom (SDOF) system, capacity curves can be directly bi-linearised using the energy equivalence criterion, preserving the area under the original capacity curve. These idealised elastic–perfectly plastic curves are then intersected with a proper set () of Acceleration Displacement Response Spectra (ADRS), as shown in Figs. 3a and 3b, corresponding to increasing Peak Ground Acceleration (PGA) levels. Each intersection produces a displacement demand value for a given material realisation and PGA level. By repeating this process across all the bi-linearised capacity curves (each representing a unique material configuration), a distribution of displacement demand is obtained for each PGA (Fig. 3c). As for capacity, KDE may be used to derive a smooth probability density function of the demand values.

Once the capacity and demand distributions are obtained, a MCS is performed, sampling n values from each distribution, where n can determined by the smallest probability to be estimated with sufficient accuracy. Figure 4a provides an example of the capacity and demand distributions. Notably, the simulation must be repeated for all PGA values in the initial set. An indicator function () is then defined to count the number of times demand () exceeds capacity (), with . The indicator function assigns a value of 1 when demand exceeds capacity and zero otherwise. By counting the number of times the indicator function equals one and dividing by the number of samplings (n), the frequentist estimate of the probability of failure is obtained (Eq. 1).

where if demand exceeds capacity, and = 0 otherwise. The result represents the probability of reaching or exceeding the selected damage state, for each i-th PGA value in the set, i.e. a fragility point conditioned to the PGA. This procedure is repeated across all PGA values, yielding a discrete set of fragility points. These are then fitted using a lognormal CDF, resulting in smooth analytical fragility curves (Fig. 4b). Eq. 1 implies that the PGA values causing the failure, or exceeding the selected damage state, are defined among the positive real numbers, consistent with the physical meaning of intensity measure itself, which cannot assume negative values and is lognormally distributed (Baker 2015; Eads et al. 2013). Moreover, the lognormal CDF can define the relationship between PGA and the probability of damage with only two es, namely the mean and the standard deviation, leading to a comparatively simple regression process widely used in the literature.

Regression-based vulnerability modelling

Following the derivation of fragility curves for each archetype, the next step is to formulate a vulnerability index that represents a church's seismic performance based on its geometric indicators and seismic intensity level. The primary objective is to move from probabilistic fragility functions toward a more operational and interpretable tool for comparing the vulnerability of churches and guiding prioritisation at a regional scale. To this end, a regression-based approach is adopted. Given the need for interpretability in regional-scale analysis, a linear model is considered, enabling the definition of a vulnerability index as a linear combination of weighted geometric indicators (). Successful applications of linear regression to develop vulnerability indices within the context of regional-scale analysis can be found in studies such as Ahmed et al. (2022) for reinforced concrete buildings or Aldemir, Guvenir, and Sahmaran (2020) for masonry structures. The output of the regression analysis is a set of coefficients, and , which quantify the influence of each input indicator, , on the probability of failure. To estimate these coefficients, the use of Generalised Linear Models (GLMs) statistical technique is proposed (Dobson and Barnett 2018; McCullagh and Nelder 1998). This technique is used to model the linear relations between a response variable and a set of explanatory variables. Compared to classic linear regression models, GLMs do not assume that the response variable is normally distributed, allowing it to assume various distributions as long as they belong to the exponential family, e.g., normal, Poisson, and gamma. Moreover, GLMs do not rely directly on a linear relationship between the dependent and independent variables. Instead, they introduce different functional transformations, referred to as link functions, that enable the modelling of a non-linear relationship as a linear one. All link functions, therefore, allow the predicted values to be transformed to a variety of distributions according to the general form:

The selection of the appropriate link function depends on the probability distribution type at different PGA levels. Once the suitable link function is defined, and coefficients are estimated by minimising the difference between the observed values and the model’s predictions, through an iterative process. The output consists of a set of fitted coefficients, one for each i-th PGA level, that quantify the effects of geometric indicators on the probability of failure (or exceeding a selected damage state).

Vulnerability index definition

The final step of the proposed framework is to define a proper vulnerability index . Having obtained a set of and coefficients for each PGA, the challenge is to derive an index that does not depend on a single intensity level but rather reflects the overall vulnerability trend across a range of seismic intensities. Different strategies can be adopted to achieve this, depending on the specific goals and data availability. One possible approach is to integrate or aggregate the and coefficients across selected intensity ranges, thereby obtaining an intensity-range-dependent index. Another possibility is to choose a reference PGA level that reflects desired conditions. The choice depends on the intended application, the available input data, and the desired balance between functionality and accuracy.

Finally, to make interpretable and bounded, the inverse transformation of the link function used is applied:

This transformation yields a value in [0,1], which can be interpreted as an approximate probability of triggering the failure mechanism, under a given combination of geometric indicators and PGA values. Thus, the index obtained can describe vulnerability as the aggregate effect of church resistance (i.e., the combination of geometric indicators) and seismic intensity (i.e., PGA values).

This section presents the application of the methodological framework to the 12 archetypes identified by the authors in Del Carlo et al. (2025) for the one-nave church typology, with the objective of defining the seismic vulnerability index. Section 3.1 describes how pushover analyses are conducted for each archetype model, considering material property variability through sampling. The resulting capacity curves capture the seismic response of the archetypes, enabling the determination of displacement capacity points. As the primary objective is to explore the relationship between geometric indicators and the probability of failure, a single limit state corresponding to the global ultimate state of the church is considered. Global ultimate state refers to the attainment of a state of damage in a substantial portion of the church such that the church's serviceability is prevented. Section 3.2 introduces the computation of seismic displacement demand using the N2 method and details the derivation of the fragility points and curves, accounting for geometric and material variability. Section 3.3 illustrates how the probability of exceeding the damage state is linked to representative geometric indicators using the GLMs and the relative influence of each indicator on the estimated probability. The result is a vulnerability index capable of conveying the seismic physical vulnerability by combining geometric configurations and seismic intensity levels.

Capacity analysis

Pushover analysis is carried out in the longitudinal (+ Y) direction for each archetype (Fig. 2a). Gravity loads are applied in the first loading step, followed by incrementally applied horizontal forces proportional to the structure's mass, until the analysis stops due to model collapse or solution divergence during the iterative process. Only material nonlinearity is considered in the analysis. Notably, while recognising that multi-directional analysis would provide a more complete picture of the seismic behaviour, the study was limited to one direction due to the exploratory nature and novelty of the proposed methodology. Given the need to implement and test the approach for the first time, the focus was on establishing a scalable and reproducible workflow rather than fully characterising the seismic behaviour of the archetypes. Moreover, several studies (Ceroni et al. 2022; Cescatti et al. 2020; Di Meo et al. 2023) analysing damage data at a regional scale have demonstrated that the most recurring mechanisms for single-nave church typology are related to the macro-element of the façade (overturning of the façade or its upper part) and the macro elements of the apse (apse overturning, in-plane apse mechanism, and apse vault mechanism). These failure modes are predominantly activated under longitudinal seismic action. Therefore, the choice of the Y-direction reliably reproduces the damage patterns observed in post-earthquake surveys. Furthermore, applying the action in the positive longitudinal direction ensures the consideration of the most critical condition for the rear façade, while the effect remains symmetric for the main façade.

The analysis employs a simplified Finite Element (FE) model approach, developed using DIANA 10.7. A FE model is developed for each archetype, with the geometric characteristics detailed in Del Carlo et al. (2025), Fig. 5. The masonry walls, apse, and roof slab are modelled using a combination of eight-node quadrilateral (CQ40S) and six-node triangular (CT30S) iso-parametric curved shell elements. These elements ensure an accurate representation of the curved geometry while effectively capturing the structure's nonlinear behaviour. The timber roof beams are modelled using three-node, three-dimensional Class-III beam elements (CL18B). These beam elements incorporate six degrees of freedom per node, enabling realistic simulation of structural response. The final mesh sizes selected are 0.30 m for the shell elements and 0.50 m for the beam elements, with both quadrilateral/hexahedral and triangular/tetrahedral mesh types applied based on geometric needs. Triangular elements are used in regions with curved or irregular shapes, such as the semicircular apse and the triangular gable, where quadrilateral elements would be less effective in accurately fitting the geometry and maintaining mesh quality. The base of the structure is assigned fully fixed boundary conditions, restricting both translational and rotational movements to simulate the foundation’s restraint conditions, given the absence of detailed information on soil–structure interaction. The roof is not connected to the gable shell elements, allowing the activation of out-of-plane mechanisms in the façade walls and therefore ensuring a more realistic response under seismic action. Hinged connections are defined between the strut beams and longitudinal walls, allowing rotation perpendicular to the walls while maintaining rigid translational constraints in other directions. This approach accurately replicates realistic wood beam-to-wall interactions, capturing the semi-rigid behaviour observed in historical masonry structures. It is worth noting that, in this preliminary phase, orthogonal walls are modelled as fully connected, thus not accounting for the actual connection quality, despite its important role in the out-of-plane stability of façades. Given the novelty of the proposed methodology and the lack of detailed, systematic information on the quality of wall-to-wall connections in the surveyed buildings, it would not be easy to make reasonable assumptions or reliably parametrise this aspect. This choiche should therefore be regarded as a simplification made to enable the overall methodological development rather than to reproduce the exact local behaviour.

Therefore, pushover analyses are performed by varying the material properties. Specifically, from performed surveys, four recurrent materials were previously identified (Table 1). Each material is assigned a code, moving from MAT1 – Rubble masonry, to MAT 4 – Squared stone masonry, as detailed in Table 1. Since no experimental data were available, their mechanical characteristics were extrapolated from the reference values provided by the commentary to the Italian National Code — NTC 2018 (D.M.17.01.2018. 2018), Table C8.5.I.

In the FE models, masonry is treated as a homogeneous material, and linear and nonlinear behaviour is modelled using the Total Strain Fixed Crack Model (TSCM) (Vecchio and Collins 1986). Based on the TSCM, the compressive strength () and the Young’s modulus () result as the two primary input parameters, from which all other mechanical properties are derived. Therefore, the tensile strength () is assumed equal to 2% of the compressive strength to represent an almost null tensile capacity while avoiding numerical instability. The value of the compressive fracture energy () is computed, according to Lourenço’s formulation (Lourenço 2010), as a function of the compressive strength:

where d = 1.6 mm. The tensile fracture energy ( is modelled with a value of 50 N/m, a widely accepted value in the literature (Dimovska et al. 2021; Endo et al. 2015; Endo, Pelà, and Roca 2017) and based on sensitivity analysis performed on a representative archetype (Del Carlo, Caprili, and Roca Fabregat n.d.). Finally, the compressive response is assumed to follow a parabolic softening curve, while the tensile response is associated with an exponential softening curve under tension.

In order to be able to treat the mechanical properties of masonry as random variables with the aim of incorporating material variability into the pushover analysis, and values are sampled from their associated Cumulative Distribution Functions (CDFs). As is commonly recognised, material mechanical properties can be modelled through a lognormal distribution to avoid negative values. Therefore, considering a generic mechanical property X, an estimate of the parameters and of its distribution can be calculated from the minimum and maximum values in the table, using the following equations:

where and are the location and scale, respectively. Figure 6 illustrates the lognormal CDF obtained for and for MAT1 (Table 1).

The Latin Hypercube Sampling (LHS) method (Helton and Davis 2003) is employed to sort a suitable number of samples from the CDFs of and . The LHS procedure, which is computationally less expensive than MCS, uses a stratified sampling method that divides the range of each variable (i.e., the and ) into disjoint intervals of equal probability. Then, one value is randomly selected from each interval. In a two-dimensional problem, as in this case, the values obtained for are paired randomly without replacing them with the values obtained for . It’s important to note that the individual variables must be independent for this sampling technique to achieve the desired results. This assumption is realistic within each material class, given the relatively narrow range of values considered. Moreover, although design codes propose a simple proportional relationship between E and for masonry, experimental evidence indicates that the ratio between these two properties can vary significantly, ranging from below 100 to over 1000, depending on the type and quality of the masonry (Cuadros-Rojas et al. 2024; Drougkas, Roca, and Molins 2015). This high variability justifies their independent sampling within the analysis. According to the literature, a proper number of values to be sorted is slightly more than twice the number of variables considered (Celarec and Dolšek 2013).

For each material category (MAT1 to MAT4), the LHS method is applied to derive, altogether, 38 samples of compressive strength and Young’s modulus data pairs. To depict material variation more realistically, a different number of data pairs is drawn from each material, according to their distribution within the sample, i.e., more frequent materials have a higher number of associated data pairs than less frequent materials.

Table 2 presents the values obtained for the 38 identified combinations, with z ranging from 1 to 38. Poisson’s coefficient (ν) and mass density (ρ) are always assumed equal to 0.20 and 1900 Kg/m³, respectively.

Therefore, for each of the 12 archetypes, 38 pushover analyses are conducted. Pushover curves are initially plotted with reference to two control points (see Fig. 7) to monitor the response of the most critical parts of the structure throughout the analysis. The choice of control points and the selection of the representative displacement for the pushover curve are crucial for the subsequent conversion into the capacity curve. In fact, the capacity curve can show significantly different displacement capacities depending on whether the monitored displacement corresponds to a wall that has reached failure and on whether the failure mode is brittle or ductile. In this case, after a detailed study of the damage evolution across the 38 analyses, in-plane failure mechanisms, such as shear failure of the longitudinal walls, are excluded, as they do not appear to characterise the collapse behaviour of the churches during pushover analysis in the positive longitudinal direction. As a result, the selected control points reflect out-of-plane mechanisms, specifically those affecting the main and rear façades of the churches. These control points enable to distinguish distinct structural behaviours under lateral loading.

Figure 8 illustrates, as a proof of concept, the pushover curves obtained for the Archetype ID1, with reference to control points 1 and 2.

Although both control points respond to the same input load increment, their displacements evolve differently. The main façade (dashed red curves) exhibits a brittle response, with a steep initial branch that abruptly transitions to failure. On the contrary, the rear façade shows a more ductile response, as evidenced by the extended horizontal branch in the grey pushover curves. This plateau corresponds to the separation of the façade from the nave and the cracking of the apse vault. Despite the onset of damage, the system retains some capacity to accommodate displacement without a significant increase in acceleration, and the response continues up to the overturning of the upper part of the façade, involving portions of the longitudinal walls and the propagation of cracks in the apse wall and vault.

These behaviours are repeated, with slight variations, across all archetypes and materials. As a proof of concept, Fig. 9 shows the principal strain (E1) at the peak acceleration capacity. As material properties improve, a general increase in peak load factor is observed, progressing from 0.24g (MAT11) to 0.64g (MAT38), confirming the expected enhancement in global capacity. In terms of damage evolution, the main façade consistently shows vertical cracking in all configurations, aligned with an overturning failure mechanism, under the assumption of full connection with the nave walls. However, the extent and intensity of strain concentrations reduce progressively as material strength increases, particularly in the façade corners and around the entrance. Similarly, for the rear façade and the apse vault, crack development is delayed and becomes more localised as material quality increases.

For capacity estimate, two distinct failure modes are therefore identified based on the behaviour observed at the control points. This distinction reflects the inherent complexity of churches as structural systems. Indeed, these buildings are composed of multiple interacting substructures, or macro-elements, whose individual responses collectively govern the global seismic behaviour. Due to this heterogeneous configuration, two separate sets of pushover curves are plotted for different substructures, which, although belonging to the same structure, exhibited two clearly recognisable and distinct failure modes.

Since it is not possible at this stage of the analysis to determine which substructure is more fragile, and thus more likely to govern the overall structural failure, both are analysed separately up to the fragility computation. One substructure is associated with control point 2. It corresponds to the onset of cracking in the apse vault and the initial separation of the rear façade from the longitudinal wall. The displacement capacity for this mode is defined as the point at which the pushover curve transitions into a horizontal branch, indicating a loss of stiffness and damage (Fig. 10). The second substructure corresponds to the activation of the internal overturning mechanism in the main façade, monitored through control point 1. The displacement capacity for this mode is again defined as the point at which the pushover curve transitions to a horizontal branch, or where a drop in horizontal acceleration is observed (Fig. 11).

The displacement capacity () is defined as the maximum displacement observed on the pushover curve. Having generated 38 pushover curves for each archetype, the obtained 38 capacity points represent the displacement capacity as a function of material variability. The KDE method is used to get the PDFs of the capacity-displacement distributions for each archetype. KDE is a non-parametric method for estimating the underlying probability density of a data set. Unlike parametric methods, which assume the data follow a known distribution (e.g., normal, Poisson, or exponential), nonparametric methods do not make such assumptions, making them useful in scenarios where the data do not fit standard models. Here, a Gaussian kernel is used to fit the capacity and demand displacement points. As an example, Fig. 12 illustrates the KDE of the probability density distributions of the , for control points 1 and 2, plotted for archetype ID1.

The PDF of displacement capacity for control point 1 exhibits a sharp peak around 0.02 m, suggesting that the displacement capacity values are more concentrated. On the contrary, control point 2 displays a broader distribution, reflecting greater variability in displacement capacity, likely due to a higher sensitivity to variations in material properties compared to control point 1. This outcome is consistent with the behaviour observed in the pushover curves (Fig. 8), where control point 2 demonstrates a more ductile response. Since ductility is highly influenced by material characteristics, a ductile failure mechanism is expected to result in a wider spread of displacement capacity values. Conversely, control point 1 is associated with a brittle failure mechanism: although the load factor required to activate the mechanism increases, the displacement capacity remains relatively stable across different material properties.

Fragility points derivation

Demand definition

Once the pushover curves are identified, the N2 method is employed to obtain the displacement demand for a set of ADRS. The N2 method is a nonlinear static procedure for seismic analysis that combines a capacity spectrum approach with inelastic demand estimate. The procedure involves transforming the multi-degree-of-freedom (MDOF) system into an equivalent single-degree-of-freedom (SDOF) system using modal properties, typically the first mode response. Seismic demand is represented by ADRS, while capacity is described by a pushover curve converted to an acceleration-displacement format. The intersection of the demand and capacity curves defines the performance point, which provides an estimate of the expected displacement response under seismic load. In the case of the analysed churches, which lack intermediate diaphragms, their structural response can be reasonably approximated as a single-degree-of-freedom (SDOF) system. Consequently, the pushover curves are directly bi-linearised following the N2 method to derive the capacity curves. A set of 24 pseudo-acceleration spectra is obtained from the Italian seismic code, selected to cover a wide range of PGA values. ADRS are obtained from the pseudo-acceleration spectra based on Eq. 6.

being T the period in [s] and and the displacement [m] and pseudo-acceleration [g], respectively. Note that the number of spectra is not fixed, but it was chosen to be twice the number of archetypes to provide a good fit for the fragility functions and cover a wide range of PGA. This choice ensures there are enough points to capture the variability in seismic demand accurately and allows robust definition of the fragility curves for each archetype.

The displacement demand () to which the SDOF system is subjected can be determined efficiently, particularly when its period (T’) exceeds , which is the transition period between the constant pseudo-acceleration and constant pseudo-velocity regions. In this case, the displacement demand can be obtained directly from the elastic response spectrum (Eq. 7). Otherwise, when the period T’ is lower than , the displacement demand exceeds the elastic spectral displacement and must be determined by amplifying the latter according to the factor (Eq. 9).

where m is the SDOF system mass and is the yield force. This process is repeated for each i-th PGA (24 in total), yielding 38 displacement-demand points per spectrum. The obtained demand points define the distribution of displacement demand as a function of material variability for each selected PGA level. Figure 13 illustrates the KDE of the probability density distributions of the , for control point 1 and 2, plotted for archetype ID1. Each curve represents the KDE derived from the intersection of the 38 bi-linearised capacity curves with one of the 24 defined ADRS, reflecting the influence of material variability on displacement demand distribution. Each KDE distribution provides insight into the most probable values of displacement demand for a given PGA level. The presence of multiple peaks in a distribution may indicate that specific variations in material properties affect demand outcomes across the analysed cases.

Fragility curves derivation

Fragility is defined as the conditional probability of failure, i.e., the conditional probability of exceeding prescribed limit or damage states for a given Intensity Measure (IM). In this context, fragility curves are built based on fragility points obtained by comparing the probabilistic distributions of displacement capacity and demand via MCS. From the previous steps, for each archetype, two displacement-capacity distributions were derived, one for control point 1 and one for control point 2, representing the two identified failure modes. As previously mentioned, these failure modes are associated with distinct substructures of the church. However, although these points are linked to specific failure modes, their displacements are governed by the global deformation and overall structural response of the building. Therefore, both control points effectively capture a global measure of vulnerability rather than purely local behaviour, making it essential to discern which of the two exhibits the most fragile response. On the other side, for each of the 24 ADRS spectra, a corresponding displacement demand distribution is computed. As previously discussed, these distributions account for the variability induced by material property differences. The MCS procedure consists of drawing n samples from both the capacity and demand distributions. Note that in this context, the term sample from a distribution refers to the single realisation x, whose probability is . This differs from the more general statistical usage of sample, which typically refers to a collection of realisations . The number of samples needed is determined based on the desired accuracy of the estimated probability of failure. The standard deviation () of the Monte Carlo estimator is given by:

where is the estimated probability of failure, n is the number of samples, and is the coefficient of variation. The approximation in Eq. 10 holds for small . Assuming a relative error of approximately 10%, i.e., , and targeting a minimum meaningful failure probability of, the required number of samples to sort can be approximated as: .

For each pair of sampled values, capacity and demand , an indicator function () is defined, as shown in Eq. 12:

where is the standard cumulative distribution function, and and are the logarithmic median and standard deviation parameters, respectively.

It is worth noting that the proposed approach, unlike traditional performance-based methodologies, allows for the inclusion of material variability within a single fragility curve by generating pushover responses from probabilistic distributions of mechanical properties, thereby providing a more realistic representation of uncertainty in structural capacity.

Figure 14 illustrates the fitted fragility curves for archetype ID1, plotted with respect to control point 1 and 2.

Being median PGA equal to , fragility of control point 1 shows a median PGA of 0.25g, while for control point 2, the median PGA is 0.35g. The result indicates that the failure mode associated with the substructure, as monitored through control point 1, governs the collapse limit state of the church. Therefore, the fragility curve associated with control point 1 is selected to represent the structural performance, as it defines the most critical failure mode.

Regression-based vulnerability modelling

The derivation of fragility points for each archetype served as a basis for the development of a seismic physical vulnerability index through regression modelling. This index is expressed as a combination of geometric indicators, and it is linked to the probability of failure via appropriate transformation functions (link functions) within the GLM technique, as detailed in Section 2.3. Since the objective is to develop an index that expresses the overall vulnerability of the church, only the failure probabilities associated with control point 1 are used to define the vulnerability index.

GLMs are based on a mathematical framework that assumes the dependent variable’s distribution comes from the exponential family; therefore, a binomial distribution is taken for the dependent variable (.) The model formulation described in Eq. 2 is applied to characterise the transformed as a linear combination of selected geometric indicators. To choose the most appropriate link function and therefore transform to fit into Eq. 2, its distributions across different PGA levels are analysed. Figure 15 presents representative failure probability distributions for low, medium, and high PGA values across the 12 archetypes analysed. Distributions are left-skewed at low PGA values, approximately symmetric at intermediate levels, and right-skewed at high PGA levels.

Among the commonly used link functions presented in Section 2.3, logit function (Eq. 14) is uniformly employed to simplify the , as illustrated in the following.

The quantification of the regression coefficients is done via the Iteratively Reweighted Least Squares (IRLS) method. The IRLS procedure is implemented in Python and used to determine the optimal values of the model coefficients and , which express the influence of each geometric indicators on the probability of failure. Four key indicators, wall thickness, span, height, and length, are considered, as illustrated in Fig. 16, which were previously used to derive the archetypes (Del Carlo et al. 2025).

Twenty-four separate GLM models are calibrated, one for each PGA level, using the 12 archetypes’ geometric data and the associated probability of failure for each PGA level. The result is a detailed characterisation of how the influence of geometric indicators varies with increasing seismic demand. This characterisation is expressed in terms of and coefficients. To fully understand the meaning of these coefficients, it is important to remember that the logit function is based on the concept of odds, defined as:

Since the odds increase when increases, a positive indicates that an increase in the associate indicator increases the probability of failure. On the contrary, a negative indicates that an increase in the associate indicator decreases the likelihood of failure. A value close to 0 implies that the indicator has a reduced effect on vulnerability.

Table 3 summarises the results of the regression analysis conducted for control point 1. Each row reports the selected function, the model performance metrics and the regression coefficients associated with the geometric indicators, for each PGA level. The Mean Absolute Error (MAE) provides a measure of the model accuracy; the smaller the MAE, the better the model’s predictions fit the actual data. McFadden's Pseudo-R², referred to as pseudo-R², quantifies the improvement of the fitted model compared with an intercept-only null model, i.e., the model containing only the intercept () and no regression variables, therefore providing a measure of the goodness-of-fit. Pearson Chi² and Deviance are both measures of residual discrepancy.

The pseudo-R² values increase progressively with PGA, reaching a maximum of 0.81 at 0.207g. This indicates that geometric indicators are highly effective at explaining the probability of failure within this intensity range. For very low (PGA < 0.10g) and very high (PGA > 0.40g) seismic intensities, the pseudo-R² values drop below 0.60, suggesting, as expected, that in these cases the probability of failure is less explained through geometric characteristics. Despite this, the MAE remains consistently below 0.05 across all PGA levels, indicating high predictive accuracy of the models. Moreover, Pearson Chi² and Deviance are reasonably close in magnitude, suggesting that the model is well-behaved, and they remain consistently low across all PGA levels (ranging approximately from 0.002 to 0.24). Given the seven residual degrees of freedom (each model has 12 observations and uses five parameters, i.e., the intercept and four predictors), the result indicates an excellent model fit with no overdispersion. It is worth noting that the Deviance and Pearson’s Chi² values tend to be higher in the range where pseudo-R² is also higher. This is not contradictory: while pseudo-R² is a relative measure of explanatory power, Deviance and Pearson’s Chi² provide absolute measures of dispersion. At low and high PGA values, predicted probabilities are typically close to 0 or 1, resulting in low variability and, consequently, low dispersion measures, even when the model's explanatory power is limited. In contrast, for intermediate PGA values (approximately between 0.10g and 0.40g), the observed probabilities are more dispersed, leading to higher Deviance and Pearson’s Chi² values.

Regarding the regression coefficients (Fig. 17a), wall thickness () and church length () remain close to zero across the full range of PGA values, indicating a limited impact on the selected failure mode. A different behaviour is observed for span () and the wall height (). The latter shows consistently positive values that increase with PGA, indicating, as expected, that taller façades are more vulnerable to the selected failure mode. Span (), on the other hand, becomes increasingly negative with increasing PGA, indicating that wider spans are associated with lower failure probability, particularly at medium PGA levels. This last result may appear counterintuitive. However, it could be influenced by the fixed door geometry across archetypes. Specifically, archetypes with larger spans, but identical door dimensions, may have a lower relative opening-to-wall area ratio, which could result in improved seismic capacity compared to archetypes with narrower spans. Notably, while the trends in wall height, wall thickness, and length align with expectations, the correlation between span and reduced vulnerability should be interpreted with caution, as it may be significantly affected by modelling assumptions and warrants further investigation. The intercept () steadily increases with PGA (Fig. 17b), independently of geometry, reflecting the general trend of growing vulnerability at higher seismic intensity levels, as expected. The shape of the curve also mirrors the fragility transition zone, where the failure probability evolves from low to high values as seismic intensity increases.

While the regression models demonstrate good performance across the considered PGA levels, it is essential to acknowledge that the analysis is based on a relatively limited number of observations (12 archetypes), which limits the regression's statistical power and reduces the model’s capacity to generalise trends. Additionally, the geometric characteristics of archetypes exhibit limited variation: two values for the length, span, and height, and three values for the wall thickness, resulting in limited variability in the input space. This homogeneity could affect the model's ability to detect strong relationships between geometry and failure probability, particularly when more complex or localised structural mechanisms are involved.

Vulnerability index definition

Despite the inherent limitations discussed above, the results obtained were considered sufficiently robust to support the derivation of a seismic vulnerability index. Indeed, while the regression models do not achieve perfect predictive performance, they consistently capture meaningful relationships between geometric indicators and failure probability. To translate these findings into a practical and scalable vulnerability index, three ranges of PGA values within the interval 0.10g ≤ PGA ≤ 0.40g are identified. The interval corresponds to pseudo-R² values over 0.60, out of which, as discussed, the probability of having or not having failure is less explained through geometric characteristics. Ranges are selected based on the trend of the logit model’s pseudo‑R² as a function of PGA (Fig. 18). Two clear ‘elbows’ appear at approximately 0.16 g and 0.28 g, where the model’s explanatory power first accelerates and then begins to recede (while remaining above 0.60). Within each of these ranges, the mean values of the regression coefficients, including both the intercept and the geometric weights , are computed from the corresponding values of the regression models (Table 4).

The use of mean values for the regression coefficients is intended to provide an estimate of the central tendency of the vulnerability index within a selected range of PGA values where the model shows strong explanatory power. From an engineering perspective, this choice does not explicitly introduce conservatism, as it does not favour worst-case or upper-bound scenarios. If a more conservative estimate of vulnerability is required, alternative aggregation strategies, such as using the 95% quantile or the maximum coefficients, could be considered. Given that this represents a first attempt at introducing the index and considering the relatively small number of observations and the modelling assumptions discussed earlier, the use of a central tendency appears to be a safe and balanced choice. The resulting PGA-specific vulnerability index , offers a simplified yet informed representation of seismic vulnerability central tendency, that can support regional-scale analysis and prioritisation processes.

Being and directly extracted from the regression, their combination can be positive or negative, making them difficult to interpret or compare across buildings. To overcome this, a sigmoid transformation is applied, leading to an intuitive and interpretable vulnerability index, bounded between 0 and 1:

The sigmoid function is the inverse of the logit function, used to transform the probability of failure into the (Eq. 14) for the regression model. Therefore, represents an indication of the probability of obtaining the failure mechanism described by control point 1, within specific ranges of seismic intensity (PGA). The proposed formulation is applied to a set of 44 one-nave churches, the is computed for each church individually using its actual geometric characteristics, namely wall thickness, span, length, and height, as collected during the in-situ surveys.

The result is illustrated in Fig. 19. Please note that the PGA values used for computation correspond to those with a 10% probability of exceedance in 50 years, on rigid ground conditions, in accordance with national seismic regulations (OPCM 28 aprile 2006 n. 3519 2006). Whether information on local soil conditions and amplification factors is available for the area where a church is located, a more accurate PGA value should be used in the computation of . This would enable site-specific analysis, thereby improving the reliability of the vulnerability estimate.

The present work was developed in response to the growing need for coherent and integrated methodologies capable of addressing the complex and multifaceted risks affecting masonry churches. The primary objective was to develop a scalable and reproducible methodology for estimating seismic vulnerability using an index-based approach, compatible with regional-scale data availability.

The regression analysis proposed prove efficient in providing a vulnerability index that expresses each church's overall seismic vulnerability as the probability of experiencing a selected failure mode, based on easily measurable geometric indicators. By avoiding the need for detailed material or structural characterisation, the index offers a simplified yet effective tool for preliminary regional-scale vulnerability analysis. This could be particularly useful in contexts where only limited information is available, enabling vulnerability classification based on a small set of basic and measurable geometric indicators. These indicators can often be acquired through tools such as GIS or existing databases, avoiding costly and time-consuming on-site survey campaigns.

The derivation of representative archetypes constituted a key step enabling the determination of geometric indicators and the NSA performing. By defining recurring typological patterns and grouping churches into archetypes, it was possible to capture the variability of the building stock while avoiding the need to model each structure individually. The 12 archetypes identified for the one-nave church typology provided a statistically representative sample of the surveyed one-nave churches, serving as a robust foundation for the definition of fragility functions and the subsequent regression analysis. With reference to the latter, the choice to aggregate the model parameters ( and coefficients) across specific ranges of PGA, using the mean value rather than conservative bounds, reflects a reasoned balance between simplicity, data limitations, and interpretability.

Despite the positive results, the index constitutes the first attempt to define a relationship between geometric indicators, the seismic event intensity, and the vulnerability of the church. Consequently, some limitations should be identified. Firstly, the pushover analyses were performed in the positive longitudinal direction. Although this choice was grounded in post-earthquake damage data from surveyed churches, it inherently limits the vulnerability index to represent few failure modes only. This choice stems from the time-intensive nature of the numerical implementation of the proposed vulnerability analysis method, which was being tested for the first time.

Additionally, the process relies on NSA that, despite being applied to a reduced set of representative archetypes rather than individual buildings, remain computationally demanding and time intensive. The substantial time required is further compounded by the decision to explicitly incorporate material variability by deriving probabilistic distributions of capacity for each archetype. While this represents a significant advancement in seismic physical vulnerability analysis, it comes at the cost of increased analytical effort.

Finally, it is important to remember that the choice to base the regression model on geometric indicators alone, while advantageous in terms of rapid applicability, inevitably omits other influential factors such as construction quality or wall-to-wall connections. Consequently, the developed seismic vulnerability index has to be interpreted as a screening tool rather than a detailed predictive model. Its primary purpose is to inform regional-scale strategies for resource allocation, rather than to serve as a substitute for engineering analysis at the building scale.

Future developments could aim to overcome the aforementioned limitations. For example, surrogate modelling techniques could significantly reduce computational time, enabling the consideration of more loading directions and failure modes. Furthermore, adopting vector regression or machine learning algorithms would allow the inclusion of a larger population of buildings without the need to group them into archetypes, while explicitly accounting for geometric variability based on their actual distributions.