The power of seismic performance assessment of masonry monuments based on simulated ground motions: A case study

The power of seismic performance assessment of masonry monuments based on simulated ground motions: A case study Ahmad Fathi * , Shaghayegh Karimzadeh , Paulo B. Lourenço University of Minho, ISISE, ARISE, Department of Civil Engineering, Guimar˜ aes, Portugal ARTICLE INFO Keywords: Arge-Tabriz Ancient masonry Cultural heritage buildings Fragility curves Stochastic finite-fault simulations Simulated records Rupture uncertainties ABSTRACT A comprehensive investigation of the seismic behaviour of the Arge-Tabriz monument, located in a region with very high seismicity in northwest Iran, is presented in this study, making use of simulated ground motions. In seismically active regions characterised by scarcity of large-magnitude events owing to a seismic gap, ground motion simulations provide an effective tool to assess potential scenario earthquakes. In this study, simulations based on the stochastic finite-fault simulation method, incorporating the dynamic corner frequency concept, are employed for scenarios with varying magnitudes. These simulations account for the rupture uncertainties of the North Tabriz Fault in northwest Iran, a region known for its large seismic hazard. The earthquake records are systematically chosen from an extensive dataset covering a broad range of peak ground accelerations. The structure is modelled using the finite element method in ABAQUS. After calibration of the numerical model in terms of the building’s natural frequencies, based on experimental measurements, a non-linear static analysis is performed. Subsequently, the capacity curve of the building is drawn from the pushover analysis using the N2 method, followed by dynamic time history analyses using the simulated ground motion records. Finally, fragility curves are derived based on the drift and Park-Ang damage indices. Thresholds of the Park-Ang damage index for different performance limit states are adopted for masonry. Crack patterns of the structure resulting from numerical modelling are in good agreement with the current situation of the building. The fragility curves, together with the crack patterns, confirm that the monument is able to withstand significant seismic events. The results show that the monument exhibits an extremely high probability of exceeding the damage limitation state, close to 100 %, at a PGA of 0.35 g. For the significant damage limit state, this probability is about 20 % at 0.35 g and increases to 40 % at 0.6 g. Notably, the near collapse limit state was not induced even under the maximum considered PGA of 0.6 g, suggesting a robust ultimate capacity. The findings of this study illustrate the power of using simulated ground motions for exceptional buildings and will guide future conservation efforts for this monument. 1. Introduction Masonry construction possibly appeared by making simple dry-stone walls without mortar, then developed with stone and mud bricks bonded with earthen mortar [1]. The advances in technology and materials, such as in cutting and carving stone, firing bricks, and lime and pozzolanic additions to make mortar, resulted in complex masonry structures, such as massive walls, palaces, and churches over the world throughout history. Unreinforced masonry (URM) buildings form the highest percentage of the built heritage worldwide [2]. Natural hazards such as earthquakes and floods often affect these buildings. More specifically, despite their proper resistance against vertical loadings, historic masonry buildings have not been designed to withstand horizontal loads and are susceptible to seismic loads [3–5]. This vulnerability is mainly due to the low tensile strength of masonry, inadequate connections between structural elements, poor capacity to dissipate energy, low out-of-plane capacity, irregularities in plan/elevation, decay of the material over time and large inertia forces [3,6–8]. Conducting a seismic assessment of masonry buildings is a key step for taking actions to minimise their vulnerability to earthquakes and to preserve them. Various modelling strategies have been proposed in the literature to model masonry structures, including block-based models, continuum models, macro-element models, and geometry-based models. * Correspondence to: ISISE, Department of Civil Engineering, Campus de Azur´ em da Universidade do Minho, Guimar˜ aes 4800-058, Portugal E-mail address: [email protected] (A. Fathi). Contents lists available at ScienceDirect Structures journal homepage: www.elsevier.com/locate/structures https://doi.org/10.1016/j.istruc.2025.110263 Received 14 February 2025; Received in revised form 14 September 2025; Accepted 18 September 2025 Structures 81 (2025) 110263 Available online 23 September 2025 2352-0124/© 2025 Institution of Structural Engineers. Published by Elsevier Ltd. All rights are reserved, including those for text and data mining, AI training, and similar technologies. Different strategies can be adopted depending on the accuracy and simplicity required [4,9]. However, micro-modelling and macro-modelling, which are block-based and continuum models, respectively, are widely used in engineering practice to study the structural behaviour of the masonry members and buildings [10–14]. Micro-models replicate masonry components, i.e., units and joints, as continuum elements, while interface elements with zero thickness represent the interface between these two components (blocks). The approach can predict all possible failure mechanisms of masonry [1] but it is mostly limited to the modelling of structural members, since it requires a high level of detail and huge computational demands [4,9]. Conversely, macro-models are not able to reproduce the detailed geometry of a masonry element, because they treat masonry as a homogenous and continuous material. This approach can only capture the overall structural behaviour of masonry as a composite material. Yet, macro-models are appropriate to study the behaviour of the entire structure since they are faster and easier to implement. A variety of applications of the approach in studying the behaviour of full-scale masonry buildings under static and dynamic loading can be found in the literature [4,5,15–25]. The seismic assessment of masonry buildings presents an engineering challenge due to the wide range of variables affecting their structural behaviour and the inherent uncertainties that rise with the intensity of earthquakes [26]. Using fragility curves is an efficient strategy for a probabilistic characterisation of the seismic response of masonry buildings [27,28]. Fragility curves are analytical methods for evaluating the extent of damage that a structure may endure during seismic loadings. These curves represent the relationship between seismic intensity and damage in terms of the conditional cumulative probability of reaching or exceeding a predetermined damage state, corresponding to a specific design value of the chosen demand parameter [29,30]. Fragility curves are common for assessing the seismic vulnerability of masonry structures. Numerous examples of applying this method to masonry buildings can be found in the literature, further highlighting its relevance and effectiveness [28,31–37]. In the development of analytical fragility curves, the use of real ground motion records from earthquakes recorded worldwide is a common practice for conducting structural analysis [38–40]. However, the availability of such records is often limited, particularly for large-magnitude events and high-intensity levels that accurately represent the specific seismological characteristics of a particular region [41]. To overcome this limitation, ground motion simulations offer an alternative approach [42]. By incorporating pertinent seismological, tectonic, geological, and geophysical information, ground motion simulations can generate region-specific ground motion records, facilitating the customisation of ground motions for specific case studies to increase their adequacy. These simulations also provide records with higher intensity levels, eliminating the need for scaling or modifying real ground motions obtained from global earthquake databases. Consequently, recent literature has witnessed the successful utilisation and validation of simulated ground motion datasets in various engineering applications [43–49]. Within this context, Otarola et al. [50] employed ground motion simulations in the seismic performance assessment of bridge structures. Smerzini et al. [51] explored the use of physics-based numerical simulations to generate earthquake ground motions, particularly in regions with limited seismic data, such as near-source areas and complex geological settings. Their study introduced a validated dataset of broadband near-source earthquake ground motions, demonstrating that their use in inelastic time-history analyses did not introduce systematic biases in engineering demand parameters compared to those derived from recorded motions. Rosti et al. [52] validated the use of 3D physics-based numerical simulations for generating ground shaking scenarios in seismic fragility studies, using the 2009 L’Aquila event as a case study, demonstrating the effectiveness of simulations in region-specific seismic vulnerability and risk assessments by comparing observed and predicted damage distributions. In the context of masonry structures, Karimzadeh et al. [33,53] utilised simulations to develop fragility curves for simplified structural models, such as equivalent single-degree-of-freedom models, specifically tailored to regions in Türkiye. More recent efforts by Bernardo et al. [54] and Karimzadeh et al. [55] have concentrated on masonry cultural heritage and establishing validation frameworks for using stochastic simulations in the seismic response evaluation and fragility curve development of structures in the Azores region, Portugal, which experiences moderate seismicity. Building on these advancements, this research investigates the seismic vulnerability of the ancient masonry building Arge-Tabriz in northwest Iran under simulated site-specific ground motions. A recent study by Hoveidae et al. [20] valuated the seismic response of Arge-Tabriz for four different scenario events due to rupture of the North Tabriz Fault (NTF) simulated using the stochastic finite-fault approach by Motazedian and Atkinson [56]. However, their study did not address record-to-record variability or generate fragility curves, key aspects often overlooked. This gap is particularly significant for historical monuments like Arge-Tabriz, which require detailed structural modelling and advanced seismic performance evaluation. In this context, the present study introduces a framework to overcome these challenges by integrating simulated site-specific ground motions, thereby addressing the scarcity of real strong motion records in the region. Unlike traditional approaches that rely on scaling existing records, the proposed method captures site-specific seismic hazards more reliably. Furthermore, fragility analysis is employed, with the major innovation being the development of fragility curves for ArgeTabriz and the proposal of new thresholds for the Park-Ang damage index tailored to historical masonry structures, offering a more reliable basis for seismic performance evaluation of historical structures. A recent study by Temiz et al. [57] developed a region-specific ground motion dataset for the Tabriz region, which was generated using the stochastic finite-fault ground motion simulation method by Motazedian and Atkinson [58]. In the present paper, the methodology for ground motion simulation and the adopted input parameters are detailed in Section ↱2 while the selection of ground motion records is discussed in Section ↱3. Section ↱4 describes the case study of the Arge-Tabriz monument and the finite element model. Section ↱5 and ↱6 discuss the nonlinear static and dynamic analyses, respectively, leading to the derivation of fragility curves. Finally, Section ↱7 summarises the key findings and conclusions, highlighting the implications for the seismic assessment of historical masonry structures. 2. Simulated earthquake dataset Alternative simulation methodologies offer a comprehensive timeseries of potential earthquake scenarios, encompassing a wide range of magnitudes, varied site characteristics, and different source-to-site distances. These methodologies include deterministic, stochastic, and hybrid approaches [42]. Deterministic methods rely on well-defined seismic sources and velocity models as inputs to numerical solutions for the partial differential equations associated with wave propagation [59]. These methods are typically suitable for the lower frequency range of seismic records due to wavelength constraints. On the other hand, stochastic approaches primarily simulate higher frequencies that exhibit incoherent and random characteristics. By combining random phases with a deterministic far-field shear wave spectrum, stochastic techniques generate the mean horizontal component of the records [60]. Stochastic simulation methods can be further categorised into source-based and site-based approaches [42]. Source-based methods explicitly model the source, path, and site effects [56,61], whereas site-based methods implicitly account for these effects through calibrated empirical ground motion models [62,63]. Both approaches have been widely employed in state-of-the-art research and engineering A. Fathi et al. Structures 81 (2025) 110263 2 practices such as [20,64]. Hybrid ground motion simulation approaches are preferred to simulate the complete frequency band accurately. These techniques employ deterministic methods for simulating the low-frequency portion of the records and stochastic methods for simulating the high-frequency portion [65]. In the subsequent sections, the ground motion simulation methodology, including the input model parameters, will be discussed. Subsequently, the results of the simulations will be presented, followed by a description of the selection process for the records in the final stage. 2.1. Methodology In regions lacking comprehensive regional source and velocity models, the stochastic method emerges as a practical and preferred approach for ground motion simulation. This methodology provides significant advantages by effectively addressing the broadband frequency range, which is of particular importance for many structural types. The stochastic method primarily focuses on simulating the incoherent mediumto high-frequency components of ground motions, typically exceeding ~1 Hz, even in the absence of a complex source model and detailed wave propagation information [56,60,66,67]. It is noted that compared with traditional scaling techniques, which adjust limited worldwide real records to match target spectra, the stochastic simulation approach employed in this research provides a physically consistent suite of motions that inherently captures record-to-record variability without introducing spectral distortions. Despite their limitations, stochastic techniques exhibit the ability to accurately simulate ground motion amplitudes within the frequency range of interest to engineers. Numerous studies have demonstrated their effectiveness in capturing the characteristics of real ground motion records [68–74]. By incorporating the stochastic approach, engineers can appropriately address the uncertainty and variability associated with ground motion, even in regions with limited available data and complex geological conditions. The simulation of scenario events involves the utilisation of the EXSIM-beta platform, which is built upon the stochastic finite-fault ground motion simulation technique developed by Motazedian and Atkinson [56]. This method employs a rectangular representation of the fault plane and treats it as a combination of smaller sub-faults functioning as individual point-sources [60]. By considering source, path, and site effects, the acceleration response spectrum is modelled in the frequency domain for each of these point-sources. Subsequently, the responses of all sub-faults are aggregated in the time domain to obtain the final response of the primary fault in the following manner: a(t) = ∑nw i=1∑ nl j=1 aij(t+Δtij)(1) where, a(t)is the acceleration of the entire rectangular fault at time t and aij corresponds to the acceleration time-series of the ijth sub-fault at time t+Δtij. Here, Δtij is the relative time delay of the ijth sub-fault to the observation point. Finally, the terms nw and nl are, respectively, the considered number of sub-faults along the width and length of the rectangular fault plane. The acceleration spectrum is calculated for every individual identical sub-fault as follows: Aij(f) = Rθφ × 2 √ 4 πρ β2 M0Sij ∑nw k=1∑ nl l=1 Skl Hij (2 π f)2 1+(f f0ij )2e− π fRij Q(f)βG(Rij)e− π KfA(f) (2) where Aij(f)is the acceleration spectrum of shear wave of the ijth subfault, Rθφ represents the radiation pattern, ρ and β are, respectively, the density and the shear-wave velocity, M0 is the total seismic moment of the rectangular fault plane and Sij is the relative slip weight of the ijth sub-fault. The term Rij is the distance of the ijth sub-fault from the observation point, Q(f)is the quality factor, G(Rij)is the geometric spreading factor, e− π Kf is a high-cut filter included to provide the spectral decay at high frequencies described with the Kappa factor (K) of soils [75], and A(f)is the site amplification. The term f0ij represents the dynamic corner frequency of the ijth sub-fault, which is calculated as: f0ij(t) = NR(t)−1/3×4.9×106β(Δ σ M0/N)1/3(3) where, Δ σ is the stress drop, NR(t)is the cumulative number of ruptured sub-faults at time t, and N is the total number of sub-faults (nw×nl). Finally, Hij is a scaling factor introduced to conserve the highfrequency spectral level of the sub-faults, which reads: Hij = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ N∑⎡ ⎢ ⎢ ⎢ ⎣ f2 1+(f f0)2⎤ ⎥ ⎥ ⎥ ⎦ ∑⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ f2 1+(f f0ij)2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ 1/2 (4) where N is the total number of sub-faults, and f0 is the corner frequency at the end of the record, which is given by: f0=f0ij(t=tend) = N−1/3×4.9×106β(Δ σ M0/N)1/3)(5) 2.2. Input-model parameters Tabriz, a historically significant city in Iran, is renowned for its cultural heritage and is home to notable monuments like Arge-Tabriz. Geographically, Tabriz is situated between eastern Anatolia, the western Caspian Sea, the southern Caucasus thrust belt, and the northern Zagros Mountain range. The city is characterised by the North Tabriz Fault (NTF) [58], which has experienced several seismic activities throughout history, including notable events in the 18th century [76]. However, detailed and reliable scientific data from those periods are lacking, leading to an incomplete understanding of the extent of destruction and ground deformations caused by historical earthquakes. Nevertheless, more recent seismic events in the 20th century, though less destructive, highlight the ongoing seismicity of the region. Given the active nature of the NTF and the seismic gap between large-magnitude earthquakes, city centre of Tabriz, which includes the Arge-Tabriz monument, has been selected as the study area, aiming to address the lack of recorded motion data for potential seismic events. Fig. 1 presents the tectonic map of northwest Iran with the historical earthquakes. A recent study by Temiz et al. [57] conducted ground motion simulations for the Tabriz region at randomly distributed nodes associated with the potential ruptures of NTF. The resulting dataset was validated using a robust and multi-step framework, including direct comparisons with regional ground motion models, evaluation of attenuation trends, and assessment of inter-period correlation to ensure the physical reliability and statistical consistency of the simulations. To generate a comprehensive dataset of ground motion records, extensive simulations have been conducted at these random nodes, encompassing scenario earthquakes of varying magnitudes (M w =7.7, 7.4, 7.1, and 6.8). The estimation of ruptured fault length and width is based on Wells and Coppersmith’s study [78]. The NTF, with an overall length of approximately 240 km [77,79], is divided into multiple alternative ruptured fault lines to account for uncertainty. For instance, the largest event (M w A. Fathi et al. Structures 81 (2025) 110263 3 =7.7) is divided into 17 alternative fault lines, each spanning 140 km. Similar divisions are made for other magnitudes, resulting in multiple scenarios. According to Temiz et al. [57], uncertainty in the hypocentre location is considered, treating the focal depth and epicentre location as random variables. Previous research suggests a median focal depth of 12.1 km with a standard deviation of 4.9 km for earthquakes in the region [80]. The focal depth is considered within a range of 6.0–18.0 km, while the epicentre location is assumed to vary along the fault’s length. Other input-model parameters are determined deterministically based on previous studies including strike [79] and dip angles [81], quality factor [82], Kappa factor [83], site amplifications [84], stress drop [85], in addition to geometrical spreading, pulsing percent, rupture velocity, density, and window type [20,86]. A summary of this information is provided in Table 1. This study utilises the input model parameters provided by Temiz et al. [57] to conduct simulations at the Arg-Tabriz site, incorporating the V S30 information for the location to deliver site-specific results. It is noted that this study considers site amplifications for a soil class characterised by a generic shear-wave velocity, as reported by Karimzadeh et al. [87]. The soil amplification factors are derived from the work of Boore and Joyner (1997) [88] using quarter wavelength approach developed for an average soil class (V S30 =310 m/s). 3. Results of simulations Through simulations conducted at the designated site, a comprehensive dataset consisting of 4134 acceleration time-series has been generated. This dataset covers a wide range of ground motion intensity levels resulting from the rupture of the NTF while considering the uncertainties associated with the rupturing process. Fig. 2a visually represents the distribution of the generated ground motion records for different scenario events using a pie chart, highlighting the relative proportions. Furthermore, Fig. 2b illustrates the relationship between the magnitude (M w ) and the Joyner and Boore distance (R JB ) through a scatter plot. It is evident that a greater number of simulations are conducted for the largest event, owing to the consideration of additional uncertainties and the inclusion of a larger number of records in the overall dataset, as this event is deemed more critical. Moreover, the dataset encompasses a variety of distances spanning roughly 10–82 km resulting from the initiation of rupturing fault planes from distinct points. The obtained records are analysed in terms of peak ground motion parameters, specifically peak ground acceleration (PGA) and peak ground velocity (PGV), yielding the following ranges for the scenario events for the events with Mw=7.7, 7.4, 7.1, and 6.8: a) PGA ranges 0.19–0.88 g, 0.09–0.59 g, 0.03–0.55 g, and 0.01–0.42 g, respectively; b) PGV ranges 14.0–62.1 cm/s, 6.1–50.6 cm/s, 2.3–35.5 cm/s, and 1.6–33.5 cm/s, respectively. In addition, Fig. 3 illustrates the distribution of PGA for different magnitudes. As the event magnitude increases, the probability distributions tend to exhibit lognormal characteristics, likely due to the extensive number of simulations conducted for this event. It is noteworthy that most of the ground motions at the site of interest have a maximum PGA not greater than 0.6 g. However, a smaller portion of the data does exhibit higher PGA values, reaching up to 1.0 g. 4. Ground motion record selection In order to provide proper input ground motions for structural Fig. 1. Tectonic map of northwest Iran with the historical earthquakes. The red star shows the location of the Arge-Tabriz monument. The figure is adapted from Azad et al. [77]. Table 1 Deterministic input-model parameters used for ground motion simulations [57]. Parameter Value M w 7.7 7.4 7.1 6.8 Length (km) 140 84 52 32 Width (km) 24 20 16 12 Stress drop (MPa) 11.0 9.5 8.0 6.5 Fault mechanism Strike-slip Strike (◦) 115 Dip (◦) 90 Shear wave velocity (m/s) 3.3 Rupture velocity/Shear wave velocity 0.8 Density (kg/m 3 ) 2800 Pulsing percent (%) 35 Quality factor 103f0.88 Geometrical spreading*R−1forR ≤85 R−0.5forR >85 Kappa value 0.035 Window type Saragoni-Hart Site amplifications Generic amplification factors for average soil [88] Sampling time (s) 0.005 * R: Source to site distance A. Fathi et al. Structures 81 (2025) 110263 4 analysis, a total of 60 records have been carefully selected to encompass a broad range of PGAs, which is used as the chosen intensity measure for deriving fragility curves in this study. The selected records span a PGA range up to 0.6 g, divided into six intervals of 0.1 g each. Within each interval, 10 records are included to capture a diverse range of seismological parameters, including uncertainties associated with ground motion records. Due to the high number of ground motions, it is not convenient to include the time histories of all records, however, a detailed list of the selected records can be seen in Table 2 and the time histories of six records (one from each interval) are indicated in Fig. 4 as Fig. 2. A graphical representation of:(a) the distribution of the ground motion dataset with respect to magnitude (M w ); and b) the distribution of M w with respect to the Joyner and Boore distance (R JB ). Fig. 3. The distribution of the peak ground acceleration (PGA) with respect to different magnitudes for the entire ground motion dataset. A. Fathi et al. Structures 81 (2025) 110263 5 examples. Each graph represents a unique seismic event, with varying amplitude, frequency content, and duration, which can provide insights into the intensity and characteristics of the ground shaking during the event. Fig. 5 illustrates the distribution of PGA versus PGV for the selected records. The records exhibit a homogeneous scattering pattern, effectively representing the full PGA and PGV ranges reaching up to 0.6 g and 70 cm/s, respectively. 5. Case study: Arge-Tabriz monument 5.1. Numerical modelling Based on the historical documentation, Arge-Tabriz, which dates back to the 14th century, was damaged due to severe earthquakes in 1721 and 1780 [20,89]. However, since the developed cracks and damages were repaired, as seen in Fig. 6, a non-damaged model of the building is adopted here to perform the seismic behaviour assessment of this monument. This simplification was assumed to minimise the computational time and to avoid the complexity of modelling cracks and details that weakly known. This would also introduce a complicated mesh with small size of the elements, which likely leads the analysis to convergency issues. However, to minimise the impact of this simplification on the final results and assure the reliability of the model, the same was calibrated based on the experimental measurements of its natural frequencies [86], as detailed in the subsection “Model Calibration”. This calibration, which achieved a maximum error of 2 %, implicitly accounts for the current structural stiffness and dynamic characteristics, thereby reflecting existing damage that may influence the global dynamic response. As a result, despite the initial simplification, the calibrated model provides a robust and reliable basis for the subsequent seismic assessment. After constructing a 3D model of the building according to the laser scans provided by the cultural heritage organisation (CHO) of East Azerbaijan province, the meshing in ABAQUS ended up with 25,864 solid elements, type C3D4, and 129,367 nodes (Fig. 7). Detailed drawings of Arge-Tabriz can be found in [20]. As stated, the structure is modelled using the macro-modelling approach, considering masonry as a homogeneous material. This method of modelling is one of the most common strategies to analyse the behaviour of the masonry structures [18,19,21,90–93], being able to replicate the main failure mechanisms of the masonry and requiring low computational time. This is important in case of large structures with a high number of elements and nodes, compared to detailed strategies such as micro-modelling [94]. 5.2. Model calibration In order to minimise the influence of the simplifications on the final results of the analyses and to validate the material properties used, i.e., to ensure the overall accuracy of the modelling process, the structural model needs to be calibrated using the real case study. In this research, calibration was performed based on the available experimental data, measurements of the natural frequencies of the structure for the first six vibrational modes, which have been obtained from the microtremor measurements on the Arge-Tabriz monument by Fallahi et al. [95]. The details of the measurements and methodology to extract the natural frequencies and shear wave velocity can be found in [95,96]. Due to the lack of experimental data on the mechanical properties of the material, reference values belonging to the nearest (in terms of distance, material type, and construction era) ancient structure, namely the Tabriz historic Bazaar, have been initially adopted from Aghabeigi et al. [97]. The study includes experimental assessments of the elastic modulus, compressive strength, and Poisson’s ratio (based on three compressive tests), as well as tensile strength determined through three-point bending tests on three flexural samples loaded parallel to the bed joints. Then, modal analysis has been performed to extract the frequencies of the model and to compare them with the experimental results [95]. Note that the corresponding experimental modes are not available, but the large number of frequencies available provides good confidence in the numerical model update. The process is carried out in the elastic range, so that the E-modulus (which was 2060 MPa initially, according to [97]) and density of the material were changed gradually until the numerical results of the first six modes were in agreement with the experimental data. Table 3 represents the comparison of the Table 2 List of the selected ground motions for Tabriz. Record ID Magnitude (M w ) R JB [km] PGA [g] P Sa * [g] PGV [cm/ s] D 5–95 [s] PGA Interval [g] 15–1–1 6.8 77.9 0.01 0.03 2.0 18 0.0–0.1 1–16–1 6.8 66.2 0.02 0.04 2.1 18 0.0–0.1 14–4–1 6.8 65.6 0.02 0.03 2.1 12 0.0–0.1 13–16–1 7.1 57.4 0.03 0.06 2.3 26 0.0–0.1 12–9–1 6.8 41.1 0.04 0.07 3.2 19 0.0–0.1 13–19–1 7.1 57.4 0.05 0.09 3.4 17 0.0–0.1 1–8–1 7.1 47.5 0.06 0.15 5.3 15 0.0–0.1 13–37–1 7.1 57.4 0.07 0.13 4.1 4 0.0–0.1 11–14–1 7.1 29.3 0.08 0.12 6.2 26 0.0–0.1 11–10–1 6.8 29.3 0.08 0.15 4.6 15 0.0–0.1 3–39–1 7.1 25.7 0.10 0.17 8.3 25 0.1–0.2 12–5–1 7.4 29.3 0.12 0.20 8.1 28 0.1–0.2 11–13–1 7.1 29.3 0.12 0.29 13.9 3 0.1–0.2 10–11–1 6.8 18.3 0.13 0.15 13.1 11 0.1–0.2 10–2–1 7.1 18.3 0.14 0.15 13.5 23 0.1–0.2 1–33–1 7.4 23.1 0.15 0.39 10.9 12 0.1–0.2 5–8–1 6.8 20.7 0.15 0.20 8.0 16 0.1–0.2 3–19–1 7.1 25.7 0.18 0.26 12.4 7 0.1–0.2 1–9–1 7.4 23.1 0.18 0.37 16.6 6 0.1–0.2 1–10–1 7.7 20.4 0.20 0.25 21.0 28 0.1–0.2 1–115–1 7.7 20.4 0.20 0.35 16.7 27 0.2–0.3 2–77–1 7.7 20.0 0.21 0.36 18.1 29 0.2–0.3 3–30–1 7.4 17.9 0.22 0.44 14.1 12 0.2–0.3 1–28–1 7.7 20.4 0.24 0.42 28.2 31 0.2–0.3 1–5–1 7.7 20.4 0.26 0.36 20.5 27 0.2–0.3 12–59–1 7.7 15.8 0.26 0.42 25.7 25 0.2–0.3 3–34–1 7.7 19.5 0.26 0.34 19.8 33 0.2–0.3 12–31–1 7.7 15.8 0.26 0.45 24.5 23 0.2–0.3 5–11–1 7.7 18.7 0.27 0.29 30.6 31 0.2–0.3 1–53–1 7.7 20.4 0.27 0.35 22.0 26 0.2–0.3 8–15–1 7.7 17.4 0.30 0.40 27.8 27 0.3–0.4 11–32–1 7.7 16.2 0.30 0.51 26.3 25 0.3–0.4 4–20–1 7.7 19.1 0.32 0.40 21.9 26 0.3–0.4 11–1–1 7.7 16.2 0.32 0.45 33.5 27 0.3–0.4 13–15–1 7.7 15.3 0.33 0.48 35.9 21 0.3–0.4 3–32–1 7.7 19.5 0.33 0.35 24.1 32 0.3–0.4 13–32–1 7.7 15.3 0.34 0.52 41.7 24 0.3–0.4 12–5–1 7.7 15.8 0.35 0.52 33.0 26 0.3–0.4 14–14–1 7.7 14.9 0.36 0.42 26.1 22 0.3–0.4 7–21–1 7.7 17.9 0.38 0.64 31.9 26 0.3–0.4 12–9–1 7.7 15.8 0.42 0.41 48.4 23 0.4–0.5 15–114–1 7.7 14.5 0.43 0.54 27.5 21 0.4–0.5 17–101–1 7.7 13.7 0.44 0.43 33.2 19 0.4–0.5 17–110–1 7.7 13.7 0.44 0.44 32.0 27 0.4–0.5 14–10–1 7.7 14.9 0.44 0.56 62.1 23 0.4–0.5 17–46–1 7.7 13.7 0.45 0.48 32.9 19 0.4–0.5 7–92–1 7.7 17.9 0.45 0.42 26.8 27 0.4–0.5 13–104–1 7.7 15.3 0.46 0.57 31.2 24 0.4–0.5 15–99–1 7.7 14.5 0.48 0.41 23.8 21 0.4–0.5 17–86–1 7.7 13.7 0.50 0.46 45.4 17 0.4–0.5 12–86–1 7.7 15.8 0.50 0.48 28.6 22 0.5–0.6 14–89–1 7.7 14.9 0.51 0.49 22.2 22 0.5–0.6 10–77–1 7.7 16.6 0.51 0.50 27.8 29 0.5–0.6 17–65–1 7.7 13.7 0.52 0.51 34.2 19 0.5–0.6 17–43–1 7.7 13.7 0.54 0.46 23.2 18 0.5–0.6 17–77–1 7.7 13.7 0.54 0.57 46.1 19 0.5–0.6 17–7–1 7.7 13.7 0.54 0.63 34.3 20 0.5–0.6 14–8–1 7.7 14.9 0.55 0.52 34.9 23 0.5–0.6 17–66–1 7.7 13.7 0.55 0.56 36.5 19 0.5–0.6 17–41–1 7.7 13.7 0.56 0.62 35.9 21 0.5–0.6 * P sa has been calculated for the fundamental period of the structure (T =0.38 s) A. Fathi et al. Structures 81 (2025) 110263 6 experimental and numerical frequencies. Finally, non-linear compressive and tensile behaviour was added to the properties using the concrete damaged plasticity (CDP) constitutive model [98]. The final material properties can be seen in Table 4. 5.3. Pushover analysis After calibration of the model, a pushover analysis was performed in the direction of interest (i.e., the weakest or transverse Y direction in Fig. 4. Examples of the simulated ground motions for Tabriz: a) 13–19–1; b) 5–8–1; c) 1–5–1; d) 14–14–1; e) 15–99–1; and(f) 17–41–1. A. Fathi et al. Structures 81 (2025) 110263 7 Fig. 7, as confirmed by previous studies on this structure [18,19]) by applying a uniform mass-proportional loading pattern aiming at a better understanding of the behaviour of the structure under lateral loading. Given the substantial influence of mass distribution on the seismic behaviour of Arge-Tabriz, attributed to its significant mass (~34,750 tons, based on the created model), the mass-proportional loading pattern was chosen for the pushover analysis. This decision was based on its ability to estimate load capacity with results sufficiently close to those obtained from nonlinear dynamic analyses, according to Endo et al. [99]. The adopted loading pattern is expected to predict the load capacity of the structure adequately close to the nonlinear dynamic analysis [99]. For this purpose, the analysis started by applying the self-weight of the structure using a “General Static” numerical step and then completed by imposing monotonically increasing lateral forces. Two separate analyses were carried out for the +Y and –Y directions, and in both analyses, the nonlinearities of the material and geometry were considered. Small (non-affecting the results) artificial damping was introduced to the model to overcome the numerical convergence problems. The pushover curves in Fig. 8 indicate that the structure shows higher force capacity and (nonlinear) stiffness in the –Y direction, i.e., it needs higher values of lateral loads to reach the same amount of drift compared to the +Y direction. The crack pattern of the structure for the weakest loading direction (+Y) is given by the maximum principal plastic strain (PE) contours of Fig. 9. Given the U-shaped geometry of the structure, the web (central) wall of the structure tends to outward overturn because of the tensile vertical cracks in the intersections of the web wall with wing walls (side walls) and horizontal cracks between two openings. Compressive behaviour is expected in the intersections of web wall with wing walls in the case of the inward overturning mechanism. As expected, the behaviour of the structure in the +Y direction is more critical than in the –Y direction. To ascertain the seismic demand of the structure, the capacity curve of the equivalent single-degree-of-freedom (SDOF) system has been obtained from the force-displacement curve of the building in the +Y direction as a multi-degree-of-freedom (MDOF) system. The Fig. 5. Distribution of PGA versus PGV for the selected records. Fig. 6. Southern view of Arge-Tabriz (Arge Alishaah). Fig. 7. A perspective of the 3D finite element model of Arge-Tabriz. Table 3 Comparison of the experimental and numerical values of the natural frequencies of the structure. Mode Number Frequency (Hz) Error Experimental Numerical Mode 1 2.60 2.60 0 % Mode 2 2.60 2.65 2 % Mode 3 3.40 3.35 1 % Mode 4 3.60 3.52 2 % Mode 5 5.00 4.98 0 % Mode 6 5.30 5.43 2 % Table 4 Material properties used in the numerical model. Density Compressive strength Young’s modulus Poisson’s ratio Tensile strength Tensile fracture energy ρ f c E ν f t G f kg/m 3 MPa MPa - MPa N/mm 1750 2.64 [97] 2350 0.2 [11, 97] 0.18 [97] 0.02 [11] Fig. 8. Pushover curves of the structure in the North-South (Y) direction. A. Fathi et al. Structures 81 (2025) 110263 8 transformation of the MDOF system to an SDOF system has been performed using the N2 method [100] following [3,8]. One of the most important tasks before the transformation procedure is determining the ultimate displacement of the structure, which remains an open point of debate. The displacement of the structure corresponding to a certain level of load drop in the post peak part of the pushover curve may be considered [101] or a certain level of drift regardless of the load level may be adopted [102,103]. Considering the shape of the pushover curves in the current analysis, which does not indicate a post peak branch clearly, the second procedure is adopted, and the ultimate displacement was considered as 3 % of the structure’s height, based on the recommendations of [102], following the work of [8]. After the definition of the ultimate displacement, mode 2 of the structure was selected as the most representative mode shape of the structure’s behaviour in the +Y direction; then, the capacity curve of the SDOF system was obtained by applying the transformation factor, or socalled participation factor, (Γ) to the capacity curve of the MDOF system. The mode shapes of Arge-Tabriz are available in the work of Hoveidae et al. [20] and a comprehensive explanation of the N2 method can be found in [100]. For the bilinear idealisation of the SDOF system, the procedure adopted by [3] and also recommended by the Iranian seismic code [104] has been used, i.e., the graph is approximated in a way to have the same stiffness in the ascending branch and the equal areas of the sections above and below the intersection of the actual and the idealised curve, up to the ultimate displacement. Then, the steps of the N2 method were followed to calculate the seismic demands of the structure for three different limit states, namely damage limitation (DL), significant damage (SD), and near collapse (NC) corresponding to a probability of exceedance of 20 % (return period of 225 years), 10 % (return period of 475 years), and 2 % (return period of 2475 years), respectively, over a period of 50 years. Table 5 lists a summary of different performance limit states. It is worth noting that the Iranian seismic code [104] presents the earthquake spectrum with a return period of 475 years, as the reference earthquake, which corresponds to the SD limit state. To obtain the spectra of the other two limit states (DL and NS), the procedure recommended by Eurocode 8 [106] has been adopted. The bilinear capacity curve is plotted against the elastic demand spectra of these three limit states, indicating the target displacement values for all three limit states in Fig. 10 with the corresponding maximum principal plastic strain (PE Max) contours presented in Fig. 11. The contour regarding the DL limit state has not been shown due to the very low values of plastic strains. As can be seen in Fig. 10, the target displacement for the DL limit state is close to the yield point of the structure, while the displacement for the SD limit state is almost twice the value for DL. The contours of PE Max for these two limit states (Fig. 11a) indicate adequate performance. The drift values regarding the first two limit states, 0.047 % for DL and 0.097 % for SD, confirm the hypothesis that the structure can withstand the reference earthquakes with return periods of 225 and 475 years. Considering the target displacement value of the NC limit state, which imposes a drift of 0.301 %, and considering the contours of PE Max, it seems that the main cracks, which will lead the web wall of the structure to overturn, start to form at this stage. The cracks are not fully developed and the drift value is not enough to force the structure to global collapse. 5.4. Time-history analysis The last step of this study is carried out using dynamic time history (TH) analyses based on the simulated ground motion records. First, adopting a static step, the structure is subjected to the vertical load of its self-weight while the earthquake loads, simulated ground motions, are applied to the structure employing a dynamic step. As mentioned in the previous section, the Y direction represents a more critical and vulnerable loading direction due to the U-shaped geometry of the structure, as confirmed by previous studies [18,19]. In this regard, the pushover analysis showed that loading in the +Y direction requires smaller lateral forces to achieve the same level of lateral displacement. This is because the toes of the wing walls and the intersections of the web wall with wing walls are subjected to tension under loading in the +Y direction, whereas they are in compression when the load is applied in the -Y direction. This tension-compression asymmetry makes the +Y direction the most critical for the structure. Therefore, the focus of the TH analyses in this study has been on the +Y direction. The non-linearity of the material and the occurrence of large deformations (non-linearity of geometry) are allowed during the analyses. The results of the TH analyses are discussed in the form of fragility curves and crack patterns. As a function of varying intensity measures, the probability of exceedance (POE) of a certain damage threshold is represented by fragility curves. Considering the past research on the Arge-Tabriz building and the usage of drift and Park-Ang damage index by Hoveidae et al. [20], these two indices have been considered as the seismic demand to calculate the fragility functions and plot the curves. Although both the drift and Park-Ang index are indicators for seismic damage evaluation, their bases are significantly different. Drift quantifies the lateral displacement of a given point in the structure normalised by the height and reflects the instantaneous global deformation state of the structure and is usually correlated with visible damages like cracking. On the other hand, Park-Ang damage index is a combined and more comprehensive seismic damage measure that considers effective parameters on the behaviour of buildings during an earthquake and integrates the damage caused by excessive deformations with the damage triggered by the cumulative hysteretic energy dissipation [20,107,108]. In the Park-Ang damage index, the deformation component works similar to drift while the energy component considers the impact of cyclic loading which is particularly important for ductile and non-brittle structures. The chosen intensity measure to derive the fragility curves is PGA, which is highly correlated with the damage patterns observed in rigid Fig. 9. A screenshot of PE Max principal contours for the pushover analysis in the +Y direction at a drift of 1 % (displacement at top of 0.36 m). Table 5 Summary of performance limit states and corresponding demand parameters. Limit state Return period (years) Probability of exceedance Description [105] DL 225 20 % in 50 years Minor structural damage; strength and stiffness mostly retained; nonstructural cracking possible; negligible permanent drift SD 475 10 % in 50 years Notable structural damage; reduced lateral resistance; non-structural components cracked; moderate permanent drift NC 2475 2 % in 50 years Severe structural damage; low residual strength and stiffness; most non-structural elements collapsed; large permanent drift; near failure A. Fathi et al. Structures 81 (2025) 110263 9