Numerical modelling for the seismic assessment of complex masonry heritage buildings: the case study of the Giralda tower

Existing heritage buildings can be especially vulnerable to earthquakes. They were designed only considering gravity loads and some of them are located in earthquake prone areas, such as the southwestern Iberian Peninsula. Besides, there is a high uncertainty in the definition of their constructive parameters and complex geometry. Due to that, it is paramount to develop accurate numerical models to obtain a reliable assessment of their seismic behaviour. Given this, the main objective of this manuscript is to analyse the seismic behaviour of the Giralda tower, located in Seville (Spain). It was declared as a Word Heritage Site of Outstanding Universal Value by the UNESCO in 1987. Seville has a moderate seismic hazard, but it has been shown that the soft alluvial strata amplify the seismic action. The tower has a brick and stone masonry structure, which was constructed in several construction phases. A 3D Finite Element Model of the tower has been developed using OpenSees software, employing a 3D CAD model. Modal analyses and nonlinear static analyses have been applied to calibrate and to assess the tower’s seismic behaviour. The results showed significant differences in function of the load pattern. It should be remarked that the boundary conditions have a notable effect upon getting a good calibration of the model. Regarding the damage, it has been found to match the historic records: the ramps would be ruined and, in the outer wall, it would be concentrated near the openings, especially close to the belfry.


