Stochastic Dynamic Analysis of Cultural Heritage Towers up to Collapse

: This paper deals with the seismic vulnerability of monumental unreinforced masonry (URM) towers, the fragility of which has not yet been sufﬁciently studied. Thus, the present paper ﬁlls this gap by developing models to investigate the seismic response of URM towers up to collapse. On mount Athos, Greece, there exist more than a hundred medieval towers, having served mainly as campaniles or fortiﬁcations. Eight representative towers were selected for a thorough investigation to estimate their seismic response characteristics. Their history and architectural features are initially discussed and a two-step analysis follows: (i) limit analysis is performed to estimate the collapse mechanism and the locations of critical cracks, (ii) non-linear explicit dynamic analyses are then carried out, developing ﬁnite element (FE) simulations, with cracks modelled as interfacial surfaces to derive the capacity curves. A meaningful deﬁnition of the damage states is proposed based on the characteristics of their capacity curves, with the ultimate limit state related to collapse. The onset of slight damage-state is characterised by the formation and development of cracks responsible for the collapse mechanism of the structure. Apart from these two, another two additional limit states are also speciﬁed: the moderate damage-state and the extensive one. Fragility and vulnerability curves are ﬁnally generated which can help the assessment and preservation of cultural heritage URM towers.


Introduction
Unreinforced masonry (URM) towers constitute an important element of architectural heritage.Their high seismic vulnerability was witnessed during recent earthquakes in Italy which resulted in many of them collapsing [1][2][3][4][5][6].Apart from the very low tensile capacity of URM, a number of additional critical features of these structures increase dramatically their vulnerability: increased slenderness, high weight/strength ratio, insufficient connection with timber diaphragms, inhomogeneous nature of materials, complex constructive stages, lack of rigid diaphragms, presence of vaulted systems, and progressive material degradation due to high compressive stresses, ageing, and environmental effects.
URM towers subjected to ground shaking usually collapse due to overturning of their most critical part, as a consequence of the insufficient connection between perpendicular walls and critical adjacent structural members [6][7][8].For instance, the disintegration of the belfry when the uppermost part hosts bells, owing to the overturn of the piers, can cause the collapse of the whole structure [9].Other possible collapse mechanisms relate to diagonal cracks or, occasionally, vertical ones [10].The formation of cracks occurs mainly in zones with high principle tensile stresses which force the construction to be split in parts.When a structure is cracked and there is no other link between the adjacent members, its dynamic response resembles rocking oscillation.This type of vibration is very frequent in URM towers and drives out-of-plane collapse [11,12].The development of such cracks and the initiation of a collapse mechanism may occur even for very low loadings and may be difficult to notice.
Vulnerability analysis of historical towers subjected to dynamic loads is a challenging task for two main reasons: (i) masonry's non-linear (NL) behaviour is hard to rigorously model due to its complexity, and (ii) accurate material idealisation is not straightforward due to variation in the mechanical properties of the structure.Hence, a number of approaches have been proposed to yield reliable numerical tools capable of precisely representing the seismic behaviour such as micro-models, e.g.[13], discrete models, e.g.[14], and macro-elements together with homogenisation approaches, e.g.[15][16][17].Traditional masonry buildings are not generally able to show sufficient seismic resistance due to ageing of materials and intrinsic inadequacies including: poor tensile strength of masonry, high weight/strength ratio, insufficient connection with timber diaphragms, and the inhomogeneous nature of materials.Evidence, in some cases, of constructions with an intention to resist seismic activity in areas of high seismicity by incorporation of timber and steel elements in the walls can be found only locally [18][19][20][21][22][23], but, even in these cases, not following a rational design.Sophisticated techniques involving dynamic monitoring allow an estimation of the actual state of masonry buildings and towers [24][25][26][27][28].
To prioritise the needs of structural strengthening for preserving architectural heritage, their vulnerability should be assessed [29] and, along with risk scenarios, possible losses should be evaluated.In this framework, our aim is a statistical treatment of the vulnerability by generating fragility curves as a fast tool for assessing the seismic capacity of masonry towers.To this end, we studied the towers' population and selected a representative group of towers to estimate their average response characteristics in terms of expected damage.The basic features of the selected case-study towers are presented in Section 2. Limit analyses were carried out to define, firstly, the collapse mechanism.Section 3 presents the theoretical framework of limit analysis as well as the collapse mechanisms considered.Once the critical mechanism was identified, finite element (FE) simulations and analyses were carried out to investigate the seismic performance of the cracked structure with a gradually increasing intensity up to collapse of several time-histories.In each step the maximum displacement at the top is measured.In Section 5 the considered time-histories are presented and a parametric analysis for the material properties is performed.Then, in Section 6, NL explicit dynamic analyses for the case-study towers are presented and the displacement capacities are estimated.Using the properties of the collapse mode at hand, we converted these displacements into spectral quantities.Proposing a new definition of the damage thresholds for rocking structures, we generated a set of fragility curves to describe damage and out-of-plane collapse expressed in terms of spectral displacement.

