The Engineering Approach to Conservation of Massive Archaeological Structures in Seismic Areas: The Apollo Nymphaeum in Hierapolis of Phrygia

ABSTRACT The safety conditions of archaeological remains in seismic areas are strongly jeopardised by a lack of structural completeness, which may trigger dynamic responses among parts and led to partial or overall collapses. However, some of those structures maintain their massive behaviour, which makes the dynamic characterisation difficult for possible predictive analyses. It is the case of the Nymphaeum of Apollo in Hierapolis of Phrygia (Turkey), an impressive stone block masonry structure of the 3rd c. CE, which lies over an active fault. The incomplete structural conditions and the presence of multi-leaf cross sections make this case study of great interest to verifying the reliability of using numerical models for assessment. The discrete element approach was applied to perform nonlinear dynamic analyses under increasing levels of seismicity. The numerical model was based on: i) the detailed shaping of blocks texture, and ii) the onsite inspections and non-destructive testing. Data was also implemented into the numerical model to perform sensitivity analyses. Results showed that the engineering approach proposed herein is able to overcome the challenges imposed by such a massive archaeological structure, facing the unknowns left over by experimental investigations with the opportunity offered by numerical methods such as discrete element modelling.


Introduction
Archaeological sites represent a fundamental testament to our most ancient history. Their remains are defined as fragile and non-renewable cultural resources that must be protected against the increasing threats worldwide (ICOMOS-ICAHM International Committee on Archaeological Heritage Management 1990). The conservation of archaeological heritage is strictly connected to its structural safety, especially in the case of architectural remains standing in seismic prone areas (Ambraseys 2006). The incompleteness of ruins, in terms of structural systems, combined with the high exposure to environmental deterioration, can increase the weakness of the constructions against dynamic actions (Autiero et al. 2020;Lorenzoni et al. 2017a;Marson et al. 2016;Valluzzi et al. 2019b). In addition, such complex configurations are commonly made of masonry, a composite material often arranged in irregular cross sections with a number of possible variables which affect its mechanical behaviour (e.g., non-linearity and anisotropy, quality of composing materials, connection among portions). Ancient materials also underwent alterations and reconstructions over time, thus making it difficult to generalise methods to evaluate their current conditions and predict their behaviour in light of possible interventions aimed at protecting and conserving them.
Contributions in literature focus on freestanding masonry components, as either single-wall models (Griffith et al. 2003) or multi-layer systems (de Felice 2011), which commonly derive from studies of historic buildings in seismic areas. Archaeological remains often concern massive constructions made of large stone blocks with dry or deteriorated mortar joints. The study of block masonry (Azevedo, Sincraian, and Lemos 2000) led to the progressive refining of the Discrete Element Method (DEM) (Lemos 2007), formerly developed as a distinct element approach to simulate rock movements (Cundall 1971). The DE approach widens the capability of the modelling of masonry, meant as continuous and homogeneous material by the Finite Element Method (FEM) (Giamundo et al. 2014;Giordano, Mele, and De Luca 2002;Lee et al. 1996;Mendes et al. 2017). DEM models masonry as a discontinuous assembly of elements interacting along nonlinear borders and has been increasingly applied to failure analyses of large components under either static or dynamic/seismic loadings, from single rigid-block structures (Peña et al. 2007) to entire masonry arches (De Lorenzis, and Ochsendorf 2007;Tóth, Orbán, and Bagi 2009) and vaults Gobbin, de Felice, and Lemos 2020;Mcinerney and Dejong 2015). However, the focus on masonry walls behaviour shows scarce contributions, which concern the modelling of solid (i.e., single-leaf only) masonry, whether in-plane (Malomo, DeJong, and Penna 2019), out-of-plane (Lemos and Campos Costa 2016;Meriggi et al. 2019), or both (Bui et al. 2017).
The application of DEM to archaeological ruins mostly refers to either multi-drum columns, composed colonnades of ancient temples and aqueduct arcades (Drei and Oliveira 2001;Mordanova and de Felice 2020;Oliveira et al. 2012;Papantonopoulos et al. 2002;Pappas, da Porto, and Modena 2016;Psycharis, Fragiadakis, and Stefanou 2013;Psycharis et al. 2003;Sarhosis et al. 2016). It seldom refers to more articulated sub-structures. Psycharis, Drougas, and Dasiou (2011) analysed the seismic behaviour of some walls of the Parthenon to evaluate the effect of restoration works; Young, Schultz, and Lemos (2015) estimated high magnitude earthquakes to make the remains of a large-sized temple in Nemea (Greece) collapse. Recently, in Mordanova and de Felice (2020) a 2D analysis focused on the causes of damage of the restored massive wall of the Colosseum in Rome; in Saygili and Lemos (2020), DEM was used to identify the dynamic behaviour of the Frontinus Gate in the archeological site of Hierapolis (the same area of this study) and predict its damage under several seismic inputs.
These last few experiences demonstrate the current difficulties in comprehensive approaches towards complex structures and the need for more cases to investigate the reliability of specific strategies applied to archaeological monumental sites. However, despite the increasing capability of modelling tools, what essentially is missing, also in recent applications, is the availability of proper strategies to make models representative of the real structural behaviour. This must be based on direct information collected onsite by the experimental diagnostics through minor or non-destructive investigations (European Committee  for Standardization 2005; ICOMOS-ISCARSAH 2003; International Organization for Standardization 2012), so that the main properties of the assemblages can be identified and used for model generation and validation (Clementi et al. 2017;Elyamani et al. 2017;Kim et al. 2021). In addition, the results from assessment of a historical structure need to be crosschecked with the real structural conditions (the so-called 'plausibility check', according to International Organization for Standardization (2012)), which can involve the effects of either deteriorated conditions or residual stability capacity. This particularly applies to the challenging field of conservation of archaeological remains and is eventually integrated in the framework of a more comprehensive approach herein proposed.