Introduction
Europe has a vast Cultural Heritage, which produces a significant economic and cultural activity. Built heritage located in urban centres are in this set of valuable cultural assets. Therefore, preserving its cultural heritage is a current and challenging issue for the European Union (EU). In that sense, all the countries are committed to preserving and to improving their cultural heritage through policies and programmes. The structural security assessment of heritage buildings is a paramount task. This is especially relevant in seismic areas such as the Iberian Peninsula (Amaro-Mellado et al. 2017a).
Seville is located in the southwestern Iberian Peninsula. It has a moderate seismic activity, but it is located on soft alluvial strata that are known to increase the effects of earthquakes (Amaro-Mellado et al. 2017b;Requena-Garcia-Cruz et al. 2022). The city is affected by large earthquakes (Mw > 7) of long-return periods, which makes the population unacquainted with them (Fazendeiro Sá et al. 2016). The seismic hazard of the region is due to the contact between the Eurasian and African tectonic plates (Amaro-Mellado et al. 2017a) and it has been affected by several historic earthquakes, such as the 1755 Lisbon (Mw = 8.5) and the 1969 earthquake (Mw = 8) (Sá et al. 2018).
The city has a great cultural heritage, which affects its historic and social identity. In this context, the Giralda tower is the most representative building of the city. Originally, it was built as a minaret for the Islamic major mosque in 1198 (Jiménez Martín 1998). The tower has undergone different construction phases and significant modifications over time. It was declared as a Word Heritage Site by the UNESCO in 1987, with the maximum level of protection: Outstanding Universal Value. In relation to its seismic risk, it is known that it was affected by several historic earthquakes in the past. For example, it was recorded that the 1755 Lisbon earthquake (Gentil 1989) seriously damaged the tower.
The structural assessment of architectural heritage buildings is a complex task that requires a multidisciplinary approach. Knowledge about the evolution of the building in its history is essential: geotechnical properties, phases of construction, rehabilitation interventions, materials and the geometry of the structure. Whit this in mind, the seismic vulnerability has been studied in a vast number of research works due to the uncertainties in the performance and in the complexity of this type of buildings (Degli Abbati et al. 2019). Seismic assessment plays an important role in the preservation and in the safety 1 3 of architectural heritage, allowing the development and implementation of appropriate retrofitting strategies. These buildings were not constructed considering seismic actions and they are very vulnerable to lateral loads. In this respect, ancient towers are very remarkable case studies because they have some geometrical features that increment the seismic vulnerability (Sarhosis et al. 2018;Valente and Milani 2018). Furthermore, past earthquake damage assessments have shown that ancient buildings are very vulnerable during a seismic event (D'Ayala and Paganoni 2011;Parisi and Augenti 2013;Cattari et al. 2014;Saisi and Gentile 2015).
The main difficulty in the structural assessment of ancient buildings is to obtain the mechanical properties of the materials, especially in masonry buildings, due to their heterogeneous properties and the numerous parameters (Kržan et al. 2015;Cattari et al. 2021). Because of the uncertainties in the mechanical parameters of ancient masonry buildings, an in situ survey campaign is required to calibrate the numerical model (Ponte et al. 2021). These surveys campaigns are composed of a Non-Destructive Test (NDT) given the high level of protection of the ancient buildings. In this regard, the NDT have been discussed and illustrated in the literature (McCann and Forde 2001;Boschi et al. 2019). The ambient vibration test is the main experimental method used to evaluate the dynamic behaviour of the structures and to calibrate the numerical models (Gentile and Saisi 2007;de Silva et al. 2018;Nohutcu 2019). Moreover, several research works analyse the different mechanical properties of masonry walls with the aim of their being able to be used as reference values. In (Boschi et al. 2019) the results of 105 in situ tests performed in a great number of brick masonry typologies were shown. In (Kržan et al. 2015), a range of mechanical properties reference parameters was provided for a great variety of stone and brick masonry, based on experimental tests and data from the literature.
The numerical modelling of historic monuments, generally characterised by large dimensions and complex and irregular geometric features, is a difficult task. In this vein, the nonlinear Finite Element Method (FEM) is one of the most advanced strategies to simulate masonry structural behaviour (Requena-Garcia-Cruz et al. 2023). Considering this, different modelling approaches and analysis methods are considered for the seismic assessment depending on the characteristics of the structure and the level of accuracy of the analysis and the results: micro-modelling, macro-modelling and homogenisation (Paloma et al. 2010;D'Altri et al. 2020;Cattari et al. 2021). In that sense, the correct definition of the modelling strategy is essential. Macro-modelling is more commonly used for the global behaviour of a full-scale structure. Moreover, this is applied in the case of solid thick walls so that the stress across or along the macro length will be fundamentally uniform (Pineda 2016). In this regard, the masonry is modelled as an isotropic or anisotropic homogeneous continuum (Degli Abbati et al. 2019). One of the most complex tasks in this type of modelling is to define the material's anisotropic inelastic behaviour (Petracca et al. 2017b). Employing a solid 3D Finite Element (FE), which represents the mechanical behaviour of the masonry in a continuous form, enables the maintenance of cost-effective computational efficiency (Tiberti et al. 2016). In this context, there are a great number of research works in which different types of modelling are applied (Milani et al. 2012;Malcata et al. 2020;Milani and Clementi 2021).
The seismic vulnerability of ancient buildings is generally performed considering nonlinear analyses: static pushover and/or incremental dynamic, e.g., (Lagomarsino and Cattari 2015a). Furthermore, the different failure mechanisms can be assessed applying the kinematic limit analysis approach (Milani 2019). Static and dynamic analyses have been widely used in research and engineering works. Furthermore, the different seismic codes propose both as valid test methods for the analysis of existing buildings (European Union 2005a;

3
Ministerio delle Infrastrutture e dei Trasporti 2018). The implementation of the NLSA in 3D numerical models is a proper, and commonly used, method to correctly predict the global damage in the numerical model and active failure mechanisms (Milani and Clementi 2021). There are some research works that compare the application and the results obtained for both methods in ancient buildings (Milani et al. 2012;Lagomarsino and Cattari 2015b;Milani and Valente 2015). The NLSA results represent the envelope of all possible structural responses and it is possible to apply in case of global behaviour analysis, taking into account the inherent limitations (Milani et al. 2012). The NLDA are more accurate methods to simulate the structural behaviour of ancient masonry buildings (Ferrante et al. 2021). It should be noted that dynamic excitation activates a higher vibration mode, which could be associated with damage in the weak upper levels of the tower, especially in towers with large hollows at the top (Valente and Milani 2018). However, the NLDA needs a high computational effort and calculation time, compared to the NLSA. Taking into account all issues and addressing the main objective of this work, the NLSA has been used to validate the upgrade FE model.
The seismic behaviour of masonry towers is an important issue due to the large number of structures of this typology in the built heritage around the word. In this respect, many works employ the FEM for the analysis of masonry towers, such as (Valente and Milani 2016), where eight masonry towers were analysed, using a detailed 3D FEM. This work showed a good agreement between the results of the nonlinear pushover and the dynamic analyses. However, the results of the pushover analysis are more conservative. In (Torelli et al. 2020), the seismic assessment methods proposed by the Italian codes for the assessment of the seismic vulnerability in cultural heritage buildings were analysed in a masonry tower. In (Shehu 2021), the limitations and the possibilities of the pushover analysis with the FEM, using three masonry towers as a case study, were discussed. In (Peña et al. 2010), three different numerical FEMs of old masonry towers were carried out with several complexity levels, to compare the results.
The different reference seismic codes provide a relevant framework for structural analysis. In this context, the European reference standard for assessing earthquake resistance structures is the Eurocode-8 Part 3 (European Union 2005b). In the same way, the Eurocode-6 (European Union 2005c) and the Italian code NTC-2018 (Ministerio delle Infrastrutture e dei Trasporti 2018) provide valuable data about masonry material, which are important to consider in the numerical definition. What is more, the Spanish seismic standards have also been considered in this work (Spanish Ministry of Public Works [Ministerio de Fomento de España] 1998, 2002Gobierno de España 2009).
The main aim of this manuscript is to numerically model the Giralda tower by applying a 3D FEM for its seismic assessment. For this purpose, the OpenSees open-source software (McKenna et al. 2000) and the STKO software (Petracca et al. 2017a) have been used. The 3D FEM has been developed applying macro-modelling elements (D'Altri et al. 2020) due to the complexity and to the size of the tower. A specific method has been applied to resolve the uncertainties which are present in the modelling of complex heritage buildings. Besides, a sensitivity analysis has been performed to define the mesh size, mechanical properties, and the pushover load pattern. The modal and the nonlinear static analyses (NLSA) have been carried out to obtain the seismic behaviour of the tower. This type of study is an important opportunity to offer useful information and knowledge related to complex heritage buildings. Despite the conservation projects carried out on the building, this is the first advanced numerical modelling work for its seismic behaviour assessment. Although the tower has a great cultural value and it is exposed to seismic hazards, this study has not yet been completed, but important highlights about the numerical modelling have already been provided, that can be useful for future research works and for the conservation projects done in the monument. Furthermore, it is worth highlighting that, according to the authors' knowledge, this is the first time that the OpenSees framework, with macro-mechanical modelling has been used for the seismic assessment of a complex masonry heritage building. This study has: (1) implemented advanced numerical methods to study the structural characteristics of the tower for the first time; (2) developed an updated FE model for future assessments on the ancient building; (3) made it possible to depict, for the first time, the damage pattern of the building.
The paper has been structured as follows. The methodology followed has been described in Sect. 2. In Sect. 3 the historic data, the geometry and the structural characteristics of the tower are presented. Then, Sect. 4 includes an analysis of the uncertainties and the in-situ test campaign. The numerical modelling of the building has been presented in Sect. 5. In Sect. 6, the modal and the nonlinear static analyses are performed. Later, the results and the discussion are set out. Finally, the conclusions are drawn (Sect. 7).

Methodology
In this section, the methodology followed to develop the research is presented. The flowchart is shown in Fig. 1.
Firstly, the building information has been gathered. In this connection, a historic analysis of the tower and its surroundings has been developed. During this phase, different research documents and projects (Jiménez Martín 1998)(Barrios Padura et al. 2012)(Jiménez Martín and Cabeza Méndez 1988) have been consulted and evaluated. This phase is important to understand the historic evolution of the building, detecting the different construction phases and interventions carried out.
Secondly, the geometrical information available of the tower (Jiménez Martín and Cabeza Méndez 1988; Almagro and Zúñiga 2007) has been compared with several in situ surveys, in order to build an accurate geometry model of the building. For this purpose, a visual inspection, a laser measurement and a digital image analysis, have been performed. Also, an ambient vibration test campaign has been done to obtain the dynamic characterisation of the building. Furthermore, an accurate analysis of the different uncertainties has been developed. In this step, different variables have been identified to be applied in the sensitivity analysis and in the calibration of the numerical model. Then, a simplified 3D CAD model has been worked out to make the numerical model for the structural assessment. Rhinoceros software v7 (McNeer 2021) has been used for this. Then, the 3D CAD model has been exported to the STKO Software (Petracca et al. 2017a) to obtain the base geometry of the building. In this software, the definition of the 3D FEM model has been performed. The solid mesh, the boundary conditions and the properties of the materials have been defined in it. Lastly, the OpenSees open-source software (McKenna et al. 2000) has been used for the calculations.
Afterwards, a modal analysis has been carried out. This has been employed to compare and to analyse the main dynamic characteristics of the tower with other similar research works and in situ surveys. In fact, the results of the in situ ambient vibration tests have been used to check that the definition of the FEM model is adequate and to calibrate it. Then, nonlinear static analyses have been done in both directions, to obtain the capacity curves of the building, applying different load patterns. In addition, a sensitivity analysis concerning the mesh size and mechanical properties of the numerical model has been conducted. Finally, the different capacity curves and the damage patterns have been analysed and compared to identify the global seismic behaviour of the tower.

Historic data
The Giralda tower is the most emblematic building of Seville. Furthermore, it is one of the most visited monuments in Spain. It is an excellent example of the different cultures and periods during which the city has lived throughout its history. Originally, it was built as a minaret for the Islamic major mosque from 1184 to 1198. It is important to highlight that it was the tallest minaret built during the Islamic period, even considering other prominent ones such as the Marrakech and Rabat mosques (Jiménez Martín 1993).
The construction of the Major Mosque of Seville was started in 1172 by Almad ibn Bāsu. Seville was the capital of Al-Andalus during the Almohad period (1147-1248). The Major Mosque was inaugurated in 1182 by the caliph Abū Ya'qūb. This caliph enclosed the building and built its minaret. The foundation of the minaret was built in 1184. Then, the construction of the alminar started in 1188 after four years of non-construction. The alminar was concluded with four golden spheres (Yamur) on its top in 1198 by the caliph Abū Ya'qūb Yūsuf (Jiménez Martín 2007) (Fig. 2).
The Christian Cathedral was built on the site of the former Islamic Mosque. It was designed in Gothic and Renaissance styles, and it is the largest Gothic building in Europe. The minaret and the yard of the ancient Islamic building are the only parts that are preserved today. The tower has undergone different construction phases and relevant modifications. The most important was the Renaissance style construction of the Bell tower in by Hernan Ruiz in 1568 (Jiménez Martín 1998). The building thus acquired a similar appearance to its current one (Fig. 1b). The tower was affected by several historic earthquakes, such as the 1356 earthquake, when the four Islamic golden spheres fell and the Lisbon earthquake in 1755 (Mw = 8,9), when it was severely damaged (Gentil 1989).
In the UNESCO declaration, it was highlighted that the Giralda tower is a masterpiece of architecture, which is an important example of cultural syncretism (https:// whc. unesco. org/ en/ list/ 383). Moreover, it influenced the construction of many towers in Spain and around the world.

Geometrical and structural features
The tower was designed with two parallel brick ceramic masonry walls. The floor has a square form of 13.60 × 13.60 m and it is 94.69 m in height. Compared to other masonry towers, the masonry wall increases its thickness at the upper levels. The outer wall is around 2.00 m thick in the lower levels, and around 2.30 m thick in the upper levels. The inner body has a square floor of around 6.15 m side, and its wall thickness changes depending on the tower's height, ranging from 1.30 m on the lower levels to around 1.50 m on the upper levels. In the core of the tower, several brick masonry vaults define seven different rooms (Fig. 3). It is important to highlight that the walls are entirely made of solid ceramic brick, as it has no filling or other materials inside. This was tested in 1985 when several boreholes were done in different areas of the walls (Pérez 1985).
The connection between both walls is made through the ascent ramps. These are comosed of several brick masonry vaults and compact limestone concrete, with incremental thickness in order to build the slope. The thickness of the ramps ranges from 1.10 to 1.40 m.
Finally, the foundation of the tower is formed by four courses of calcarenite blocks, which make a total thickness of around 2.50 m. This rests on a foundation slab with a greater thickness in its centre (around 3.00 m in the centre and 1.00 m in the perimeter (Romero ).

Analysis of the uncertainties
The analysis of a building of cultural heritage has several uncertainties due to numerous issues: its complex geometry, the boundary conditions, the localisation of the materials and the characterisation of their mechanical parameters. In this work, these uncertainties have been faced by applying several methods and tools that are described hereafter.

Complex geometry
Firstly, this issue has been developed by an accurate historical analysis of the literature (Jiménez Martín and Cabeza Méndez 1988), the images and the blueprints (Almagro and Zúñiga 2007). Then, the information gathered was completed and compared with several in situ measurements (a laser measurement, images analyses and a visual inspection of the building). A precise geometry model of the building has been elaborated with this information.
Then, the geometry has been simplified to compute a 3D CAD model (Fig. 7). This model has been developed considering several issues for the exportation of the FEM. For this purpose, the Rhinoceros software (McNeer 2021) has been used. In this regard, it should be observed that the Giralda´s outer wall is formed by the structural wall and a non-structural wall in the openings to carry out the decorations (Fig. 4). Only the hollows in the structural wall have been considered in this analysis. The definition of the different openings has been done after field surveys and an analysis of the hollow configuration. Generally, the towers have some main typological features, such as slenderness, inclination, and presence of openings and belfry. These geometrical properties could affect in its structural behaviour. The influence of some geometrical properties on the seismic response of the towers is quantitatively investigated in (Valente and Milani 2018). These geometrical features have been considered in the analysis of the modal and the NLSA results.

Different materials and location.
Ancient buildings were built using several types of materials with irregular locations. In this context, an accurate analysis of the historic documentation and the literature is important. In this phase, the different construction stages and the materials employed have been identified and analysed. Besides, the past characterisation project executed in the building (Pérez 1985) has been consulted to identify the type of materials, their physical properties, the composition of the structural elements and their characteristics. Finally, the information collected has been analysed and compared with a visual inspection of the building and the image analysis.
Hence, four samples were extracted in different locations: three in the brick masonry walls and one in the fourth ramp (vertical sample). The retrieved samples had a diameter of 10 cm and a length of 1.20 m. These samples helped to identify the composition of the walls and the ramps. Also, several laboratory tests were carried out with the aim of identifying the materials. In this vein, the composition of the ceramic brick masonry and the lime mortar was analysed (Pérez 1985).
The raw material of the ceramic bricks is composed of mica, quartz, feldspar and calcium carbonate. This is the composition of the alluvial clays of this area. There are two types of ceramic brick, each with a different concentration of illite and carbonates. Different tonalities (reddish and light) have been identified because of the concentration of the chemical elements and the baking temperature of the brick. A small amount of gypsum has been identified in both. Lime mortar is composed of lime (53%) and sand (47%). It should pointed out noted that the lime is 80.8% pure.

Constructive details of the structural elements.
It is important to simplify some aspects to build an accurate numerical model with a costeffective computational efficiency. In this regard, the first step was to identify the different structural elements for their modelling to reproduce their adequate structural behaviour. In this sense, it is important to organise an NDT campaign (McCann and Forde 2001) to identify the composition of the structural elements. For this purpose, thermography and georadar can collect some valuable information about wall composition (Ponte et al. 2021). In that case, an accurate analysis of the literature (Esbert et al. 1988.;Jiménez Martín and Cabeza Méndez 1988;Tabales Rodríguez 1998;Barrios Padura et al. 2012) and the previous characterisation project (Pérez 1985) has been done. Furthermore, this information has been compared and completed with a visual inspection of a building and a digital image analysis.
In this context, the samples analysed in the past project (1985) (Pérez 1985) provided a fundamental information about the composition of the brick masonry walls and ramps. It should be remarked that the masonry walls are made up of a single leaf of solid clay brick and lime mortar. Their composition has been shown in Sect. 4.2. The dimension of the clay bricks is 32 × 15x5.5 cm. It is important to note that the thickness of the bricks varies from 3 to 6 cm. The dimension of the mortar joints ranges between 2.5 and 4 cm. This dimension is higher in the core of the wall, with the aim of regularising it.
It should be mentioned that the wall was built entirely of brick masonry. It presents a homogeneous composition in the visual inspection, due to which its structural behaviour is considered uniform. Furthermore, the wall has a great rigidity due to its high thickness (2.00-2.30 m outer wall and 1.30-1.50 m inner wall). Therefore, the presence of out-ofplane failures is very difficult.
The base of the ramps is a brick masonry vault, on which the slope formation was built (Fig. 5). The upper layer of the ramps is made up of clay bricks and lime mortar with a thickness of 10 cm. Below, there is a layer of compact limestone concrete with a thickness of 12 cm. This concrete is composed of gravel, slice clay bricks, rubble, clay, etc., all of them conglomerated with lime. The next layer is a less compact limestone concrete. This concrete is made up of ceramic elements, pipe clay and organic material, all of them mixed with lime. The thickness of this layer is 15-17 cm. The next layer is composed of silty sand with lots of organic material. This last layer has different thickness with the aim of creating the slope of the ramps. Finally, the vault was created with brick masonry.

Mechanical parameters of the materials
The definition of accurate and reliable mechanical properties of the masonry material is a challenging task. In the case of ancient structures, this becomes even more complicated because destructive and semi-destructive in situ testing is limited or often prohibited. In that sense, the characterisation of the mechanical properties of the materials has been carried out based on reference values proposed in codes, in other similar research works and in the results of previous tests done in the building.
In this work, the reference parameters for the mechanical properties have been obtained from different sources. First, the values of a previous characterisation project, performed in 1985 in the tower (Pérez 1985), have been collected. In this project, several brick masonry samples were analysed by applying compressive laboratory tests. These values have been used as a reference to compare with other research works and codes (Table 1). Furthermore, masonry units (ceramic brick and lime mortar) were independently identified and analysed. In addition, the compressive strength has also been determined by means of the equations provided in the Eurocode 6 (European Union 2005c), . Where f ck is the characteristic compressive strength, K is a constant depending on the type of unit and mortar, f b and f m are the compressive strength of the masonry unit and the mortar, respectively. The tensile strength (f t ) has been obtained in function of the compressive strength f cp (f t = 6%f cp ), where f cp is the peak compressive strength. In Table 1, the reference preliminary mechanical properties of the materials used, based on other research works and codes, are shown. These are mechanical properties of the masonry, which have been defined first in the numerical model.

Dynamics of the building
Due to the uncertainties in the mechanical parameters of ancient masonry buildings, it is necessary to carry out an in situ survey campaign to calibrate the numerical model. An ambient vibration test is one of the main experimental methods used to evaluate the dynamic behaviour of large-scale structures (Ponte et al. 2021). The extraction of modal parameters (natural frequencies and modal shapes) from the experimental data was done by applying an Operational Modal Analysis (OMA). The ambient vibration test is preferred for analysing historic buildings because it is non-destructive. In carrying this out, the vibration of the structure is caused by the environmental excitations without the necessity of equipment to force the vibration. This type of test has been widely used to assess the dynamic behaviour of ancient buildings (Gentile and Saisi 2007; Nohutcu 2019; Malcata et al. 2020). The free ambient vibration has been carried out on the rooftop, in the belfry and in the 6 th floor interior room. The accelerometers were located at the corners of these levels in order to obtain the main vibration modes of the tower. The reference accelerometer was located on the rooftop to correlate the set-ups (Fig. 6 a blue mark). A forcebalance triaxial accelerometers GMSplus measuring system from GeoSIG Inc. has been employed for this test (Fig. 6). The sensor has a high sensitivity and can reach low-noise individual 24-bit Δ − Σ ADC per channel.

Adjacent buildings iteractions
The interactions with adjacent buildings are an important issue in the seismic behaviour of ancient buildings (Degli Abbati et al. 2019). Normally, these types of buildings have undergone several construction phases, and so have the buildings around them. Originally, the tower was built as an isolated building. It was built on the north side of the old Islamic Mosque. There was a structural joint between both buildings. However, this condition changed with the construction of the Gothic Cathedral. In this respect, two walls of the Cathedral are in contact with the west and south façades of the tower (Fig. 7a). This iteration would affect the tower's dynamic behaviour.
This issue has been analysed in this work in order to check the change in the dynamic behaviour of the tower. With that in mind, some equivalent portions of its adjacent walls have been introduced into the model (Fig. 7b and c), like (Castellazzi et al. 2018). It should be observed that for modelling the boundary conditions that there is a structural joint between the Tower and the Cathedral (Fig. 7). The only structural elements in contact are the two walls modelled. These walls have been modelled with the same thickness and hight. However, they have been developed with an equivalent length (Fig. 7a). This fact allows reproducing the boundary condition, keeping a low computation effort and a simple model. Nevertheless, the stiffness of these adjacent walls has been calibrated, applying an iterative process. In this regard, the elastic modulus (E) of this masonry material has been changed in order to fit the result of frequencies (Hz) between numerical model and in-situ measurements (Sect. 6.1). The results are analysed and compared with the isolated tower condition.

The numerical model of the building
In this section, the steps followed to perform the numerical model of the Giralda tower are described in detail.

Geometrical CAD model in Rhinoceros
The 3D CAD model of the building geometry has been developed according to the available information and in situ measurement surveys (Fig. 8). This model has been done with the Rhinoceros v7 software (McNeer 2021), which is a NURB-based 3D modelling tool.
Several items have been considered to develop an in-detail geometry model. As a first step, an accurate geometrical analysis of the tower has been performed, considering the historic and available information about the tower. Furthermore, the historic drawings published in (Jiménez Martín and Cabeza Méndez 1988) and the architectural atlas of Seville Cathedral (Almagro and Zúñiga 2007) have been used. Then, in situ measurements and image analysis have been performed to complete and to adjust the geometry survey. In this regard, several items were considered: the wall thickness, the openings dimension, the different parts of the tower described above, the areas of different materials, the location of the ramps, etc.
The complex geometry of the tower was simplified to prepare the geometrical model. The 3D layout (Fig. 8) has been defined as follows: (1) the inner body, considering its levels and its openings; (2) the ramps, which have been modelled as horizontal planes, that (3) the outer wall, considering its openings and its thickness (increases in height); (4) the belfry and the upper body level have been modelled considering the different structural elements and the different materials. In this connection, the brick masonry vaults have been modelled as horizontal elements with equivalent thickness. This simplification is accepted in the study (Castellazzi et al. 2018). In the literature, some simplification of the complex geometrical features is commonly accepted with the aim of obtaining efficient FE models without excessive complexity (Milani et al. 2012;Valente and Milani 2018).
With the aim of improving the exportation and the definition of the FEM in the preprocessor and the analysis of the results in the post-processor of the STKO software, the model has been divided into different parts. First, it has been divided into four symmetric parts in plan to easily access the central nodes of the different floors. Second, it has been divided according to the different levels of the inner rooms (Fig. 8).

Finite Element Model (FEM) generation
The definition of the 3D FEM ( Fig. 9) has been done using the pre-and post-processor STKO for OpenSees. The simplification of the geometry has been performed to limit the computational effort and the meshing problems. In past research works, it is adopted and accepted to do some model simplifications to carry out a global behaviour assessment of a structure (Lagomarsino and Cattari 2015a;D'Altri et al. 2020). The X and Y axes of the numerical model correspond to the west-east and to the south-north cardinal directions, respectively. The 3D FE meshes have been generated from the CAD-based geometry model developed in Sect. 3.2.1. The 3D FEM has been modelled as a macro-modelling approach with solid elements, due to the complexity and to the size of the tower. These The macro-modelling approach has been meshed as a non-structured tetra mesh, applying four-node tetrahedron solid elements (Fig. 9). The mesh has been developed to reproduce the global behaviour of the tower. A mesh-sensitivity analysis was analysed, by modelling the tower with the same solid models but different structured mesh sizes (Fig. 10). The structural element size and the number of tetrahedron solid elements in the thickness of the outer wall have been considered to define these mesh sizes. The coarse mesh size has been defined to deal with at least three tetrahedron elements in the outer wall thickness, as in (Castellazzi et al. 2018).
The materials have been defined in the numerical model by applying the specific material implemented in the OpenSees software (Sect. 3.2.3). As mentioned above, four different materials have been considered: the unreinforced ceramic brick masonry wall, the stone masonry located at the base of the tower, the stone masonry of the belfry and the upper body, and another for the ramp's elements.

Material model
The definition of the non-linear material model for the masonry is an important and complex issue in the seismic analysis of ancient masonry buildings. The masonry is a composite heterogeneous material with an anisotropic behaviour, which depends on several features, such as, texture, typology of units, regularity, presence of multi-layer walls, etc. Despite that, in the macro-modelling approach, the material is defined applying advanced homogeneous and damage models to reduce the computational time of the analyses. Furthermore, the use of this type of material models, adjusting the parameter for the masonry behaviour is commonly accepted in the literature (Valente and Milani 2016). In the latter approach the mechanical behaviour of the masonry has two main features: its low tensile strength and its compression behaviour.
To that effect, the nonlinear mechanical behaviour of the masonry has been defined using the 3D Tension-Compression damage model (Petracca et al. 2017c, b), which is implemented in the STKO and the OpenSees software. This constitutive model, based on continuum damage mechanics, is used to simulate the behaviour of quasi-brittle materials, such as masonry, in the nonlinear analyses. It reproduces the nonlinear behaviour of the masonry material under shear and controls the dilatancy. This material model develops two independent failure criteria for the tensile and for the compressive parts, assuming different damage parameters. Besides, it has independent values of the ultimate stress and strain. In this context, this material reproduces two damage parameter indexes d + and d − : one for the positive part of the strain tensor (tension) and the other for the negative part (compression) (Eq. 1). Each of these affects the positive ( + ) and the negative ( − ) component of the effective stress tensor ( ). The results of these damage indexes are plotted in the range of 0 (intact material) to 1 (completely damaged material).
The failure criteria of this material model have been described in (Petracca et al. 2017b). In this regard, equivalent stresses τ + and τ-have been defined to identify the reciprocating loads conditions. The failure criteria, in compressive and in tensile conditions, are calculated in Eqs. 2 and 3, respectively. These are based on the plastic-damage concrete model defined by (Lubliner et al. 1989).
The variables and are obtained from: where: • I 1 is the first invariant of the effective stress tensor. • J 2 is the second invariant of the effective deviatoric tensor. • max is the maximum effective principal stress. • K 1 changes the size of the compression surface to control the shear behaviour of the material. This ranges from 0 (Drucker-Prager's criterion) to 1 (Lubliner's criterion). • K b is the ratio of the bi-axial strength to the uniaxial strength in compression.
• > 0 define the strength increases, in case of the stress state of the triaxial compression.
In Fig. 11a, the surfaces of the tensile and the compressive damage have been plotted in the 2D plane-stress scenario. It should be pointed out that the compression and the tensile surfaces are defined for any stress state, so, they must be deactivated under certain conditions. The compressive surface can be only developed when at least one of the main stresses is negative. The tensile surface can only evolve if at least one of the main stresses is positive. These conditions are implemented in the definition of each surface with the Heaviside function H (x) (Eqs. 2 and 3) (Petracca et al. 2017b). Furthermore, the tensile and the compressive surfaces might be active when the principal stresses have different signs. In that case, the positive and the negative stresses are affected by the tensile damage and the compressive damage, respectively. In this damage model, an anisotropic scaling of the effective (elastic) stress tensor is applied to obtain the damage stress tensor, as demonstrated in a previous work (Pelà et al. 2011).
In the case of pure shear distortion, the d + /d-damage model's damage stress σ has a different evolution law, in the case of compressive and tensile stress. The increase in compression damage is smaller than that of the tensile damage. Furthermore, the tensile and the compressive damage indexes do not have the same value due to the existence of two different failure surfaces.
This material model has implemented a formulation that includes the possibility of calibrating the dilatancy behaviour under shear stress. In this regard, the size of the compressive surface in the tension/compression quadrants (Fig. 11a) is controlled by the parameter k1 (Eq. 2). Therefore, this parameter implicitly controls the dilatancy of the  (Petracca et al. 2017b) model, considering that the larger compressive surface (compared to that of the tension surface), the higher the dilatancy. Figure 12 shows the stress-strain diagrams for this material model in uniaxial tension (a) and compression (b). This constitutive law is like other damage-plasticity material implemented in other research works about ancient towers (Valente and Milani 2018;Milani and Clementi 2021). On the one hand, the tension uniaxial law is characterised by the exponential softening curve shown in Fig. 12a. This constitutive law shows a linear elastic until the ultimate tensile strength (ft) is reached. Then, a softening phase is developed. On the other hand, the compressive uniaxial law is defined by a hardening/softening law (Fig. 12b). This curve is composed of several parts: (1) a linear elastic part, (2) a hardening part until the peak strength (f cp ) point is attained, (3) parabolic softening, (4) softening transitions to the residual strength (f cr ) and (5) a final constant residual stress part. Parts 2, 3 and 4 are quadratic Bezier curves, each having three control points. These control nodes define the shape of the curve and its tangents at the end position (Petracca et al. 2017c).
It is important to highlight that the compressive fracture energy (G c ) is not the total area under the uniaxial compressive curve. Hence, Gc is only the part after the peak-strength, that needs to be regularised (Petracca et al. 2017b). In this work the following equations have been used to obtain the real fracture energy values of the material (G t = 0.073·f cp^0 .18; G c = ((f cp /f t )^2) · G t )) from (Comité Euro-International du Béton 1990), where Gt is the tensile fracture energy, f cp is the compressive strength and f t is the tensile strength. Furthermore, the different consideration about that in (Lourenco 2013) has been applied.
The total energy dissipation is proportional to the FE size (Petracca et al. 2016). Therefore, the energy dissipated by the elements is dependent on their size. In this vein, the energy dissipated is higher when the size of the mesh element is larger. However, this material model has implemented the fracture energy-based regularisation to consider the characteristic length of the FEM. In that sense, it is important to define the regularisation of the material in order to make the numerical model be mesh size independent. In this model, the regularisation of the fracture energy is calculated by dividing the input real fracture energy by the characteristic length of the element (l ch ) (Petracca et al. 2017c). Therefore, the stress-strain curves of the material are adapted based on the mesh size of the FEM. With this, the global response of the material will be independent of the mesh size.

Modal analysis
First, a modal analysis has been performed to obtain the main vibration modal shapes and the periods to analyse them. In this respect, three fundamental global modal shapes have been identified in the numerical analysis. Modes 1 and 2 are translational in the X and the Y direction, respectively. The third vibration mode is rotational (Fig. 13 and 14). The 3D FEM of the tower has been calibrated using the results of the in situ survey. With this aim in mind, the values of the periods and frequencies of the numerical FEM have been compared with the ambient vibration test, described in Sect. 4.4.
The calibration of the FEM model has been carried out by applying an OMA test. For this purpose, the Artemis Modal Pro software (Structural Vibration Solutions A/S 2019) has been used. The experimental dynamic parameters of the tower have been extracted employing the Enhanced Frequency Domain Decomposition (EFDD) and the Stochastic Subspace Identification (SSI) technique. The vibration modes are determined by selecting peaks in the response spectral density functions for the set of performed tests. The different modal shapes obtained in the OMA test have been compared with the numerical ones in order to calibrate the 3D FEM (see Figs. 13 and 14). Three global modal shapes have been identified in both analyses: two bending modes and a torsional mode. With a view to calibrating the numerical model, the three vibration modes have been compared (Figs. 13  and 14).
An iterative process has been applied to calibrate the numerical model. For this, Young's modulus (E) and the density (ρ) of the different masonry materials have been modified iteratively to match the results with the in situ measurement. The main objective is to minimise the difference in the natural frequencies (or the periods) (Fig. 13). Furthermore, several variations of the Elastic modulus have been considered in the calibration process, based on the results of the previous characterisation project (Pérez 1985): (1) applying to the brick masonry an average value of f cp , and (2) considering several f cp in the brick masonry according to the localisation in the tower (bottom, medium and top).
The E values have been obtained based on the f cp for each masonry material ( Table 2). The masonry stiffness changes significantly in line with its compressive strength, as shown by the values proposed by several codes (European Union 2005c; Ministerio delle Infrastrutture e dei Trasporti 2018) and the literature (Kaushik et al. 2007;Paulay and Priestley 1992;Mohamad et al. 2007;Costigan et al. 2015). As Eurocode-6 (European Union 2005c) described, the elastic module (E) is linearly related to the masonry compressive strength (f cp ) (E = k·f cp ). With that in mind, codes and literature works have applied this equation with different K values. Several equations to calculate the E, in terms of f cp , have been tested in the calibration (Table 3). Thus, a set of different models is defined, taking into account several modifications in the parameters explained above. The value of E has been varied from 700 to 1000f cp and Poisson´s ratio has been set constant (ν = 0.2). What is more, different configurations in the model have been tested (Fig. 13). Based on the results reached, it was verified that the numerical model defined with 800f cp , and boundary conditions is the closest to the experimental model. The first modal shape is the same (0.66 Hz). The second is very similar, with only a 2.9% difference. In addition, the torsional modal shape (mode 3) is also very similar, with an 11.9% difference. Figure 14 shows that the calibrated and the experimental frequencies are closely related. Furthermore, it is important to highlight that the three modal shapes obtained from the calibrated FEM and the OMA deformed in the same direction (Fig. 14). The first two modal 550·fcp (Kržan et al. 2015) 82-231·fcp (Kaushik et al. 2007) 250-1100·fcp (Paulay and Priestley 1992) 750·fcp shapes have very similar periods due to the square-shaped plan of the tower. Furthermore, the third vibration mode is torsional with a period smaller than the first two. These results are similar to those obtained under several masonry towers in (Valente and Milani 2018). The updated FEM of the tower has been modelled with the boundary conditions and an elasticity modulus, obtained from E = 800·f cp , in the masonry brick.
After the calibration of the FE model had been done, the final mechanical properties used in the updated FE model for each material have been listed in Table 3.

Nonlinear static analysis
The seismic assessment of the tower has been evaluated applying NLSA (Nonlinear Static Analysis). This method has been widely used for the seismic analysis of masonry buildings (Lagomarsino and Cattari 2015a;Vuoto et al. 2022) and it is proposed in the EC8 (European Union 2005a). It is considered an effective method to analyse the global seismic behaviour of existing buildings (Lagomarsino and Cattari 2015a). It is important to properly define several parameters for its execution: the location of the control node, the horizontal load pattern, and the representative displacement in the pushover curve (Lagomarsino and Cattari 2015a). The NLSA has been carried out in the main direction of the tower (X and Y) in positive and negative senses, due to its non-symmetric configuration, which is caused by the ramps and the openings distribution.
In this research, the control node has been set in the centre of the mass in the X and Y coordinates, and on the top of the building (Z-coordinate). Seismic codes propose locating the node on the top floor because it is necessary to place it above the level where the collapse occurs (Lagomarsino and Cattari 2015a). Regarding the horizontal load pattern, several options have been considered, keeping in mind the recommendations of the seismic standards (European Union 2005a) and of the literature (Lagomarsino and Cattari 2015a): proportional to the masses and to the height of the inner levels (pseudo-triangular) and proportional to the masses (uniform).
Firstly, a vertical analysis has been done considering different gravity loads. To that effect, the material loads and the live loads have been applied. Then, horizontal loads have been gradually applied to push the structure up to its collapse. The analysis was performed until the base shear drops by 20%, which usually represents the collapse state of the structure (Castellazzi et al. 2018). The model has been calculated with the parallel option in OpenSees, using multiprocessors. Due to the characteristics of the computer used, the model has been divided into 24 parts. It has been calculated on a computer with an AMD Ryzen 9 5900 × 3.70 GHz processor. Firstly, a mesh-sensitivity analysis has been done. Several mesh sizes have been defined, as mentioned in Sect. 3.2.2. Figure 15 plots the capacity curves in the X-direction with different mesh sizes. The numerical models with mesh sizes between 30 and 60 cm have a stable behaviour. It can be observed that the difference in the results among these is negligible. The 70 and 80 mesh sizes are unstable. Moreover, their capacity has increased proportionally to the mesh size because of its higher dimension. Furthermore, the computation time increases considerably with the refinement of the mesh size (Fig. 15b). In addition, the damage distribution has been analysed and compared. As can be seen in Fig. 15, the damage distribution considering the 40 cm mesh (Fig. 16b) is more accurate than the 60 cm mesh (Fig. 16a). In this regard, the location of the damage in the 60 cm mesh is more irregular, very dispersed damage appearing near the belfry. Due to that, the 40 cm mesh has been adopted for this analysis for its balance between accuracy and computational effort. Furthermore, in the case of the 40 cm mesh, there are, at least, three tetrahedron elements in the thickness of the inner masonry wall. The average calculation time for the NLSA with the 40 cm mesh is 1 h and 9 min.
Secondly, different pushover analyses have been carried out in the X and the Y direction in the positive and negative ranges, applying different load patterns. Also, the tower model isolated and the model with the boundary conditions have been analysed and compared. All the capacity curves obtained in the NLSA analysis have been plotted in Fig. 17.
Significant differences have been found between the different load patterns. The models with a uniform horizontal load have a higher capacity than the pseudo-triangular models. This is due to the fact that the horizontal load is applied uniformly throughout Fig. 15 Sensitivity analysis of mesh size in the numerical model. a Pushover curves, b time required the cane that is more massive than the upper bodies. Furthermore, the resultant of the triangular load is higher. Generally, the capacity is greater in the X direction than in the Y direction. This increase is likely caused by the irregular distribution of the hollows in the outer and the inner walls. In addition, the irregular distribution of the ramps can also affect it. Despite the difference in the maximum strength between the different load patterns, it is worth mentioning that in both load patterns, the maximum deformation capacity is similar (i.e., deformation at the collapse). As can also be observed, the models with boundary conditions have a higher stiffness in the elastic range in the positive analysis. Anyway, the capacity is similar in both models (I and B). Nevertheless, in the case of + X with both load patterns (T and U), the capacity is lower in the case of the  . This may be due to the irregularity in the damage distributions provided by the boundary conditions. However, the capacity is greater in the B models than in the I model in the negative range analysis. Thus, the Cathedral walls provide a higher rigidity in this direction. In these models, the ultimate displacement is higher than in the isolated models.
The different damage patterns obtained from the FEM in the X and Y directions have been analysed and compared. In this vein, in this manuscript the tensile damage models have been plotted because they are more severe than those that are compressive (Figs. 18,19,20 and 21). These damage patterns have been plotted in the ultimate displacement. It is possible to identify the areas of damage concentration. Besides, the damages in both FEM (isolated tower and with boundary conditions) have been compared to identify which is the most accurate numerical model.
Generally, the greatest damage is concentrated near the belfry because it is weaker than the cane of the tower. In addition, the damage is located near the openings in the outer walls as it is its weakest area. The presence of openings is an important issue for local collapse, as demonstrated in (Valente and Milani 2018). The severe and the extensive damage are concentrated mainly near the openings. In that sense, the main geometrical features of the bell tower are the large hollows at the top level. This issue causes considerable damage near the openings of the belfry, suggesting that this is the most vulnerable part of the structure where a local failure could happen. It should be observed that the ramps would be seriously damaged given that they are weak elements between two rigid bodies. Furthermore, there are greater damages near the belfry in the models with boundary conditions. These walls are not symmetric (being only located in the south and west façades). Therefore, they introduce more irregularity into the model. In addition, there is great damage in the contact between these walls and the tower. As can be seen in Figs. 18, 19, 20 and 21, the upper body does not show significant damage. However, in the model with boundary  It is important to highlight that these damage patterns are similar to those described in the testimonies of past earthquakes. In that sense, several testimonies of the 1755 earthquake described that the tower presented the most serious cracks in its upper levels. Besides, a great number of ramps were in a state of ruin.

Conclusions
This manuscript discussed the numerical modelling for the seismic assessment of a complex heritage masonry tower. To do so, the numerical model of the Giralda tower was carried out considering different uncertainties. To this end, a historical analysis and in situ surveys were developed. The FEM was calibrated using the results of the modal analysis and the ambient vibration test. Moreover, a sensitivity analysis of the mesh size in the FEM was done in order to verify the stable behaviour of the numerical model. The numerical model calibrated was used to perform the NLSA.
These types of studies are an opportunity to provide useful information and knowledge about the seismic security of these complex ancient structures. Yet, this useful data is important to be applied in future conservation projects and in the scope of scientific research activity. In this respect, the digital information generated (geometrical and numerical FE model) will be very important for the collection of the data in future works. Additionally, the numerical FE model calibrated can be used in future analysis of the tower. Furthermore, this study showed that advanced numerical analysis can provide valuable information for assessing the damage of existing cultural buildings. This work also aims to It is worth highlighting that the analysis of the uncertainties provides outstanding information for the analysis of its seismic behaviour. Field surveys carried out in this work were important for the characterisation and definition of the numerical model. With that in mind, the ambient vibration test is fundamental for calibrating the numerical model because of the difficulty of carrying out destructive or semi-destructive in situ surveys. The calibration of the numerical model is a crucial step to proceed to the seismic assessment of the ancient building. In addition, it is important to calibrate the mesh size in the FE model due to the complexity of these ancient buildings. Although the energy regularisation of the material is defined, it is essential to check its behaviour according to the mesh size. In addition, the tetrahedral element is better adapted to this type of irregular structures. Furthermore, the results of this study underlined that the best mesh size must be analysed based on stable capacity behaviour, the damage distribution, and the number of elements in the thickness of the masonry wall. It seems advisable to adopt four-node tetrahedron solid elements and a mesh size of 40 cm for this type of structure, and never higher than 60 cm.
It should be noted that, although the nonlinear analysis of complex masonry ancient buildings usually requires a high computational effort, the OpenSees software has been able to perform this analysis in an acceptable computational time. This open-source software significantly reduces the computation time.
The results achieved in this work emphasise the importance of analysing the boundary conditions of ancient buildings to obtain an accurate and updated numerical FE model. It is important to highlight that the boundary conditions have also modified its global dynamic behaviour, as has been observed in the modal analysis (see Fig. 14). It should be mentioned that more precise results have been obtained in the calibration process with the modelling of the equivalent boundary conditions. The frequencies obtained in the modal analysis of the FE model are very close to the in situ measurements (OMA) (difference 2-11.9%). Furthermore, the three modal shapes obtained from the calibrated FEM and the OMA have deformed in the same direction. In addition, the damage patterns of the tower have varied slightly by considering the effects of the perimetral wall of the cathedral. The damage located near the bell tower level has increased minimally (Figs. 18,19,20 and 21). Likewise, the cathedral walls introduce some irregularities in the tower, which should be taken into account in its seismic assessment.
Tensile damage plots have shown that the main geometric features of the bell towers increase its seismic vulnerability. Openings are a weak area in which damage is concentrated. Therefore, the belfry level is the most vulnerable region due to the large hollows that may cause the local collapse mechanism.
The NLSA was used to validate and analyse the global behaviour of the Giralda tower due to its feasibility and limited calculation time. However, future studies should aim to perform a nonlinear dynamic analysis applying real accelerograms to analyse the damage distribution.
Although the Giralda tower of Seville is the highest and most emblematic Islamic tower, there are many other Islamic towers around the word, the most relevant being the Marrakech and the Rabat towers (Jiménez Martín 1993). Therefore, the method developed in this work can be also applied to the analysis of other ancient Islamic towers.