Selection of Case-Study Towers
Mount Athos, in Greece, has enjoyed self-administrative rights for more than one millennium.It was formally established in the territory of the Byzantine Empire in 972 A.D. Owing to its uniqueness and its abundance in monumental structures covering several centuries, it has been declared a world heritage site by UNESCO [30].The monastic compounds are composed of-apart from the churches-fortification walls, buildings with various uses (e.g., living, storage, etc.), and towers.These towers were erected as early as the 10th century A.D., mainly as part of the defence system or to host bells [31].Today, there exist more than a hundred towers and eight of them were selected for a thorough study as representatives.
The case-study towers include three defence towers and five campaniles (Figure 1) whose basic characteristics are listed in Table 1.The defence walls follow a standard style and the selection was based on their height and their slenderness ratio.The bell-towers' architecture varies more and, for better representation, five of them were selected.
The case-study towers include three defence towers and five campaniles (Figure 1) whose basic characteristics are listed in Table 1.The defence walls follow a standard style and the selection was based on their height and their slenderness ratio.The bell-towers' architecture varies more and, for better representation, five of them were selected.
The tallest tower under investigation is that of Caracallou Monastery, being 27.75 m tall with a slenderness ratio of 3.08.Caracallou Monastery was established in the beginning of the 11th century A.D. The investigated tower (Figure 1a) was erected in the mid-15th century and is equipped with embrasures (meurtrières) and loopholes while the openings along the height of the structure are very small, forming a robust structure.The horizontal cross section of the tower is square with sides 9.0 m long.The bell tower of Philotheou Monastery (Figure 1d) is among the oldest towers studied herein, built in the 11th century A.D. It has suffered several incidents and, as a consequence, has undergone interventions and partial reconstructions over the centuries.It is 24.90 m tall and its width measures 5.5 m.Hence, its slenderness ratio is 4.5, lower than that of the Vatopaidion tower.The stone masonry was carefully constructed with mortar joints as thick as 7 cm locally.The walls are approximately 1.15 m thick at the base, decreasing only slightly (by around 10 cm) at the top.The belfry is 4.5 m high but the openings are quite small, resulting in a robust structure.The tallest tower under investigation is that of Caracallou Monastery, being 27.75 m tall with a slenderness ratio of 3.08.Caracallou Monastery was established in the beginning of the 11th century A.D. The investigated tower (Figure 1a) was erected in the mid-15th century and is equipped with embrasures (meurtrières) and loopholes while the openings along the height of the structure are very small, forming a robust structure.The horizontal cross section of the tower is square with sides 9.0 m long.
The second tallest tower of the study group is that of the Koutloumousiou Monastery (Figure 1b) dating back to the 14th century.Its height is 26.90 m, while its almost square base has sides equal to 5.24 m.The square section of the tower decreases with height presenting an inclination of 10‰ towards the centre of the square.The wall thickness is 1.25 m.
A slenderer tower is the campanile of Vatopaidion Monastery with a slenderness ratio of approximately 5.7 and height 25.55 m (Figure 1c).It dates back to 1427 A.D. The height of the masonry structure without the roof is 21.0 m.The footprint of the tower is approximately square with sides 4.5 m long.
The bell tower of Philotheou Monastery (Figure 1d) is among the oldest towers studied herein, built in the 11th century A.D. It has suffered several incidents and, as a consequence, has undergone interventions and partial reconstructions over the centuries.It is 24.90 m tall and its width measures 5.5 m.Hence, its slenderness ratio is 4.5, lower than that of the Vatopaidion tower.The stone masonry was carefully constructed with mortar joints as thick as 7 cm locally.The walls are approximately 1.15 m thick at the base, decreasing only slightly (by around 10 cm) at the top.The belfry is 4.5 m high but the openings are quite small, resulting in a robust structure.
Protaton tower (Figure 1e) is one of two towers with the smallest horizontal area (20.25 m 2 ) amongst the case-study towers and was first erected in 1534 A.D. A characteristic of this tower is the presence of a bi-part foundation structure forming an arch and serving as an entrance gate.
The third defence tower under investigation is Dionysiou tower (Figure 1f), most probably also erected during the 14th century.It is shorter than the previous two with a height of 23.16 m.The tower geometry is dominated by strict symmetry with a footprint of 49 m 2 and very few openings.
Xenophontos Monastery was established in 998 A.D. but the campanile construction started around 1820-1830 and was completed on 1 March 1864.It is a more complex structure than the others with a square footprint of 81 m 2 and a height of 19.50 m.The walls are made of stone masonry.Internal levels are made of timber beams and planks forming five storeys (Figure 1g).Tie rods exist at the upper two levels which host the bells.The hip roof is also wooden.
The campanile of Iveron Monastery (Figure 1h) was destroyed several times due to invasions or natural catastrophes such as fires, etc.The last restoration took place around the end of the 19th century and the beginning of the 20th.Its height is 20.52 m.
The main architectural characteristics of the case-study towers are summarised in Table 1.