Research aim
This study aims to propose an integrated procedure that can improve the representativeness of numerical models for the seismic assessment of archaeological masonry remains with particular focus on massive structures. This method was applied to the Nymphaeum of Apollo's sanctuary in Hierapolis of Phrygia (Denizli region, Turkey), based on a multidisciplinary knowledge path, which involved a research team of structural engineers, geophysicists, archaeologists and historians (Lorenzoni et al. 2020;Marson et al. 2019;Ricci et al. 2020). The site is noteworthy for understanding how relics of ancient monuments can still undergo past and possible future archaeoseismic damage (Hancock and Altunel 1997). The onsite collection of actual structural parameters and the definition of missing data through sensitivity analyses allowed numerical 3D DE simulations to be carried out, and predictive analyses and damage scenarios to be evaluated.

Material and methods
The engineering approach proposed here is based on an experimental-numerical integrated procedure which can identify reliable parameters for material characterization and to validate representative models for structural simulations. Direct (e.g., from onsite surveys) and indirect (e.g., hailing from seismic codes, knowledge available in literature) data was combined with hypotheses and assumptions on unknowns, which must be confirmed on onsite investigational basis and/or compared with sensitivity analyses on pre-validated models.
In particular, ambient vibration testing for the structure, as well as geophysical investigations for the soil, prove essential in this context. The dynamic characterisation can be used not only for model updating and calibration but also to understand whether the building shows (or not) an overall behaviour, thus guiding the choice of the suitable modelling strategy (e.g., continuous vs. discrete approach). Finally, data collected from geophysical investigations defines soil properties, identifies possible site effects and allows the adjustment of the input for seismic analyses. Figure 1 shows the layout of the integrated approach to assess the seismic vulnerability of archaeological structures. The identification of the 'plausible' model of the Nymphaeum discussed in this paper follows the elements connected by the continuum branches in Figure 1, according to the plausibility check requirement (International Organization for Standardization 2012).

Experimental campaign and preliminary data processing
Architectural heritage is recognised as a challenging legacy to protect, due to limitations both in diagnosis and the application of modern safety codes (ICOMOS-ISCARSAH 2003; International Organization for Standardization 2012). An accurate investigation plan based on non-destructive experimental procedures is mandatory to support either the assessment or the possible intervention phases (Binda, Saisi, and Tiraboschi ). The knowledge path applied to the Nymphaeum aimed at characterising both the structure (whose mechanical behaviour is strictly based on the geometrical configuration of the block masonry) and the soil, to define the parameters for the numerical model implementation.

