Seismic vulnerability assessment against rocking and sliding failure using nonlinear dynamic analysis: Application to the temples of Bagan, Myanmar

Structures 74 (2025) 108584 2352-0124/© 2025 The Author(s). Published by Elsevier Ltd on behalf of Institution of Structural Engineers. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Seismic vulnerability assessment against rocking and sliding failure using nonlinear dynamic analysis: Application to the temples of Bagan, Myanmar Dario Vecchio a,* , Georgios Vlachakis a , Arun Menon b , Paulo B. Lourenço a a ISISE, ARISE, University of Minho, School of Engineering, Department of Civil Engineering, Portugal b Department of Civil Engineering, IIT Madras, India ARTICLE INFO Keywords: Cultural heritage Overturning Shear-sliding Discontinuous interfaces Nonlinear time history analysis Seismic assessment ABSTRACT Built cultural heritage anchors the cultural identities of peoples worldwide and constitutes capital for sustainable development. Maintenance and adequate conservation actions are often absent, increasing the exposure of historical structures to natural hazards such as earthquakes, floods, and the long-term effects of climate change. On August 24, 2016, a M w 6.8 earthquake occurred in Chauk, Myanmar, with an epicentral distance of approximately 40 km away from the UNESCO archaeological site of Bagan. Despite past interventions to strengthen temples in the area, recurring critical damage such as overturning and shear-sliding were observed in situ. This study presents a numerical methodology for a vulnerability assessment on the site, using four temples as examples of different structural typologies. The methodology aims to capture and interpret the observed collapses in situ, as well as to validate a computationally feasible modelling strategy. The numerical models are based on elastic bodies and discontinuous interfaces, aiming to simulate rocking and sliding failure. The temples are large, so two simplifications reducing computational cost are proposed and examined: first, use of three single-degree-of-freedom (SDOF) oscillators as substitutes for the dynamic behaviour of the ground floor; and second, neglect of the local disintegration mechanisms of masonry. The temples were assessed against 11 ground motion records including the 2016 Chauk earthquake, a simulation of the same, and nine physics-based ground motion simulations for a rupture scenario of M w 7.3. The results of the 44 analyses are given in terms of rocking angles and sliding displacements. Collapses and other serious vulnerabilities were found in most of the slender structures of the temples, including the central spires and the small corner stupas placed at different terrace levels. 1. Introduction Built cultural heritage links present and past by preserving ancestral technological solutions and keeping local traditions and cultural values. Moreover, cultural heritage grounds the identity of each society and fosters sustainable development. Lack of maintenance has allowed losses of built heritage to increase due to natural hazards, such as floods, fires, earthquakes, tsunamis, and climate change effects. Therefore, the built cultural heritage stock is subjected to increasing risk [1,2]. Despite the efforts of the international community, more attention needs to be drawn to the conservation of the built cultural heritage and, in particular, to masonry structures [3]. Several studies have focussed on the evaluation of the seismic vulnerability of monuments, temples, and archaeological sites using numerical models. Previous attempts include investigations into Buddhist temples [4], classical ancient temples [5-8], towers [9-12], pagodas [13], mosques [14], and cathedrals [15], among others. The focus of these studies has varied depending on context, scope, and need, in addition to the state of a given monument. Studies have been done on seismic vulnerability [4,5,10,14], the design of strengthening solutions [13], post-earthquake damage observations [16], and the interpretation of existing damage [15,11,12]. Most studies employed advanced numerical models for the seismic dynamic response of monuments [17], such as the finite element method (FEM) [4,15,13,14] or the discrete element method (DEM) [5,10-12,7,8]. The complexity and details of the models varied significantly, and they were usually constrained by computational cost [9]. Small structures with simple geometry permitted mesoscale modelling approaches [10,11], while large structures with complicated geometrical details required macroscale models [15,16]. * Corresponding author. E-mail address: [email protected] (D. Vecchio). Contents lists available at ScienceDirect Structures journal homepage: www.elsevier.com/locate/structures https://doi.org/10.1016/j.istruc.2025.108584 Received 16 January 2024; Received in revised form 22 February 2025; Accepted 24 February 2025 Structures 74 (2025) 108584 2 Due to the variety of construction materials and techniques found in existing temples and monuments, the choice of appropriate constitutive laws is challenging. Whenever feasible, calibration appears to be a sound approach to validate a model [4,15,10,11,14]. It is accepted that the seismic assessment of monumental structures is complex at the engineering and scientific levels. Past damage observations after seismic events have revealed two important factors: the high vulnerability of structural and non-structural elements [18,19] and the low efficiency of strengthening interventions [20-22]. These factors have led researchers to address non-structural elements, such as parapets, chimneys, pinnacles, and archaeological remains, weakly connected to the rest of the structure [23]. It is also common to encounter (slender) elements in a monumental structure resting on soft bases that may show low resistance to partial local mechanisms and experience collapse. This partial failure may cause cultural, material, and economic losses, as well as fatalities. On August 24, 2016, a moderate seismic event occurred in Chauk, Myanmar and affected hundreds of monuments in the archaeological site of Bagan. The site had recently been listed as a World Heritage Site [24], with a tourist industry of 7.5 million visitors per year (before the current political turmoil), according to the Asian Development Bank [25]; Rich &Franck, [26]. The earthquake struck west of Chauk city with a moment M w 6.8, killing 3 people and damaging 10 residential buildings and almost 400 Buddhist monuments [27]. Although the seismicity of the region was well known and despite retrofitting interventions having been implemented by locals or international experts, many of the buildings suffered local failure, damage, or collapse, confirming their deep vulnerability [25,26,24]. Two large projects were initiated to recover the cultural loss and gain information for seismic risk reduction, one led by UNESCO [28] and one led by the Getty Conservation Institute [29] in collaboration with the Myanmar Department of Archaeology and National Museum, part of the Ministry of Culture and Religious Affairs. This paper describes a numerical methodology to investigate four temples within the Bagan archaeological site, aiming for an adequate assessment of seismic vulnerability. These temples were selected as representative of the different structural typologies found in the site, according to seismic behaviour criteria, namely plan symmetry and vertical continuity. In Section 2, the temples are first described through their historical background and their responses to two of the most destructive past seismic activities (the 1975 and 2016 earthquakes). Recent in situ observations have revealed overturning and shear-sliding failures in non-structural elements placed at different heights of the temples, while the central cores have responded elastically to the seismic actions. Based on the damage, Section 3 describes the numerical strategy adopted to conduct the seismic assessment of overturning and shear-sliding failures. To reduce the computational cost associated with the use of nonlinear time history analysis, two main assumptions were introduced. The assumptions included the substitution of the lower part of the structure aiming to reduce the number of degrees of freedom (DOFs), modelling of the structural components as elastic bodies separated by discontinuous contact interfaces, and neglecting masonry disintegration. Prior to the extensive seismic analyses, these simplifications are validated in Section 4. Each temple was next subjected to 11 ground motion records, including the 2016 Chauk earthquake record, a simulation of the same, and nine additional accelerograms created with a broadband ground motion model. The results described in Section 5 capture the occurrence of rocking, overturning, and sliding mechanisms, which align with in situ observations. Section 6 identifies lessons from the analysis of these vulnerable cultural heritage sites. 2. Bagan archaeological site 2.1. Historical background The site at Bagan was founded in the early second century A.D. under the flourishing Pyu culture. Starting in the ninth century, the site was the capital of the Pagan kingdom, which built thousands of religious temples. As many as 3595 of these monuments survive still standing [30] (see Fig. 1). The monumental heritage includes temples, stupas (spires), monasteries, libraries, palaces, and fortifications [30]. According to a UNESCO inventory [31], Bagan religious buildings have been classified as temples and stupas, with their main differences arising from geometry. The temples were built on a square/rectangular plan configuration, with an inner shrine and terrace levels, crowned with solid structures, either square (shikhara) or round plan with curved profiles (bell-shaped dome and spire). Stupas, or spires, were built in a typical conical shape, but varied in size. Standing on a square/rectangular plan, they could be relatively tall (up to 10 m). Unlike temples, stupas may be both structural and non-structural elements. It is common to encounter monuments at the Bagan archaeological site that include smaller stupas sitting above the ground level. Most Bagan temples were constructed with three-leaf peripheral masonry walls of rectangular bricks varying in dimension [31,32]. The external leaf is made of semiregular masonry of running bond, with varying thicknesses of fired bricks and mortar. Three types of mortar are found in Bagan: mud mortar, lime mortar, and organic mortar [33]. The latter denotes a mortar mixture with organic additives, such as starch, natural fibres, proteins, and fatty acids. Despite the existence of a study on the state of conservation of masonry [32], a lack of in situ information persists regarding the mechanical properties of the constituent materials [33]. The present research focusses on four representative monuments of the Bagan archaeological site. The numbering, naming, and sources of historical information on the buildings were primarily taken from the Inventory of Monuments of Bagan [31]. 2.2. Seismicity of the area and past structural damages Myanmar is located south of the Himalayan chain, bordering with the Indian Ocean and several countries of Mainland Southeast Asia. The territory has a complex configuration, and the presence of the Alpide belt exposes the country to the hazards of severe earthquakes and moderate tsunamis [34]. Detailed tectonic and seismic maps can be found in Wang et al. [35] and Thein et al. [34]. According to these findings, the seismicity of the area is defined by the eastern flank of the Arakan Yoma belt, along the Sagaing fault, and along the northern Shan Plateau [36]. Recently, the slip rate of the Sagaing fault has been estimated as 20 mm/year [37-39,36,35].Table 1 collects the largest magnitude earthquakes occurring in Myanmar since the 20th century, where many earthquakes have a moment magnitude of M w ≥7.0. According to the Seismic map of Myanmar, the Bagan archaeological site belongs to Zone IV (Severe), with an expected intensity between VIII and IX using the Modified Mercalli (MM) scale and expected ground accelerations between 0.3 g and 0.4 g. The 1975 earthquake (highlighted in Table 1) affected several buildings, even causing collapses, which resulted in vast economic and cultural losses [4]. A retrofitting project of the monuments was started [41], but never completed, and many monuments were left under partial reconstruction. With the 2016 Chauk earthquake, several previously reconstructed portions were newly damaged or collapsed. A post-earthquake photographic survey revealed that most failures occurred by overturning and shear-sliding at the base of non-structural elements placed at higher elevations [16,41]. These elements were mainly the central bell-shaped dome and spire, the corner stupas, as well as internal statues, pediments, finials, and other similar types of elements, which, due to their geometry, experienced rocking at their base or at a higher plane of weakness. Differently from these non-structural elements, limited damage was observed in the massive structure of the temples, mostly exhibiting minor cracks in the massive masonry walls and vaults. These observations led to the conclusion that the ground floor structure of these temples responded mostly elastically to the D. Vecchio et al. Structures 74 (2025) 108584 3 seismic actions during the 2016 Chauk earthquake. Fig. 2 presents an aerial view of the location of the four temples of interest, together with the rest of the temples in the Bagan Archaeological Zone. The next section presents a detailed description of those temples, with special attention to previous seismic damage. 2.3. Description of the four temples 2.3.1. Temple No. 844: Tha-mu-ti-hpaya Temple No. 844 was constructed in 1260 A.D. It is a large singlestorey temple (Fig. 3) with a solid core and a high-vaulted corridor. The core structure includes a barrel vault over the central corridor and one hipped vault at the east end, over the entrance hall. There are also Fig. 1. Aerial view of the Bagan archaeological site: Source from [29]. Table 1 The major seismic activities that occurred in Myanmar during the 20th century. The two earthquakes that hit the Bagan archaeological site are in bold type. Source from [40]. Date Earthquake name Lat. [◦N] Long. [◦E] M w Focal depth [km] Intensity Epicentral distance*[km] 23-05-1912 May Myo 21 96.8 7.9 15 VII 198 22-06-1923 Matman 22.8 98.7 7.3 25 VII 432 08-08-1929 Swa 19.2 96.2 7.0 - - 254 12-03-1930 Pyu 18.2 96.3 7.5 10 VIII 360 05-05-1930 Bago 17.7 96.7 7.4 35 VI 428 27-01-1931 Kamaing 25.6 96.8 7.3 35 IX 531 16-08-1938 - 22.7 93.2 7.0 75 VII 242 12-09-1946 Tagaung 24.1 95.5 7.1 15 VII 333 16-07-1956 Sagaing 22 96 6.8 - - 147 08-07-1975 Bagan 21.8 94.7 6.5 157 V 72 29-05-1976 - 24.5 98.7 7.0 10 VIII 538 05-01-1991 - 23.6 95.5 7.0 19.7 VII 278 11-11-2012 - 23 95.9 6.8 13.7 VI 230 13-08-2016 - 23 94.8 6.9 136 VI 204 24-08-2016 Chauk 20.9 94.6 6.8 82 VI 40 * The epicentral distance is computed with respect to the closest of the four temples. Fig. 2. Location of the temples: Satellite view of the Bagan archaeological site and surroundings. D. Vecchio et al. Structures 74 (2025) 108584 4 two internal staircases in the east corners of the hall. The temple rises on three levels, including a square terrace with corner stupas, two additional square terraces without corner stupas connected by a staircase, and a square tower. Temple No. 844 is made of brick masonry with an average brick size of 33 cm ×17 cm ×4 cm. After the 1975 earthquake, the upper part of the square tower was destroyed (Fig. 3a). Furthermore, the small stupas at the corner of the first terrace suffered collapse. Due to the extensive damage, the temple underwent a massive reconstruction project, encompassing all collapsed parts (Fig. 3b). However, due to the 2016 Chauk earthquake, the temple again suffered severe damage, and it was tagged as “red”during the post-seismic survey. The most critical damage was the collapse of the central spire shown in Fig. 3b, which caused falling debris to destabilise smaller stupas and led to the overturning of pediments after impact [42]. An aerial photographic survey was taken during the post-earthquake assessment (see Fig. 3c). Due to the extent of the damage, Temple No. 844 was chosen as the reference for the numerical analyses in the present study. 2.3.2. Temple No. 558: Zan-Thi Temple No. 558 Zan-Thi dates to 1233 A.D. It is a small single-storey temple with a square central shrine (3.80 m ×3.83 m) (Fig. 4). The ceiling of the shrine is supported by a cloister vault, while the vestibules and the porches are topped with lower barrel vaults. The temple develops on three levels, featuring two square terraces and four corner stupas each, and a third terrace for the circular bell-shaped dome crowned on the top by a central conical spire. The masonry was built with bricks of average dimensions 30 cm ×15 cm ×4 cm. In the 1975 earthquake (Table 1), Zan-Thi suffered a shear failure in the conical spire, which was later partially reconstructed. Other reported types of damage were crown cracks in the north, west, and south entrances, specifically in the barrel vaults [42].Fig. 4 depicts two views of the temple, the first from 1992 and the second after the 2016 Chauk earthquake, respectively. The total collapse of the corner stupas and their subsequent reconstruction is evident. 2.3.3. Temple No. 1219: Kya-zin-hpaya Temple No. 1219 is a three-storey medium-sized temple with an estimated construction period of the 13th century A.D. (Fig. 5). The plan configuration is based on a rectangular central shrine. On the ground floor, the central shrine is topped by a barrel vault hipped at both ends. Two barrel vaults are also located above the entrance hall, the central corridor, the vestibules, and the porches. The temple develops on two inner levels, including an entresol level, and it features three external terraces. The lower level has a terrace with four corner stupas that is connected to the central terrace with an external staircase. The higher level shows a square terrace with corner stupas and a small temple. On top of the higher level, a square tower is crowned with a central conical spire. The bulk material of the building is made of masonry with bricks having average dimensions of 40 cm ×20 cm ×5 cm. According to the available findings [31], the temple was already repaired after the 1975 earthquake, suggesting that it suffered damages from past earthquakes. Fig. 5a is dated from 1930 and shows the extent of damage in the main spire, in the finial and in the small corner stupas. Subsequently, reconstruction work was conducted to rebuild the collapsed parts, as shown in Fig. 5b. However, according to the 2016 post-earthquake survey, the building suffered significant damage again in the upper-storey shrine, resulting in a “yellow”tag [42]. 2.3.4. Temple No. 1493: Myin-pya-gu Temple No. 1493 dates to the late 11th century A.D., the largest single-storey temple depicted herein (Fig. 6), with a solid core of 15.33 m ×15.55 m. The core is surrounded by a corridor with niches opening in the inner wall. Barrel vaults stand atop the shrines, corridors, Fig. 3. Temple No. 844 Tha-mu-ti-hpaya photographic documentation: (a) after the 1975 earthquake, (b) before the 2016 Chauk earthquake, (c) after the 2016 Chauk earthquake. (a) extracted from the ´ Ecole Française D ′ Extrˆ eme-Orient virtual photo library website, which allows free reproduction for study and research purposes. Fig. 4. Temple No. 558 Zan-thi East photographic documentation: (a) in 1992, (b) before the 2016 Chauk earthquake. (a) from [31]. Fig. 5. Temple No. 1219 Kya-zin-hpaya photographic documentation: (a) after the 1930 earthquake, (b) before the 2016 Chauk earthquake. (a) from [31]. D. Vecchio et al. Structures 74 (2025) 108584 5 vestibule, and porches. The three upper levels include three square terraces, of which the highest is topped by a circular bell-shaped dome and a central conical spire. The bricks have an average dimension of 36 cm ×18 cm ×5.5 cm. According to Pichard [31], the 1975 earthquakes damaged the temple at the central conical spire, which was reconstructed in 1978 (see Fig. 6). Despite the reconstruction, damage was again observed in the 2016 Chauk earthquake. For this reason, Temple No. 1493 was tagged as “yellow”in the post-earthquake survey [42]. 3. Assessment methodology 3.1. Modelling strategy The modelling strategy of this study aimed to assess seismic safety in Bagan temples No. 844, 558, 1219 and 1493, against overturning and shear-sliding. As is representative of the typologies of the region, the temples were modelled as assemblies of distinct continuous FEM multibodies that interacted with discontinuous contact interfaces, allowing separation, contact, and sliding phenomena [43]. Specifically referring to the lack of continuity among the different bodies and the presence of contact interfaces that connect the interacting surfaces, the position of the discontinuities with respect to the structural components were used to simulate failure planes where cracks are prone to form, allowing the separation of the temples into an assembly of multiple bodies. Due to the size of the temples and the relatively low stiffness of masonry, the flexibility of the structural components was considered explicitly, while less attention was given to local disintegration. Notably, following past seismic events, masonry disintegration was observed in some temples of the Bagan archaeological site due to the low quality of masonry. Such a failure mechanism is often detrimental to the integrity of the structure, as portions of masonry crumble into pieces [44]. Typically, an accurate evaluation of this failure mode requires very detailed and complex numerical models, where the masonry texture is simulated ad hoc together with the complex material constitutive laws of its constituents [17]. Given such complexity, the present study primarily investigated the dynamic rocking and sliding motion of the structural and non-structural components of the temples and, by selecting a texture-based model, attended less to the disintegration phenomena, as further described in Section 4. Damage observations of the past earthquakes showed that the most vulnerable components of the temples against overturning or sliding failures are located on the upper parts of the structure (i.e., above the roof terrace of the ground floor). Unlike the terraces, the ground floor did not reveal critical damage or the separation of the masonry leaves. Hence, the ground floor was not expected to rock, or slide but only to amplify the seismic signal. A simplified and convenient way to simulate the (elastic) dynamic properties of the ground floor was using independent single-degree-of-freedom (SDOF) oscillators. Instead of explicitly modelling the complete structure, only the vulnerable parts of each temple were considered. The dynamic behaviour of the ground floor was represented using three dynamically equivalent SDOF oscillators (two in horizontal directions and one in the vertical direction) to amplify the seismic signal at the interface with the modelled portion. This simplification reduced the number of DOFs of the numerical models: for example, for temple 844, total DOFs were reduced about four times (see Section 4.1). A validation of this simplified, yet efficient, modelling strategy was essential, and Section 4 corroborated this statement before the same was adopted in Section 5. 3.2. Numerical model A numerical model was constructed for each temple, including a system of three SDOF oscillators together with the discontinuous FEM model of the upper structure of the temples. The latter is composed of an assembly of distinct continuous FEM multi-bodies interacting among them only when in contact [43]. Next, each numerical model was subjected to nonlinear time history analysis (NLTHA) to assess the seismic vulnerability of the temples. The system of three SDOF oscillators idealised the elastic behaviour of the ground floor of the temples. The response of each SDOF oscillator was described using the equation of motion of an elastic oscillator. The equation of motion was expressed in a normalised form, using the frequency fand the damping ratio ξ[45]. The natural frequency fof each SDOF oscillator was selected following the proposal by Gavrilovic and Pichard [46]. These authors conducted in situ ambient vibration tests on 15 different temples of the Bagan archaeological site and measured the natural frequencies and damping ratios. Using regression analyses on the experimental natural frequencies, Gavrilovic and Pichard [46] proposed an equation to estimate the natural frequency of such temples by using only the plan dimensions and the height. Therefore, the horizontal natural frequencies fof the temples were estimated using simple geometrical considerations. Because the damping ratio ξis difficult to estimate from ambient vibration tests, it was assumed to be 5 %, a value within the range of the ambient vibration results. However, Gavrilovic and Pichard [46] provided no estimation of the vertical natural frequency, which was thus extracted from numerical simulations conducted by the Getty Conservation Institute [42].Table 2 collects the Fig. 6. Temple No. 1493 Myin-pya-gu photographic documentation: (a) after the 1975 earthquake, (b) after the reconstruction and after the 2016 Chauk earthquake. (a) from [31]. D. Vecchio et al. Structures 74 (2025) 108584 6 natural frequencies used for the system of SDOF oscillators for the four temples. The equation of motion was solved using the implicit numerical integration method of Newmark [45], and the response of each SDOF oscillator was imposed as an acceleration boundary condition at the base of the FEM model of each temple. The numerical model of the upper structure of each temple was developed using the commercial FEM software Abaqus/CAE [47]. Each model consisted of an assembly of elastic bodies and discontinuous interfaces which, under lateral action, may open and close, therefore allowing rocking or sliding. All analyses were run with the explicit solver of Abaqus/CAE, based on the central difference rule to integrate the equations of motion over time. Table 2 presents the main features of the numerical model for all the temples, including the average mesh size, the number of finite elements and the total number of nodes, while Fig. 7 shows the original (Fig. 7a–d) and the modelled geometry of the upper structure of the temples (Fig. 7e–h). Each ground floor was modelled using three SDOF oscillators to simulate its dynamic response and was omitted in Fig. 7e–h. Note that the original geometry was available from previous geometrical in situ surveys [16].Table 3 summarises the mechanical properties adopted for the homogenised masonry material and the contact interfaces. The masonry was assumed to have a density ρ =1620 kg/m 3 , homogenised elastic (or Young’s) modulus E=0.5 GPa, and Poisson’s ratio ν =0.2 [4,16,42]. An elastic modulus of 0.5 GPa was chosen according to previous in situ studies on other similar temples of the Bagan archaeological site [4,16,48]. In these studies, direct and indirect sonic tests were conducted on the Loka-Hteik-Pan temple, and the results showed an average value of 0.46 GPa. The study conducted by Bianchini et al. [4] also included in situ dynamic identification tests, and the outcomes were used to calibrate a numerical model, resulting in a Young’s modulus of 0.57 GPa. The contact interface properties were modelled with readily available constitutive laws of the Abaqus/CAE library. In the normal direction, the interfaces were modelled linearly with normal stiffness k n acting upon contact closure and no resistance upon separation. In the tangential direction, the contact properties followed the penalty formulation, defined by the tangential stiffness k s controlling the behaviour prior to sliding and the Mohr–Coulomb criterion using the friction coefficient μ [49]. Unlike for the masonry material properties, no experimental campaign for the evaluation of the interface properties of the temples was available. Accordingly, the contact interface properties were selected based on recommended values found in the literature. The friction coefficient was assigned the value of 0.75, as suggested in [50]. A preliminary investigation of temple No. 844 (omitted for brevity) showed that for a friction coefficient of 1.00, the response in terms of rotations and interface sliding of the structural elements varied on average by only 21 %, while the collapsed elements remained the same in all cases. The normal interface stiffness of contact was assigned a value of 100 MPa/m, in the lower range of values found in the literature [51,52,49] and references therein. A preliminary parametric investigation carried out by the authors was made using the “hard contact”model of ABAQUS/CAE, the upper bound of the normal interface stiffness. In this case, the response in terms of rotations and interface sliding of the structural elements varied on average by only 26 %, while the collapsed elements remained the same in all cases. In the tangential direction, the interface stiffness was assigned the value of 40 % of the normal interface stiffness, as suggested by Lourenço and Gaetani [50]. Regarding the energy loss at the contact interfaces, the Coulomb friction model captures the frictional hysteresis explicitly, acting only in the tangential direction of contact. Because during the rocking motion the displacements and velocities at the contact interfaces develop only in the normal direction, a viscous damping dashpot was assigned in this direction [52]. The adopted value of the damping ratio was ξ=2 %, also measured experimentally as a lower bound on interfaces [49]. Table 2 Natural frequencies and numerical features of the temples. Frequency Numerical properties Temple No. f E-W [Hz] f N-S [Hz] f V [Hz] Mesh size # elements # nodes 844 3.33 4.00 7.35 0.2 m / 0.4 m 13,338 16,280 588 4.20 4.20 6.35 0.2 m / 0.3 m 12,002 18,373 1219 3.33 3.84 9.89 0.2 m / 0.4 m 14,104 24,925 1493 3.33 3.33 9.10 0.4 m / 0.75 m 28,786 38,665 Fig. 7. Original geometries of the temples (a–d) and partitioned models of the upper structure (e–h): (a, e) Temple No. 844, (b, f) Temple No. 558, (c, g) Temple No. 1219, (d, h) Temple No. 1493. Table 3 Mechanical properties of the numerical models. Mechanical properties Value Units E - Young’s modulus 0.5 GPa ν –Poisson’s ratio 0.2 - ρ –Specific weight 1620 kg/m 3 k n –Contact normal stiffness 100 MPa/m k s –Contact tangential stiffness 40 MPa/m μ –Contact friction coefficient 0.75 - ξ – Contact viscous damping 2 % D. Vecchio et al. Structures 74 (2025) 108584 7 3.3. Seismic signal Eleven ground motion records were used as base inputs for the dynamic analyses. The suite included the real record of the 2016 Chauk earthquake, recorded at the Nyaung-U station, a simulation of the real record using a ground motion model, and nine additional simulations based on different rupture scenarios with a moment magnitude of M w 7.3 [53]. The nine simulated records were generated using a broadband physics-based seismological model combining different methods for the low and high-frequency ranges. The seismological model was initially benchmarked and tuned against the recorded 2016 Chauk earthquake to generate a “simulated Chauk”record. Next, the model was employed to generate hypothetical fault rupture scenarios of probable earthquakes. Basu and Raghukanth [53] explain that the Yenang-Chauk fault was selected because it is the closest known fault to the temples, and nine rupture initiation scenarios of M w 7.3 were simulated along its length. These rupture scenarios provided nine distinct seismic signals for each temple, in addition to the real and the simulated Chauk records. The seismological model was deterministic, and it did not correspond to a specific return period. The generated peak ground acceleration (PGA) (approximately 0.6 g) was close to a return period of 2475 years according to the probabilistic hazard model of G¨ ogen et al. [54]. Considering the above, the present work adopted 11 ground motions, being the most representative seismic scenarios available that could strike the temples under investigation. Fig. 8(a–l) shows the response spectra of the 11 ground motions for their three components: East–West (E–W) (Fig. 8a–d), North–South (N-S) (Fig. 8e–h) and Vertical (V) (Fig. 8i–l), while further details regarding the simulated ground motions can be found in Basu and Raghukanth [53]. Given the large number of analyses, the ground motions were labelled as follows: the recorded 2016 Chauk Fig. 8. Response spectra (a–l) and epicentral distances (m–p) for all temples in the E–W (a–d), N–S (e–h), and V (i–l) directions, including the recorded 2016 Chauk earthquake, the simulated Chauk and the nine M w 7.3 records. The dotted black vertical lines indicate the fundamental period of the SDOF oscillators in each direction. D. Vecchio et al. Structures 74 (2025) 108584 8 earthquake was abbreviated as “C”, while “SC”represented simulated Chauk, referring to the simulation of the 2016 Chauk record, and the nine rupture scenarios were denoted as “R#N”, where R represented the rupture and where #N =1–9. In addition, Fig. 8 indicates the adopted period of the SDOF oscillators of Table 2 to highlight the amplification caused by the ground floor. Fig. 8(m–p) shows the epicentral distances of the four temples from each rupture scenario. 3.4. 3D-motion and limit states This study represented the 3D orientation and motion of the structural components using intrinsic Euler angles to facilitate the discussion of the results of the numerical models. Such representation allowed a simple description of the 3D motion based on the principal axes of each body. Consequently, rocking was referred to as the case of rotation about the initial horizontal axes (i.e., x-x and z-z), and torsion was referred to as the rotation about the initial vertical axis (i.e., y-y). Note that this is not an elastic torsional distortion within the body, as often considered, but rather an integral rotation of the body around the vertical axis, also known as spinning. However, given both the complexity and the lack of interest in the torsional motion, only rocking and sliding were considered in this study. The Kabsch-Umeyama algorithm was adopted to compute the Euler angles [55]. For this purpose, the 3D displacements of at least three points of each structural component were recorded and used by the algorithm to estimate the rotation matrix over time. Next, the Euler angles were extracted from the rotation matrix using the XYZ notation, generating the so-called Tait–Bryan angles. While the rotations were calculated through the described algorithm, the sliding motion between two bodies was simply computed by subtracting the displacements of the two surfaces of each interface. To assess the vulnerability, the seismic demand of each temple was compared with its corresponding capacity. The comparison was made by adopting the instability angle α of each body as an estimator of the capacity against rocking and overturning. Given the Italian Code NTC 2018 [56], three limit states were adopted: (i) damage limit state θ/ α =0.4 (LS1); (ii) life safety limit state θ/ α =0.6 (LS2), and occurrence of the collapse when θ/ α >1 (LS3). The collapse caused by excessive sliding is also identifiable based on simple geometrical considerations, but the literature lacks a consensus on the definition of sliding limit states. The sliding collapse of each pair of bodies was assumed to occur when the relative displacement exceeded the width of a body that served as a base for another body placed on top of it (LS3). 4. Validation of the modelling strategy The proposed modelling strategy adopted two marked simplifications: i) the ground floor substitution with three SDOF oscillators for each direction; and ii) the neglection of masonry disintegration. Both simplifications, and particularly their combination, reduced and made viable the computational cost of the numerical analyses. For example, for the medium size temple No. 844, each of the simplifications alone reduced about four times the number of DOFs of the FEM model, and their combined influence resulted in a model with about 16 times fewer DOFs. However, the implications of the structural response had to be examined and validated. In this section, these simplifications are discussed based on the comparison between the simplified and the detailed models for the representative temple No. 844 subjected to the 2016 Chauk earthquake. 4.1. The SDOF oscillator simplification The validation of the simplification using SDOF oscillators was performed by comparing it with a complete model of the temple, in which the lower part of the temple was modelled as linear elastic. Fig. 9a illustrates the complete model of the temple No. 844 together with its mesh discretisation. This model was composed of 49,388 elements and 61,565 nodes, whereas the simplified model presented in Section 5 included 13,338 elements and 16,280 nodes. The comparison of the two models was performed by analysing the signal amplification of the ground acceleration at the roof terrace height of the ground floor. In the complete model, acceleration time histories were extracted from two nodes, one at the centre (P 1 ) and one at a corner edge (P 2 ), as shown in Fig. 9a. The model of the SDOF oscillators provided a single acceleration time history for each direction. The bar plot in Fig. 8b shows the peak acceleration values of the nodes in each direction for both models, including the ground motion. First, it is possible to observe the amplification of the signal at higher levels of the structure in comparison with the ground motion. Moreover, the simulation of the complete temple showed that the corner node was subjected to lower peak acceleration with respect to the centre node, while the response of the SDOF oscillator showed a peak value between the two peaks of the complete model. The results showed that the SDOF oscillator (orange bar) experienced an acceleration of 0.23 g and 0.28 g for the E–W and N–S directions, respectively, while in the complete model, P 1 reached peaks of 0.28 g and 0.33 g, and P 2 reached peaks of 0.15 g and 0.16 g. However, the V direction was characterised by a higher scatter. The corner node was subjected to an acceleration of 0.15 g, and the centre node to 0.72 g in the complete model, while the response of the SDOF oscillator showed an acceleration of 0.46 g. Overall, the two nodes of the complete model showed notable differences despite being at the same height. The discrepancy was due to the spatial configurations of the model and their associated mass and Fig. 9. (a) Complete model of temple No. 844 with the location of the nodes P 1 and P 2 for comparison, and (b) amplification of the signal and comparison between simplified and complete model. D. Vecchio et al. Structures 74 (2025) 108584 9 stiffness. The corner node P 2 stands on the solid substructure of the ground floor, while the central node P 1 stands on a hollow vaulted corridor. As a result, central node P 1 showed lower stiffness in the vertical direction when compared to corner node P 2 . Additionally, the masses standing atop the corner node P 2 are the small stupas, while the central node P 1 was placed below the central spire of the temple. These differences resulted in dissimilar dynamic properties, particularly in the vertical direction, and central node P 1 experienced a lower frequency of vibration when compared to corner node P 2 . Based on these considerations, the intermediate acceleration amplification given by the simplification of the SDOF oscillator was considered acceptable for the purpose of this engineering application, which involved numerous simulations. It is also widely accepted that, as compared to the two horizontal ones, the vertical component has a relatively small influence on the rocking and sliding motion [57,58]. 4.2. Meso-modelling Masonry disintegration occurs due to the poor quality of masonry, which cannot resist horizontal actions [44]. Especially when dealing with rubble masonry often characterised by multiple leaves and poorly connected stones (e.g., insufficient overlap of masonry units), such failure is detrimental [22]. Numerical models can simulate such a response only when the texture of the blocks is explicitly modelled, namely by using the so-called micro-modelling or meso-modelling approaches. Considering the above, a meso-model was constructed for the upper part of temple No. 844, and the results of an additional analysis were compared with the numerical reference model in Fig. 7. Given the lack of in situ information, the meso-model was characterised by 3D elastic blocks and contact interfaces, which represented the masonry pattern of temple No. 844 in a simplified fashion. The mechanical properties of the blocks and the interfaces are reported in Table 4. The normal and tangential interface stiffness differed from those of Table 3 due to the transition from the macroto the mesoscale. The initial values in Table 3 were representative of the homogenised global behaviour of the masonry. However, in the meso-model, additional flexible interfaces at the brick-and-mortar scale were introduced to explicitly represent smaller portions of the masonry. To ensure equivalence in stiffness between the macroand meso-models, the equation proposed by Lourenço [59] was employed. This equation relates the normal interface stiffness of the joints k n,j with the height of the brick h, the elastic modulus of the wall E wall and the brick E brick , as follows: kn,j=1 (h(1 Ewall −1 Ebrick)) (1) The friction coefficient was decreased to a lower-bound value to allow disintegration, while other failure modes (e.g., splitting or crushing) were disregarded. The updated mechanical properties of the meso-model are shown in Table 4. Fig. 10a depicts the geometry of the meso-model. Here, the central spire of the temple was modelled with 16 horizontal courses, while the square tower was split into 16 outer courses and corresponding inner infills. In total, the meso-model (Fig. 10a) was composed of 25,669 elements and 61,989 nodes, so the size ratio with the model presented in Fig. 7 was approximately 4:1 (Table 2). Evidently, this approach led to a significant increase in the computational effort. Fig. 10b compares the results of the two models in terms of the displacement profiles up the height in the two directions when subjected to the 2016 Chauk earthquake. The left part of the plot in Fig. 10b compares the displacement profiles for the N–S direction, while the right part of the graph compares the E–W displacement profiles, respectively. The red solid lines represent the displacement profile in the macro model, computed at the different heights of the interfaces noted with square markers (also indicated in Fig. 10a). Instead, the black dotted lines show the displacement profile of the meso-model, extracted at all the levels of the different block courses noted with cross-markers (the courses shown in Fig. 10a). Fig. 10b illustrates that the displacement profiles of both models increased gradually up the height of the structure, starting from approximately 5 cm displacement at the elevation of 8.5 m and increasing on average up to 70 cm at the top of the spire. Moreover, Fig. 10b illustrates that the increase of the displacement profiles was steeper at higher elevations, due to the larger rotations experienced by the upper parts of the structure. In general, both directions had a similar trend, where the displacement profiles of the two models differed slightly for the square tower, while a notable disagreement was apparent for the central spire. The meso-model showed less displacement at lower heights than the macro-model, while the contrary occurred at the higher courses of the tower and the central spire. Such an outcome is often the case for dynamic analysis of slender structures due to the effect of higher modes [60]. The last course of the spire of the meso-model experienced collapse, which was not found in the macro-model, highlighting the importance of local disintegration. This additional outcome may justify the partial collapse found in the 2016 Chauk earthquake. The overall comparison of Fig. 10b illustrates that the two models had a reasonable agreement in terms of displacement profiles along the height, allowing the use of the more computationally efficient macromodel to study the rocking and sliding response of the temples for preliminary simplified analyses, or in case measures are taken to ensure that disintegration is unlikely (e.g., repointing of mortar joints, masonry grout injection, or application of a reinforced render). This validation justifies subsequently adopting this approach, omitting any attempt to compare the results with the damage experienced in the 2016 Chauk earthquake, as it was assumed that measures to avoid disintegration remained necessary. 5. Results 5.1. Temple No. 844 Fig. 11 collects the numerical results for temple No. 844, subjected to the 11 ground motion records, following the colour notation shown at the bottom of the figure. The results are given for the main structural components in the elevation of the temple (i.e., the central spire and its base), and one representative corner stupa with its corresponding base. Fig. 11a shows the maximum rocking responses θ xx and θ zz of each analysis for both horizontal axes (x-x and z-z) and normalised by the instability angle α of each structural component. Fig. 11a also highlights the three limit states (i.e., LS1, LS2, and LS3) of the rocking motion with vertical dashed lines, as introduced in Section 4. The total collapse is shown by the points crossing the LS3 vertical line. Similarly, Fig. 11b presents the sliding response at the interfaces between the structural components, where only the ultimate limit state LS3, corresponding to collapse, is depicted. Fig. 11a shows that temple No. 844 experienced minor rocking motion, except for the corner stupa, which overturned due to the ground motion R4. Aside from the collapse case, none of the structural components reached the LS1 limit state, with the central and the corner stupas, together with their bases, undergoing small rotations up to θ/ α =0.1 for all the ground motion scenarios. The stockier base of the central spire Table 4 Mechanical properties of the numerical meso-model. Mechanical properties Value Units E - Young’s modulus 1.5 GPa ν –Poisson’s ratio 0.2 - ρ –Specific weight 1620 kg/m 3 k n –Contact normal stiffness 1900 MPa/m k s –Contact tangential stiffness 800 MPa/m μ –Contact friction coefficient 0.4 - ξ – Contact viscous damping 2 % D. Vecchio et al.