Material Properties of the Towers
The construction date of the case-study towers varies substantially from as early as the 14th century to the 19th century.Despite the large time span, masonry materials for the towers were extracted from common quarries due to the abundance of sources in the area and, as a result, common material properties are shared [32].The inspection of masonry also revealed similar construction techniques among them.It is noted, though, that despite similar construction techniques, a different response to similar loading conditions should be expected; this is due to the fact that each towers' mechanical characteristics are largely affected by a number of factors; geometry of blocks and bonding, type of mortars, history of each tower, humidity, etc., to name a few.
A detailed elastic analysis should account for the inhomogeneous nature of masonry, normally assumed to be an orthotropic material.However, for large structures with very thick walls, the assumption of a homogeneous material can often be sufficient, with an elastic modulus E = 2.5 GPa, a Poisson ratio v = 0.15, and a specific weight γ = 23 kN/m 2 being realistic assumptions [23,[32][33][34].
A parametric analysis was carried out to investigate the influence of the material properties (Section 5.1).

Failure Mechanisms
The critical collapse mechanism of a structure can be identified by limit analysis.The out-of-plane safety factor λ is estimated by applying the lower bound theorem for all the possible collapse mechanisms simulating the inertia effect with equivalent static forces (static approach) [8,33,35,36].The critical ratio λ of the horizontal actions needed to overturn the tower will be the smallest for the studied collapse mechanisms and can be evaluated through the principle of virtual works, as the ratio of the work of the internal actions over the work of the inertia forces [32].The latter can be expressed as the product of a load multiplier λ on the gravitational forces which results in the expression: In Equation ( 1), W i,x is the virtual work of the action i in the horizontal direction (x) and W j,y is the virtual work of the action j in the vertical direction (y).The actions included in Equation ( 1) are either gravitational forces or frictional ones developed during the sliding of the parts.The gravitational forces take part in both the numerator and the denominator while the frictional forces are only included in the numerator (internal forces).
The possible local collapse mechanisms for these types of towers [32] are presented in Figure 2: (i) overturning of the tower due to a horizontal crack at the base of the structure (Figure 2a); (ii) separation of the walls along a vertical line from bottom to top (Figure 2b); (iii) diagonal cracking of masonry and overturning of the critical triangular part (Figure 2c); and (iv) overturning of the piers of the belfry and local dislocation (Figure 2d).
turn the tower will be the smallest for the studied collapse mechanisms and can be evaluated through the principle of virtual works, as the ratio of the work of the internal actions over the work of the inertia forces [32].The latter can be expressed as the product of a load multiplier λ on the gravitational forces which results in the expression: In Equation ( 1),  , is the virtual work of the action  in the horizontal direction  and  , is the virtual work of the action j in the vertical direction  .The actions included in Equation ( 1) are either gravitational forces or frictional ones developed during the sliding of the parts.The gravitational forces take part in both the numerator and the denominator while the frictional forces are only included in the numerator (internal forces).
The possible local collapse mechanisms for these types of towers [32] are presented in Figure 2   To compare the load multipliers λ of the various mechanisms an equivalent singledegree-of-freedom (SDOF) system [37] should be defined as each of the mechanisms involves a different mass.Using the properties of the SDOF system of each mechanism, the load multiplier k is converted into spectral acceleration.The equivalent mass involved in the considered mechanism is given by the following equation [38]: where M * is the modal mass of the considered mechanism, P i are the local weights of the parts, and the modal shape components are ∆ x,i = ϕ 1,i , with ϕ 1,i being the modal displacement normalised to 1. Assuming that the deformation of masonry is negligible compared with the displacement due to rotation, all points are tilting by the same angle θ (see Figure 3) and the horizontal projection of the displacements ∆ x,i can be assumed linear along the height of the wall.The effective mass ratio is the ratio of the modal mass to the total one given by:

Rocking Response
In the previous section the identification of the collapse mechanism was achieved by limit analysis.Once this is defined, the dynamic response to a ground displacement timehistory should be estimated.If the inclined rigid block of Figure 3 with height h and thickness b is let free to oscillate, its response can be described from Equation ( 6) [46], assuming that it will not slide: where θ is the angle of inclination, W its weight, R the distance of the rotation point from the centre of mass, I0 the polar moment of inertia, and the tangent of α is the ratio b/h.The weight from the upper structure is P while the horizontal forces are assumed inertial forces and thus, connected to W and P by the coefficient λ.The response can be given in terms of displacements δx and δy, or in terms of rotation θ.Angle β and radius R are defined from the points of rotation and the centre of mass as shown in Figure 3.Typically, the rotation that a URM block can undergo does not exceed 20°.In such a case, Equation ( 6) can be simplified to: where p = (WR/I0) 1/2 .The solution of this 2nd order non-homogeneous differential equation is given by: ( ) where θ0 is the initial rotation.The error introduced by approximating Equation ( 6) through Equation ( 7) and its exact solution given by Equation ( 8) is small (about 2%).By assuming θ = 0, i.e., passing through the equilibrium point, then the time will be t = T/4 where T is the period of the oscillation as a total period would be the time required to The spectral acceleration at the base point of the equivalent SDOF system needed to initiate the rocking of the rigid body is given by the maximum value of coefficient λ normalized by the effective mass ratio according to the following expression: If the level of the base point of the collapse mechanism is not at the base of the tower, as in the case of belfries (Figure 2d), then the spectral acceleration capacity should be estimated considering the filtering effect of the substructure.The underlying structure acts as a filter of the base motion modifying the frequency characteristics of the base excitation given by the transfer function H(T s , T 1 ).A number of expressions have been proposed to estimate the floor response spectrum (FRS) given the base response spectrum and the dynamic characteristics of the structure (e.g., [39][40][41][42]).
A simplified expression of the transfer function H(T s , T 1 ) related to the modal period of the rocking system T s and the main system T 1 to estimate the filtering effect of the underlying structure is given by the next expression [43]: In Equation ( 5), ψ 1 (Z) is the shape of the fundamental mode of vibration of the building in the direction considered, normalised with the displacement at the top of the building, and ξ s and ξ 1 are the modal dampings of the rocking and the main system, respectively.For purely rocking systems, the damping was estimated to equal 3% [38].As here friction is present, the damping is assumed to be 5%, equal to that of typical masonry buildings [44] and, therefore, their ratio ξ s / ξ 1 is equal to 1.A reasonable approximation of the fundamental mode of masonry in regular buildings is given by the ratio (Z/h), in which h stands for the height of the structure measured and Z stands for the elevation of the centre of gravity of the mechanism.In Equation ( 5), γ 1 is the corresponding modal participation factor which depends on the number of internal floors n [45] as γ 1 = 3n/(2n + 1).The results of the analysis are presented in Section 6 with those from FE simulations.

Rocking Response
In the previous section the identification of the collapse mechanism was achieved by limit analysis.Once this is defined, the dynamic response to a ground displacement time-history should be estimated.If the inclined rigid block of Figure 3 with height h and thickness b is let free to oscillate, its response can be described from Equation ( 6) [46], assuming that it will not slide: where θ is the angle of inclination, W its weight, R the distance of the rotation point from the centre of mass, I 0 the polar moment of inertia, and the tangent of α is the ratio b/h.The weight from the upper structure is P while the horizontal forces are assumed inertial forces and thus, connected to W and P by the coefficient λ.The response can be given in terms of displacements δx and δy, or in terms of rotation θ.Angle β and radius R are defined from the points of rotation and the centre of mass as shown in Figure 3.Typically, the rotation that a URM block can undergo does not exceed 20 • .In such a case, Equation ( 6) can be simplified to: where p = (WR/I 0 ) 1/2 .The solution of this 2nd order non-homogeneous differential equation is given by: where θ 0 is the initial rotation.The error introduced by approximating Equation (6) through Equation ( 7) and its exact solution given by Equation ( 8) is small (about 2%).By assuming θ = 0, i.e., passing through the equilibrium point, then the time will be t = T/4 where T is the period of the oscillation as a total period would be the time required to return to the initial inclined position.The period vs. the ratio θ/α is plotted in Figure 4 where it is seen that the rocking response varies with θ/α-showing a dependence on the geometric properties (α = b/h)-as well as with θ.
The solution provided by Equation ( 8), however, does not account for either the elastic deformation of the block or for the possibility of sliding during the oscillation.A more precise estimation is carried out using interface elements in a finite element (FE) analysis [47].In Figure 5, the sequential positions are shown and a permanent sliding appears at the end of the first period, assuming a friction coefficient of 0.5 at the interface.This friction coefficient is slightly higher than the 0.4 value proposed by EC6, e.g., in [48,49].It also shown that the impact of its passing through the equilibrium position results in loss of energy, which leads to a gradual decrease in the oscillation amplitude [38].The solution provided by Equation ( 8), however, does not account for either the elastic deformation of the block or for the possibility of sliding during the oscillation.A more precise estimation is carried out using interface elements in a finite element (FE) analysis [47].In Figure 5, the sequential positions are shown and a permanent sliding appears at the end of the first period, assuming a friction coefficient of 0.5 at the interface.This friction coefficient is slightly higher than the 0.4 value proposed by EC6, e.g., in [48,49].It also shown that the impact of its passing through the equilibrium position results in loss of energy, which leads to a gradual decrease in the oscillation amplitude [38].The solution provided by Equation ( 8), however, does not account for either the elastic deformation of the block or for the possibility of sliding during the oscillation.A more precise estimation is carried out using interface elements in a finite element (FE) analysis [47].In Figure 5, the sequential positions are shown and a permanent sliding appears at the end of the first period, assuming a friction coefficient of 0.5 at the interface.This friction coefficient is slightly higher than the 0.4 value proposed by EC6, e.g., in [48,49].It also shown that the impact of its passing through the equilibrium position results in loss of energy, which leads to a gradual decrease in the oscillation amplitude [38].The comparison of the analytical solution of Equation ( 8) and the numerical one results in Figure 6: the elastic body continues to oscillate further as the bounding impact leads to a lower loss of energy.The difference in the response of the analytical and the numerical model is due to two factors: (i) the numerical model considers an elastic body, and (ii) the coefficient of restitution estimated by Housner [46] overestimates the energy loss due to bouncing [50][51][52].The comparison of the analytical solution of Equation ( 8) and the numerical one results in Figure 6: the elastic body continues to oscillate further as the bounding impact leads to a lower loss of energy.The difference in the response of the analytical and the numerical model is due to two factors: (i) the numerical model considers an elastic body, and (ii) the coefficient of restitution estimated by Housner [46] overestimates the energy loss due to bouncing [50][51][52].The comparison of the analytical solution of Equation ( 8) and the numerical one results in Figure 6: the elastic body continues to oscillate further as the bounding impact leads to a lower loss of energy.The difference in the response of the analytical and the numerical model is due to two factors: (i) the numerical model considers an elastic body, and (ii) the coefficient of restitution estimated by Housner [46] overestimates the energy loss due to bouncing [50][51][52].

Variations in Material Properties
A parametric analysis was carried out varying the material properties of masonry used in Section 2.2.In Table 2, various combinations of masonry units and mortars with varying mechanical properties are taken into account, resulting in different masonry elastic effective properties.This range was based on the relevant literature, e.g., [53][54][55][56][57][58][59][60].The characteristic strength of masonry fk is evaluated applying the Eurocode 6 empirical formula [44] and using the coefficients suggested by [61] for units from stone masonry:  = 0.2 . .

Variability Analysis 5.1. Variations in Material Properties
A parametric analysis was carried out varying the material properties of masonry used in Section 2.2.In Table 2, various combinations of masonry units and mortars with varying mechanical properties are taken into account, resulting in different masonry elastic effective properties.This range was based on the relevant literature, e.g., [53][54][55][56][57][58][59][60].The characteristic strength of masonry f k is evaluated applying the Eurocode 6 empirical formula [44] and using the coefficients suggested by [61] for units from stone masonry: Introducing the material properties of Table 2 (mortar joint f m and units f k strengths), a parametric analysis was conducted regarding a generic tower whose model is shown in Figure 7a.The tower was assumed to be founded on rocky soil and its footprint is a square 4.5 × 4.5 m 2 while its height is 27.5 m.The tower was finely discretised in Abaqus [47] with square plane elements 0.3 × 0.3 m 2 fixed at the base.The material properties of masonry are assumed to be linear; the interaction between the cracked parts at the belfry is ruled by a friction coefficient equal to 0.5.The collapse mechanism is shown in Figure 7b.The influence of the material properties is illustrated in Figure 8 in terms of maximum horizontal displacement at the control point (top) vs. the homogenised (effective) elastic modulus: for an increase 50% in the elastic modulus the maximum top displacement is only changed by approximately 25%.Introducing the material properties of Table 2 (mortar joint fm and units fk strengths), a parametric analysis was conducted regarding a generic tower whose model is shown in Figure 7a.The tower was assumed to be founded on rocky soil and its footprint is a square 4.5 × 4.5 m 2 while its height is 27.5 m.The tower was finely discretised in Abaqus [47] with square plane elements 0.3 × 0.3 m 2 fixed at the base.The material properties of masonry are assumed to be linear; the interaction between the cracked parts at the belfry is ruled by a friction coefficient equal to 0.5.The collapse mechanism is shown in Figure 7b.The influence of the material properties is illustrated in Figure 8 in terms of maximum horizontal displacement at the control point (top) vs. the homogenised (effective) elastic modulus: for an increase 50% in the elastic modulus the maximum top displacement is only changed by approximately 25%.

Variations in Seismic Excitation
A variation in the characteristics of the seismic shaking can appear due to an activation of different seismic faults.To this end, eight different seismic recordings and earthquake spectra were assumed.The displacement spectra of the earthquakes are presented

Variations in Seismic Excitation
A variation in the characteristics of the seismic shaking can appear due to an activation of different seismic faults.To this end, eight different seismic recordings and earthquake spectra were assumed.The displacement spectra of the earthquakes are presented in Figure 9.The set of the considered earthquakes covers a wide area of frequencies with a fundamental period of as little as 0.125 s, for the Whittier Narrows shaking, to 0.675 s for the Landers recording.Their original PGAs (peak ground accelerations) also vary from 0.14 g for the Thessaloniki recording to 0.40 g for the Northridge one.

Variations in Seismic Excitation
A variation in the characteristics of the seismic shaking can appear due to an activation of different seismic faults.To this end, eight different seismic recordings and earthquake spectra were assumed.The displacement spectra of the earthquakes are presented in Figure 9.The set of the considered earthquakes covers a wide area of frequencies with a fundamental period of as little as 0.125 s, for the Whittier Narrows shaking, to 0.675 s for the Landers recording.Their original PGAs (peak ground accelerations) also vary from 0.14 g for the Thessaloniki recording to 0.40 g for the Northridge one.

FE Models
Fragility analysis comes at a high cost.To maintain this cost at an affordable level, 2D models are developed using plane elements for masonry.The façade with the most critical mechanism (Figure 10b,f) or the one with the highest percentage of openings (Figure 10a,c,d,e,g,h) was simulated.As previously mentioned, the main focus is the rocking response due to the presence of cracks.The simulation of the behaviour of towers considers cracks as predefined interfaces following the classical Mohr-Coulomb friction law.These cracks were defined from the limit analysis presented in Section 3. Obviously, these interfaces can open, rotate, and slide.The material laws assumed for masonry are purely

Seismic Fragility of Towers 6.1. FE Models
Fragility analysis comes at a high cost.To maintain this cost at an affordable level, 2D models are developed using plane elements for masonry.The façade with the most critical mechanism (Figure 10b,f) or the one with the highest percentage of openings (Figure 10a,c,d,e,g,h) was simulated.As previously mentioned, the main focus is the rocking response due to the presence of cracks.The simulation of the behaviour of towers considers cracks as predefined interfaces following the classical Mohr-Coulomb friction law.These cracks were defined from the limit analysis presented in Section 3. Obviously, these interfaces can open, rotate, and slide.The material laws assumed for masonry are purely elastic as masonry deformation with respect to the opening of cracks is very small.In order to achieve mesh-size independence of the results, a sensitivity analysis was carried out leading to a discretisation of each tower with approximately 2000 elements.The discretisation and the stress fields for gravitational loads for the towers are shown in Figure 10.

Collapse Analysis
The incremental dynamic analyses (IDAs) were performed by scaling up the recordings presented in Section 5.2 with a step of 0.1 g.NL explicit dynamic analyses were performed in Abaqus [47].In Figure 11 the collapse mechanism of the final step of each tower is shown: there are three towers (a, b, and f) collapsing in diagonal failure (see Figure 2), two in belfry failure (c and e), and the rest (d, g, and h) in a combination of diagonal failure and belfry dislocation.Belfry construction can be from masonry pillars or solid stone columns.Protaton belfry's columns (Figure 11e) overturn before the final collapse.As already stated, the collapse mechanism was defined by performing a limit analysis (see Section 3).Once the mechanism collapses, the analysis stops.
elastic as masonry deformation with respect to the opening of cracks is very small.In order to achieve mesh-size independence of the results, a sensitivity analysis was carried out leading to a discretisation of each tower with approximately 2000 elements.The discretisation and the stress fields for gravitational loads for the towers are shown in Figure 10.The incremental dynamic analyses (IDAs) result in the capacity curves compiled from the maximum displacement vs. PGA, i.e., the intensity measure (IM) of the ground motion for each step.The acceleration of the seismic vibration is transformed into spectral acceleration using the dynamic characteristics of the rocking part and the results are plotted in Figure 12.As can be observed, the IDA capacity curves present a maximum spectral capacity followed by a steep negative slope.This behaviour, with the initial amplification of the capacity, was also noticed experimentally [62].In Figure 12, the capacity curves derived from limit analysis [32,38] are also plotted.The maximum displacement capacities estimated with IDA and limit analysis are generally close.
formed in Abaqus [47].In Figure 11 the collapse mechanism of the final step of each tower is shown: there are three towers (a, b, and f) collapsing in diagonal failure (see Figure 2), two in belfry failure (c and e), and the rest (d, g, and h) in a combination of diagonal failure and belfry dislocation.Belfry construction can be from masonry pillars or solid stone columns.Protaton belfry's columns (Figure 11e) overturn before the final collapse.As already stated, the collapse mechanism was defined by performing a limit analysis (see Section 3).Once the mechanism collapses, the analysis stops.The incremental dynamic analyses (IDAs) result in the capacity curves compiled from the maximum displacement vs. PGA, i.e., the intensity measure (IM) of the ground motion for each step.The acceleration of the seismic vibration is transformed into spectral acceleration using the dynamic characteristics of the rocking part and the results are plotted in Figure 12.As can be observed, the IDA capacity curves present a maximum spectral capacity followed by a steep negative slope.This behaviour, with the initial amplification of the capacity, was also noticed experimentally [62].In Figure 12, the capacity curves

Fragility Analysis
Fragility curves describing the evolution of damage can be based on a (cumulative) lognormal distribution [63][64][65][66].The cumulative lognormal fragility curves are dependent on two parameters: the mean value m i of the spectral displacement and its standard deviation σ i of each limit state i.Statistical simulations such as Monte Carlo analyses can be used to determine the mean value and the standard deviation [8].Lagomarsino [37] suggested four damage states defined on the capacity curves.
The damage states are defined on the capacity curves from IDAs (Figure 12).Four damage states are adopted: (i) slight damage, (ii) moderate damage, (iii) extensive damage, and (iv) collapse.The general capacity curve contains some distinct points (Figure 11): (i) a peak value signifying the onset of rocking, and (ii) an ultimate linear branch which has a slope very close to that of the limit analysis.In accordance with the capacity curve, the damage states are defined as (Figure 13): (i) the onset of DS1, coinciding with the onset of rocking, (ii) the first point of the last linear branch of the capacity curve, marking the onset of DS3, and (iii) DS3, coinciding with the overturn of the structure.The moderate damage state is defined with respect to the previous (DS1) and the next damage state (DS3) as their mean value [37].To derive the capacity curve the mean and the standard deviation of each damage state should be estimated from IDAs.
The standard deviation represents the uncertainty coming from the definition of the damage states β T,ds,i , as well as the variability due to the ground motion β D and the variability due to the response of the structure β C .Such fragility curves are shown in Figure 14 for each tower separately, while Figure 15 shows fragility curves for towers and campaniles with small and large area of openings.In Figure 14, the spectral displacement in the x-axis of the fragility curves has a constant range for comparison.
Buildings 2021, 11, x FOR PEER REVIEW 14 of 22 derived from limit analysis [32,38] are also plotted.The maximum displacement capacities estimated with IDA and limit analysis are generally close.

Fragility Analysis
Fragility curves describing the evolution of damage can be based on a (cumulative) lognormal distribution [63][64][65][66].The cumulative lognormal fragility curves are dependent on two parameters: the mean value  of the spectral displacement and its standard deviation  of each limit state i.Statistical simulations such as Monte Carlo analyses can be used to determine the mean value and the standard deviation [8].Lagomarsino [37] suggested four damage states defined on the capacity curves.
The damage states are defined on the capacity curves from IDAs (Figure 12).Four damage states are adopted: (i) slight damage, (ii) moderate damage, (iii) extensive damage, and (iv) collapse.The general capacity curve contains some distinct points (Figure 11): (i) a peak value signifying the onset of rocking, and (ii) an ultimate linear branch which has a slope very close to that of the limit analysis.In accordance with the capacity curve, the damage states are defined as (Figure 13): (i) the onset of DS1, coinciding with the onset of rocking, (ii) the first point of the last linear branch of the capacity curve, marking the onset of DS3, and (iii) DS3, coinciding with the overturn of the structure.The moderate damage state is defined with respect to the previous (DS1) and the next damage state (DS3) as their mean value [37].To derive the capacity curve the mean and the standard deviation of each damage state should be estimated from IDAs.The standard deviation represents the uncertainty coming from the definition of the damage states  , , , as well as the variability due to the ground motion  and the variability due to the response of the structure  .Such fragility curves are shown in Figure 14 for each tower separately, while Figure 15 shows fragility curves for towers and campaniles with small and large area of openings.In Figure 14, the spectral displacement in the x-axis of the fragility curves has a constant range for comparison.The standard deviation represents the uncertainty coming from the definition of the damage states  , , , as well as the variability due to the ground motion  and the variability due to the response of the structure  .Such fragility curves are shown in Figure 14 for each tower separately, while Figure 15 shows fragility curves for towers and campaniles with small and large area of openings.In Figure 14, the spectral displacement in the x-axis of the fragility curves has a constant range for comparison.The fragility curves shown in Figure 14 adhere to the studied towers or towers with similar characteristics.To this end, three more general sets of fragility curves were generated for towers and campaniles with characteristics within the range of the studied typologies (see Figure 1 and Tables 1 and 2).A distinction is made between towers and campaniles, with the latter further divided into two categories (Table 3): (i) campaniles with small openings' area (openings <10% of external surface) and, (ii) campaniles with large openings' area (openings ≥10% of the external surface).The fragility curves of these three sets are presented in Figure 15: towers appear to have slight or moderate damage for lower intensities, but for higher intensities they can be more vulnerable than campaniles, whose vulnerability increases with the size of openings.The fragility curves shown in Figure 14 adhere to the studied towers or towers with similar characteristics.To this end, three more general sets of fragility curves were generated for towers and campaniles with characteristics within the range of the studied typologies (see Figure 1 and Tables 1 and 2).A distinction is made between towers and campaniles, with the latter further divided into two categories (Table 3): (i) campaniles with small openings' area (openings <10% of external surface) and, (ii) campaniles with large openings' area (openings ≥10% of the external surface).The fragility curves of these three sets are presented in Figure 15: towers appear to have slight or moderate damage for lower intensities, but for higher intensities they can be more vulnerable than campaniles, whose vulnerability increases with the size of openings.The analysis shows that the towers are less vulnerable, due to their less slender geometry, than the bell-towers, whose fragility increases partly due to their belfry.The analysis shows that the towers are less vulnerable, due to their less slender geometry, than the bell-towers, whose fragility increases partly due to their belfry.Moreover, the openings, despite resulting in a less robust structure due to the reduced mass, improve the vulnerability.

Conclusions
A group of representative URM towers was investigated herein to study their seismic performance and provide statistical fragility curves which are useful for their assessment and the prioritisation of interventions.The case-study towers are located on the Athos peninsula (Greece) over an area of approximately 300 km 2 .The proximity of their locations and the use of the same material sources justifies the adoption of similar elastic properties for all of them, despite the differences in their construction dates.The analysis showed that small variations in the mechanical properties do not substantially affect the global response, as the rocking type of failure dominates.The basic assumption of the procedure is a two-step analysis where the collapse mechanism is first identified and then an incremental dynamic analysis (IDA) of the FE model is performed.
The critical collapse mechanism of the towers was estimated by applying limit analysis.For each collapse mechanism, an FE simulation was developed which was subjected to an IDA to derive the capacity curve.The model involves not only NL interactions of the cracks but also linear material properties.
The capacity curves of the towers with cracks are only meaningful in terms of spectral acceleration and displacement.The capacity curve comprises some distinct features: (a) an initial peak denoting the onset of rocking, and (b) a final descending line leading to overturn very close to the line of the limit analysis.Damage states are defined in terms of these features: the onset of rocking is assumed to coincide with the threshold of slight damage, while the two extremes of the final line denote the extensive and complete damage states.Following this procedure, the mean values of the damage state thresholds were derived for each tower and each earthquake excitation, as well as their standard deviations.Further to this analysis, three main tower types were considered: (i) towers, (ii) campaniles with small openings' area, and (iii) campaniles with large openings' area.The mean value and the standard deviations were evaluated.
Using the mean and the standard deviation, a set of vulnerability curves expressed in terms of spectral displacement are proposed which can be used for cultural heritage structures with the characteristics of the three groups.By using these curves, fast and meaningful conclusions can be deducted regarding the risk management of masonry towers.

Figure 1 .
Figure 1.Case-study towers in Mount Athos: (a) Caracallou, (b) Koutloumousiou, (c) Vatopedi, (d) Philotheou, (e) Protaton, (f) Dionysiou (g) Xenophontos, and (h) Iveron.The second tallest tower of the study group is that of the Koutloumousiou Monastery (Figure 1b,) dating back to the 14th century.Its height is 26.90 m, while its almost square base has sides equal to 5.24 m.The square section of the tower decreases with height presenting an inclination of 10‰ towards the centre of the square.The wall thickness is 1.25 m.A slenderer tower is the campanile of Vatopaidion Monastery with a slenderness ratio of approximately 5.7 and height 25.55 m (Figure 1c).It dates back to 1427 A.D. The height of the masonry structure without the roof is 21.0 m.The footprint of the tower is approximately square with sides 4.5 m long.The bell tower of Philotheou Monastery (Figure1d) is among the oldest towers studied herein, built in the 11th century A.D. It has suffered several incidents and, as a consequence, has undergone interventions and partial reconstructions over the centuries.It is 24.90 m tall and its width measures 5.5 m.Hence, its slenderness ratio is 4.5, lower than that of the Vatopaidion tower.The stone masonry was carefully constructed with mortar joints as thick as 7 cm locally.The walls are approximately 1.15 m thick at the base, decreasing only slightly (by around 10 cm) at the top.The belfry is 4.5 m high but the openings are quite small, resulting in a robust structure.
: (i) overturning of the tower due to a horizontal crack at the base of the structure (Figure2a); (ii) separation of the walls along a vertical line from bottom to top (Figure2b); (iii) diagonal cracking of masonry and overturning of the critical triangular part (Figure2c); and (iv) overturning of the piers of the belfry and local dislocation (Figure2d).

Figure 3 .
Figure 3. Rocking response of a URM block.

Figure 4 .
Figure 4. Rocking response of a URM block.

Figure 4 .
Figure 4. Rocking response of a URM block.

Figure 4 .
Figure 4. Rocking response of a URM block.

Figure 6 .
Figure 6.Comparison of the rocking response of URM blocks: rigid and elastic blocks.

Figure 6 .
Figure 6.Comparison of the rocking response of URM blocks: rigid and elastic blocks.

Figure 7 .Figure 7 . 22 Figure 8 .
Figure 7. Simulation of the generic tower and stresses in kPa: (a) gravity loads, and (b) collapse.(a) (b) Figure 7. Simulation of the generic tower and stresses in kPa: (a) gravity loads, and (b) collapse.

Figure 8 .
Figure 8. Variation of the rocking max displacement at the top for varying effective material properties.

Figure 8 .
Figure 8. Variation of the rocking max displacement at the top for varying effective material properties.

Figure 13 .
Figure 13.The definition of the limit states on the capacity curve for Protaton bell-tower and Northridge earthquake.

Buildings 2021 , 22 Figure 13 .
Figure 13.The definition of the limit states on the capacity curve for Protaton bell-tower and Northridge earthquake.

Table 1 .
Architectural characteristics of the case-study towers.

Table 2 .
Variation of material properties.

Table 3 .
Grouping of towers with similar characteristics.

Table 3 .
Grouping of towers with similar characteristics.