Structure description
The Nymphaeum of Apollo's sanctuary has scarce sources about its history. It dates back to the beginning of the 3 rd c. CE and was probably erected on the site of an older fountain, based on the investigations performed between 1958 and 1971 (Campagna and Scardozzi 2013). The Nymphaeum is located approximately in the middle of the western edge of the lower terrace of the sanctuary (D'Andria, Caggia, and Ismaelli 2016; Semeraro 2012), in front of the main buildings, and with the front facing west towards the theatron and the underlying terrace Spanò et al. 2018).
The structure is today an imposing C-shaped relic (Figure 2 The foundation layer is made of small stones and fragments of bricks bonded by lime mortar (see also Figure 2(d)); two courses of travertine blocks form the substructure's base. Multi-leaf masonry walls compose the superstructure, whose outer leaves are made of dry blocks of travertine, typical of the Denizli region. In the rebuilt and last original layers, blocks have average width (on facade) of 78 cm or 95 cm, respectively. Visual inspections on accessible elements of the superstructure (mainly at top and near niches) gave data on blocks' depth ( Table 1). The highest values of the coefficient of variation (σ) could suggest the presence of longer elements, whose effectiveness, however, cannot be proven.
The structure presents widespread conditions of scarce conservation; the south wall is the most damaged, with partial collapse of blocks, and leaves separation and tilting of the inner leaf (about 5 degrees) (Figure 3(b)). However, this damage made the wall cross-section observable: the walls' top has leaves highly disconnected from the filling, which appears thicker than the outer strata (about 1.5 m on wings and less than 3 m on the main body) and as a layer of reused stone units, made incoherent by damage and binder deterioration due to environmental phenomena. Other inspections revealed instead a core composed of a sort of conglomerate of mortar and reused fragments and units from neighbouring structures: spot from a niche in the north wing is reported in Figure 3(c), while infill top face on the north wing-main body corner is shown in Figure 3 (d). Therefore, a unique definition of the composition of the infill was not possible, and its presence along the entire structure, as well as the quality of bonding to the block leaves, can be only supposed based on data collected on site. The most probable layout concerns a coherent infill for the so-called 'last original phase', locally damaged within the south wing. This revealed to be a key factor for the model reliability. Figure 3(e) shows the main dimensions of the Nymphaeum.

Geophysical investigations.
A seismic survey was carried out to characterise the soil shear wave velocities (V S ). Twenty-four 4.5 Hz vertical geophones, placed 1.5 m apart each other, were installed (total survey length of 34.5 m) and a 5 kg hammer was used as seismic source; data was acquired with 0.25 ms sampling. The multichannel analysis of the surface waves technique (Park, Miller, and Xia 1999) was applied to evaluate the dispersion of the surface waves and retrieve the velocity model of the subsoil.
Altered rock was identified in the first 20 m depth with a harder bedrock corresponding to higher velocities in the deeper part. V S ranged from 338 m/s at the surface to 800 m/s in the deepest layer. According to Eurocode 8 (European Committee for Standardization 2005), a reference value of V S at a depth of 30 m V S,30 = 528 m/ s was found, which corresponds to a soil B class (i.e., a deposit of very dense gravel or very stiff clay at least several tens of metres in thickness characterised by a gradual increase of mechanical properties with depth).

Dynamic characterisation and modal analysis
Ambient vibration tests were carried out as nondestructive method to identify the dynamic response of the remains (Ramos et al. 2010; Ubertini, Comanducci, and Cavalagli 2016; Lorenzoni et al. 2017b). The measurement and analysis of vibrations induced by environmental noise (e.g., wind, traffic, human activities, micro-tremors) allows modal parameters such as natural frequencies, damping ratios and mode shapes to be identified. This procedure has been extensively validated (Foti et al. 2012) to detect damage (Mosavi et al. 2012), assess vulnerabilities (Gentile and Saisi 2007) and calibrate behavioural models for structural assessment and seismic analyses (Aras et al. 2011;Lorenzoni, Valluzzi, and Modena 2019).
The equipment (National Instruments PXI) consisted in 13 piezoelectric seismic accelerometers (PCB 393B12 type) connected to a data acquisition unit that recorded vibrations and response signals in the time domain. Sensors were arranged in one setup according to the expected mode shapes and installed on the corners, the edge of the north wing and the middle of the other walls ( Figure 4). Multiple sets of data were recorded from accelerometers (sampling frequency of 100 Hz) with different operational conditions: i) ambient noise produced by wind and human actions (e.g., trucks and crane operations nearby) and ii) random dynamic impulses in space and time using an instrumented hammer and mechanical impacts.
Three structural modes were extracted through output-only modal analysis techniques in accordance with the experimental measurements. The signals were processed by ArTeMIS software (Structural Vibration Solutions 2006) through two identification methods, i.e., the enhanced frequency domain decomposition (EFDD, Brincker, Ventura, and Andersen 2001) and the stochastic subspace identification principal component (SSI-PC, Peeters and De Roeck 1999). The structural modes were identified through the peak-picking technique from the singular values of spectral densities ( Figure 5a); then, the stable poles corresponding to the structural modes were selected from the stabilisation diagram calculated by the SSI-PC algorithm ( Figure 5(b)). The first three mode shapes referred to the frequency range 7-13 Hz ( Figure 5(c-e)); Table 2 reports their dynamic parameters.
No overall dynamic response was identified, due to the massive and stiff three-leaf walls: out-of-plane mode shapes activated for the southeast walls (first and second modes) and the north wing (third one). Several local modes were also extracted in the south wing, due to the high level of damage and disconnection. These outcomes also suggested that the use of a continuum model to represent the structural behaviour of the Nymphaeum is not appropriate; the use of a discontinuous approach could instead provide more reliable results.

Numerical modelling strategies
According to the modal analysis, the Nymphaeum has predominant local modes. This behaviour is typical among archaeological remains or heavily damaged structures and can be adequately described through the geometrical equilibrium theory (Heyman 1966;Housner 1963). The interlocking of blocks achieved through the particular construction technique of such structures leads to displacements controlled mainly by rocking and sliding modes. Since the discrete element method (DEM) accounts for both, it is suitable for carrying out predictive analyses on the seismic response of structures of this nature.
The model of the Nymphaeum was built in 3DEC 5.0 (Itasca Consulting Group Inc 2013) a 3D DE code which has been proposed since the 1980s  (Cundall 1988;Hart, Cundall, and Lemos 1988). Nonlinear dynamic analyses were carried out, implementing velocity time-histories (time-history analysis, THA) compatible with the onsite geophysical investigations ( §2.1.2). In particular, the geometry of the structure was constructed by site measurements and material properties were deduced by analytical estimations and literature references. Common model updating techniques (i.e., adjustment of model mass and stiffness to match target mode shapes) could not be applied here, owing to the scarce conservation status. Furthermore, sensitivity analyses were carried out to crosscheck the assumptions made, as described in the following and discussed in ( §3.1).   Figure 6. Outline of material groups implemented in DE model of Nymphaeum (a); position of control points (S1 to S3) for dynamic THAs (b).

Geometry and preliminary material properties
The geometry of the numerical model represented the current shape and conditions of the Nymphaeum, described by five groups of elements (Figure 6(a)): i) substructure (base and foundation); ii) last original phase (low superstructure); iii) rebuilt phase (top superstructure); iv) south wing (heavily damaged condition); v) infill. Figure 6(b) shows the control points assumed for the dynamic analyses of the numerical models (see §3.2). The Nymphaeum was discretised with rigid blocks connected by deformable interfaces ruled by Mohr-Coulomb law (Coulomb 1776;Mohr 1900), which lump both linearities and nonlinearities. Mechanical properties of stone blocks taken from literature were used for the calculation of joint stiffnesses, as follows.
Interfaces among blocks are defined by joint stiffness, j k , friction angle, φ, cohesion, c, and tensile strength cutoff limit, f t , and were calculated (groups refer to Figure 6(a)): -for group i) and v), based on the linear properties of the homogenised masonry, where contribution of both units and mortar were taken into account; -for groups ii)-to-iv), based on the properties of travertine disjointed blocks.
The linear properties of joints were calculated according to the analytical approach suggested by (DeJong 2009): where j kn and j ks are the normal and the shear joint stiffnesses, E and G are the Young's and the shear moduli of masonry, and L is the length of rigid material in the direction perpendicular to the joint (i.e., layers for the substructure, the single blocks for the superstructure).
Based on the onsite survey ( §2.1.1), the wall leaves of the superstructure were modelled as squared blocks (0.6 m thick, as average value from Table 1); existing damage (e.g., gaps in block arrangements, niches) and uneven profiles (i.e., deviations among blocks, tilt of the inner portion of the south wing) were also included. The unit weight of travertine blocks resulted from the geomechanical data of the Denizli region, according to studies on neighbouring structures (Yagiz 2006). A Young's modulus of 12000 MPa was chosen for the group ii) (last original phase), based on travertines from the Antalya region (Yagiz and Akyol 2005). Values of the reconstructed layer of the main body and south wing were slightly decreased (of 1000 and 2000 MPa, respectively) to take into consideration their status of scarce conservation and damage.
The substructure (group i) in Figure 6(a)) was modelled through three superimposed layers of rigid blocks forming a C-shaped support, whereas a larger block underneath was used to bear the model and apply the seismic input.
In the absence of specific values, E, G and γ of groups i) and v) were taken from the Italian seismic code (Ministero delle Infrastrutture e dei Trasporti 2019), which correspond to squared stone block masonry and irregular masonry (made of pebbles and ashlars), respectively, for the substructure and the infill. Table 3 lists the mechanical and geometrical properties for the five groups of materials (the size of the infill blocks depends on the discretisation assumed and underwent sensitivity analyses, as discussed in §3.1.1). Table 4 shows the linear properties of the interface joints calculated according to Equations (1) and (2) (values of blocks of a portion of the superstructure were also crosschecked by sensitivity analyses, see §3.1.2). Cohesion, dilatation angle and tensile strength were set to zero, while a typical value of 37° of the friction angle was considered, according to (Papantonopoulos et al. 2002;Psycharis et al. 2003;Psycharis, Papastamatiou, and Alexandris 2000;Pulatsu et al. 2017).

Seismic input
Hierapolis is located in the highest seismic hazard zone among the five subclasses of Turkey (Kayabali and Akin 2003; Ministry of Public Works and Settlement Government of Republic of Turkey 2007). Moment magnitudes Mw 6.0 to 7.0 (Richter scale) represent the range of seismic potential in the region, and Mw 6.3 is recognised as the reference for the most probable earthquake (Erdik et al. 1999;Inel et al. 2008); expected values for peak ground acceleration (PGA) of the response spectrum range from 0.2 g to 0.4 g. The target horizontal spectrum was made using FEMA-368 guidelines (Building Seismic Safety Council 2001) for soil class B ( §2.1.2): it corresponds to a PGA = 0.35 g, an average return period Tr = 475 years, a damping of 5%, and minimum and maximum distances from the epicentre of 0 and 25 km in the range of the expected magnitudes.
The modelling of motion was defined by seven natural accelerograms selected with one component from the European strong motion database, fitting the target spectrum for a soil B site and then scaled to the expected PGA of 0.35 g (Inel et al. 2008;Kayabali and Akin 2003) (a tolerance range of 30%-10% was used in a period range T = 0.15s-2s) (Figure 7(a)).
According to Eurocode 8 (European Committee for Standardization 2013) and the Italian code (Ministero delle Infrastrutture e dei Trasporti 2018), time-histories should be applied along both horizontal directions (and in some cases which, however, do not apply to the case of the Nymphaeum, also in the vertical one). Since the preferential earthquake direction for the site is unknown, seven couples of two time-histories (among the seven generated) with same PGA were built for the X and Y horizontal directions of the model. The couples of horizontal time-histories were scaled (per each target PGA) according to an SRSS (Square Root of the Sum of the Square) spectrum defined as per the following equation (3): where Sa i,x (T) and Sa i,y (T) are the unscaled spectra in the X and Y directions, respectively, PGA i,x and PGA i,y are the PGAs of records in X and Y, Sa target (T) is the target spectrum and PGA target the target PGA, α is a corrective factor (ranging from 1 to ffi ffi ffi 2 p , as suggested by Ministero delle Infrastrutture e dei Trasporti (2019)). In this study, the conservative value of ffi ffi ffi 2 p was assumed (case of component amplitudes equal each other).
Based on the seven accelerogram couples, corresponding velocigram ones were extracted (SeismoSoft 2016), to be used as input for the dynamic analyses. The whole set of velocigrams (TH1 to TH7) and their combination in the seven couples (A1 to A7) are reported in the supplemental material (SM) adjoined to this paper. Figure 7b shows the velocigram n. 4 used to perform sensitivity analyses ( §3.1), as it was the signal with the least noise. Zero damping was deemed the best approach for the conservative analysis of the Nymphaeum, as this provides the best matching results during the strong motion of an earthquake (Papantonopoulos et al. 2002;Psycharis, Papastamatiou, and Alexandris 2000;Salvalaggio et al. 2021). Energy dissipation was based on frictional sliding only.

Optimisation of archaeological masonry properties
The model of the Nymphaeum was preliminarily optimised to evaluate the influence of uncertain data on the overall dynamic behaviour of the structure. Two sensitivity analyses were performed, focusing on: a) the infill discretisation (three models) and, b) the joint stiffness evaluation at interfaces among leaf blocks (four models). The latter served as a further validation of the calculated values of Table 4, focusing the analysis on a wall portion of the last original phase of the main body.

Infill
Three models were examined with respect to the infill discretisation: i) absence of infill; ii) infill modelled as random blocks, comparable to the size of the leaves (called 'small' blocks henceforth); iii) infill modelled as large blocks (dimensions at least 10 times higher than the leaves). To evaluate the predominant out-of-plane response of the main body wall and of the two wings of the Nymphaeum the reference velocigram (Figure 7(b)) was applied as seismic input in both the X and Y directions for 25 s. Figures 8(a-c) show the state of the Nymphaeum at the final step of each nonlinear dynamic analysis. In all cases, local collapse of the damaged part of the south wing occurred. In the i) model, where infill was not modelled, an additional collapse of the biggest part of the outer leaf of the north wing took place (Figure 8(a)).
In the ii) model, the discretisation of the infill as 'small' blocks induced the collapse of most part of the main body inner leaf (Figure 8(b)); this was due to the higher deformability of the infill joints compared to those of the outer leaves (Table 4) and to the non-restraining infill-to -leaves connections. In the iii) model, no additional damage was observed (Figure 8(c)).
The out-of-plane displacements of eight control points (Figure 8(d)) were drawn for the three models. As a representative example, Figure 8(e) shows the response at point n. 2 on the north wing. Rocking occurred for all the three models.
Infill discretisation has shown to be a remarkable parameter affecting the dynamic behaviour of the Nymphaeum. Based on the comparison of the damage simulated by the three models ( Figure 8) with that of the real structure (Figures 2 and 3), which identifies the south wing as the only part clearly damaged, the model iii) (i.e., discretisation of infill through large blocks) was assumed and implemented in the numerical analyses of the Nymphaeum. Such a choice was also confirmed by the outcomes of the onsite geometrical survey and the inspections of the infill, which seems to be a coherent mixture of ashlars kept together by a large amount of mortar (Figure 3(c-d)).
The resulting 'plausible' model was composed of 993 blocks, with each one having six degrees of freedom (three translations and three rotations). This model was used to refine the assumptions on block joints and to assess the seismic vulnerability of the structure.

Block joints
A further sensitivity analysis was carried out on a portion of the main body superstructure, aimed at evaluating the contribution of joint stiffness on the behaviour of the structure to possibly correct the interface properties outlined in Table 4.
The reference velocigram was applied in the Y direction (i.e. the out-of-plane of the main body) for the strongest motion interval of the record (i.e., between 10 and 20 s, Figure 7b). The Y-displacements of four control points at the top of the wall were recorded (Figure 9(a)). Three block groups with different material properties were considered: a) infill, b) leaves of the last original phase and c) substructure. The properties of the infill and substructure were those calculated analytically in §2.2.1 and were kept constant in the sensitivity analysis. The normal and shear stiffnesses of the joints of the leaves of the last original phase varied in four configurations: the first one derived from the analytical calculation (Table 4) and the others were fixed (both j kn and j ks ) as 10 9 Pa/m, 10 10 Pa/m and 10 11 Pa/m, according to typical values found in literature (Lemos, Campos Costa, and Bretas 2011;Psycharis, Drougas, and Dasiou 2011;Psycharis, Fragiadakis, and Stefanou 2013;Psycharis, Papastamatiou, and Alexandris 2000).
The behaviour of all control points was similar; as an example, Figure 9(b) shows the curves for point n. 1. The only case where collapse occurred referred to the lowest value of 10 9 Pa/m. Looking at all cases, the wall underwent a maximum displacement at the PGA = 0.32 g (the closest to the expected one in the site) and then returned to the original position, with some residual displacements. It has therefore been concluded that the values calculated analytically (Table 4) are suitable for use. The increase in the joint stiffnesses does not significantly affect the rocking behaviour, whereas the further decrease can cause the inner leaf to collapse; this is considered improbable (plausibility check) since that portion of the building faced some significant seismic events, and does not show a significant damage pattern.

Dynamic analyses of 'plausible' model
Nonlinear THAs were performed to evaluate the structural response under seismic actions, using the representative horizontal accelerograms of the ground motion ( §2.2.2). Analyses were carried out according to the Eurocode 8 (European Committee for Standardization 2005), by combining the seven 25-s velocigrams compatible to the ground conditions of the region in the X and Y directions. For the three walls (main body and two wings) of the Nymphaeum,  the leaves and the infill were checked separately, by setting control points along and through three cross sections (S1 to S3, Figure 6(b)).
Sets of analyses were run at PGA = 0.32 g, in accordance with the target horizontal spectrum for Hierapolis ( §2.2.2), and at PGAs of 0.4 g and 0.5 g (the available accelerograms and velocigrams were scaled accordingly). Figure 10 and Table 5 summarise the state of the Nymphaeum at the end of each analysis for all PGA values in terms of collapse (collapse or noncollapse status is related to the structure leaves).
At PGA = 0.32 g, damage was limited to the collapse of part of the inner leaf of the south wing, which is currently damaged (Figure 10(a)). This confirms the reliability of the model in representing the actual state of damage: with the exception of the south wing, no other damage occurred due to the expected seismic input. Α higher level of damage occurred at PGA = 0.4 g, with additional overturning of the inner leaves of the mean body and the north wing ( Figure 10 (b)). The level of damage at PGA = 0.5 g was significantly higher compared to the two previous sets of analyses: the damaged part of the south inner leaf collapsed, as well as both leaves of the mean body and the north wing. In all cases, corners of the structure were still intact, proving they are the least vulnerable parts compared to the long slender walls between them (Figure 10(c)). The ultimate damage scenarios of the entire set of analyses are reported in SM.
Displacement versus time curves showed rocking and sliding behaviour throughout the analyses at all PGAs (i.e., the collapse of a few blocks due to sliding, or of entire leaves due to overturning). Selected diagrams among the most significant results are reported in SM. Those obtained for PGA = 0.32 g are discussed here. The rocking  Figure 11. Results of nonlinear THAs of Nymphaeum for PGA = 0.32 g: out-of-plane displacement vs. time (S1 top inner leaf) for all analyses of north wing (a); relative displacement of control points along height of north wing (S1 outer leaf) at end of analyses (b) and of main wall (S2) at end of A7 (c).
magnitude of control points varied between 5 cm and 10 cm, as well as residual displacements (Figure 11(a-b)). In all cases, top blocks showed the most significant deviation compared to the substructure, where relative displacement was minimal (Figure 11(b-c)). The relative displacement along the height of the walls was also plotted for the leaves and the infill on the same graph for all analyses, as an indicator of the overall behaviour of the walls. In general, the infill showed minimal displacements compared to the leaves. However, the sliding of blocks after the end of the analyses caused a significant deviation between the top two control points (e.g., S2 for A7 in Figure 11(c)).

Conclusions
The discrete element approach proved a reliable method of representing complex archaeological structures made of massive masonry, as in the case of the Nymphaeum of Apollo in Hierapolis of Phrygia (Turkey). DEM was able to describe the non-unitary behaviour resulting from the discontinuous features of the structure, and to evaluate the uncertain or unknown parameters through sensitive analyses, the results of which were checked with onsite observations. The 'plausible' model of the Nymphaeum underwent nonlinear dynamic analyses, performed through seven velocigrams sets per three PGA values (0.32 g, 0.4 g, and 0.5 g). Results showed that the structure can withstand the expected earthquake without losing global integrity, with only minor local collapse of few blocks of the inner leaf of the south wing, which is seriously damaged today. The rest of the structure can bear the seismic action without damage, other than residual displacements of the blocks, which become more significant at the top sections of the structure, as can be expected. Indeed, this part, for both scarce conservation and seismic exposure, remains the most vulnerable of the Nymphaeum, and should consequently be protected, e.g., by increasing the connection to infill and standing walls.
Numerical modelling provides increasingly advanced tools to simulate extensive scenarios of structural behaviour for existing heritage in seismic areas. This study contributed to the clarification that experimental nondestructive testing (e.g., with ambient vibration) could be used not only within model updating techniques, which are more suitable for more recent or, at least, better preserved structures, but also to guide the model strategy. Furthermore, the reliability of models must be checked on the basis of onsite investigations and surveys, to make analyses representative of the current behaviour. Further applications of this engineering approach will give useful directions for more comprehensive models to be applied to larger assemblages and more complex constructions.