Investigation of the Structural Dynamic Behavior of the Frontinus Gate

Abstract

The Western Anatolia Region of Turkey is an important region of high seismic activity. The active dynamics of the region are shaped by a compression and expansion mechanism. This active mechanism is still ongoing and causes strong seismic activity in the region. The Frontinus Gate is a monument in the Roman city of Hierapolis of Phrygia located in southwestern Anatolia. The aim of this study is to investigate the seismic behavior of this stone masonry structure using discrete element modeling. For this purpose, nonlinear dynamic analyses were performed to simulate the structural response of the gate under seismic excitation. Deformation, damage, and failure patterns induced in the masonry gate for different levels of seismic action are evaluated and discussed. An earthquake with a return period of 475 years is expected to cause some damage, but no collapse, while for a return period of 2475 years, the models indicate collapse of the monument.

1. Introduction

The Hellenistic and Roman city of Hierapolis in Phrygia is one of the most important archaeological sites known from the Roman world. Hierapolis of Phrygia is located in southwestern Anatolia about 200 km east of Izmir and the Aegean coast. The urban area stretches over a travertine shelf looking onto the broad and fertile valley of the Çürüksu River, the ancient Lykos [1,2,3,4]. The city was probably established by Eumenes II of Pergamum in 190 BC [1,2,3]. In Roman times, Hierapolis was an important Asian city characterized by many buildings in travertine and marble; during the Early Byzantine period the city became the metropolis of Phrygia, a very important site due to the presence of the tomb of the Apostle St Philip [5]. The region has a high seismic activity due to a fault line running right below the city for more than a kilometer, applying a NE-SW horizontal stretching to the area [6]. In the first century BC, several earthquakes are known to have caused heavy damage. However, the city was rebuilt during the reign of the Roman emperor Tiberius in AD 14–37. After a major earthquake in the middle of the seventh century AD, the city was severely damaged, but it survived until a new major earthquake in 1354. During the following centuries, Ottoman houses and small farms were constructed in the ruins of the urban area [5]. Between 1653 and 1889, the region faced six devastating earthquakes which caused serious damage [7,8,9,10].

Under the auspices of the Flavian proconsul Sextus Iulius Frontinus, a monumental gate was built at the northern edge of the city, in the first century AD, and the adjacent part of the main street was extended and framed by a Doric colonnade [11]. The gate has three openings in squared travertine blocks, with masonry arches, flanked by two round towers. Under the directorship of Paolo Verzone from the Turin Polytechnic Institute, the Missione Archeologica Italiana (MAIER) was set up in 1957, and restorations on the Frontinus Gate progressed [12]. Archaeological surveys have been carried out in the urban area and in the territory surrounding the ancient city to collect information [13]. The Frontinus Gate, which is the monumental entrance to Hierapolis (Figure 1), suffered damage from earthquakes in the fourth and seventh centuries, after which some reconstruction interventions were applied [12,14]. One of the towers of the opening is still in good condition today. Almost half of the other tower collapsed and was damaged from destructive earthquakes and other natural disasters. Since 1988, the site has been included in the UNESCO World Heritage List with emphasis on the extraordinary natural conditions, the Greco-Roman thermal installations, and the Christian monuments [15]; thus the Frontinus Gate and other monuments at the site have great importance in providing the present generation a connection to its past.

The aim of this study is to investigate the dynamic characteristics of the Frontinus Gate and simulate the seismic behavior of the structure under the seismic excitations that can be expected at the site, in order to predict the potential damage levels. The tool employed is the discrete element method (DEM), a numerical technique increasingly applied to the analysis of masonry, in particular of historical heritage structures [16]. The analysis of masonry structures may be performed either by equivalent continuum models or by discrete block models [17]. The former, based on the material homogenization approach, are typically preferred for large complex structures, while the discrete (or discontinuum) methods, where DEM is included, were initially applied to components or simple structures, but are presently capable of addressing more complex systems [18]. DEM represents the structure as a system of discrete blocks, thus closely reproducing the physical nature of masonry. It is capable of simulating the nonlinear phenomena of slip and separation along the joints, which are typically observed in masonry structures under seismic loads, allowing the analysis to proceed into the large displacement range to follow the processes of damage and progressive collapse [19]. DEM models have been frequently applied to tall, slender structures, from the early work on the Parthenon columns [20], to the more recent analyses of towers [21], obelisks [22], and stone minarets [23]. More complex geometries are now being addressed, such as arch bridges [24], domes [25], various types of buildings [26,27], or archaeological structures [28]. Comparisons of DEM results with the behavior observed in shaking table tests of scaled models have provided a progressive validation of its performance. Experiments with marble drum columns provided the earlier assessment of DEM representations [29]. Extensive shaking table tests of single stone blocks rocking under harmonic and seismic records have confirmed the good performance of the numerical method [30]. Comparisons with tests of more complex structures have also been undertaken, namely masonry walls under out-of-plane loads [31,32], or the scaled model of a mosque [33]. DEM has been applied to many historical structures, namely the Parthenon columns [20,29], the Colosseum walls [26], a Roman temple [27], an obelisk [22], historical minarets [23], or the dome of Florence cathedral [25]. In the present study, the effect of the wall curvature on the stability of the round towers is an issue to be examined. The DEM models presented herein are intended not only for safety assessment under large seismic actions, but also as a means to quantify the damage to the structure that can be expected from lower intensity, more likely events.

2. Seismicity of the Region

The present tectonic structure of Turkey has been formed by active interactions developed by the continental collision of the Arabian plate with the Eurasian plate in the east, and subduction of the African plate beneath the Aegean [34,35]. These tectonic motions caused the development of two significant active faults, the North Anatolian Fault (NAF) and the East Anatolian Fault (EAF), and the westward extrusion of the Anatolian plate. Another result of the motion is the ongoing lithospheric scale extension caused by trench roll-back in the Hellenic subduction zone [36,37]. As a consequence of the active tectonic motions, strong extensional deformation developed in western Anatolia.

Before instrumental recording (1900), only earthquakes having magnitudes greater than 4.5 were reported. One of the reported destructive earthquakes occurred in Hierapolis and Phrygia in 494 BC. Then, another devastating earthquake occurred in Laodicea and Hierapolis in 60 BC. Both events caused heavy damage in the regions. Throughout history, similar destructive earthquakes have occurred in this region and surrounding areas. Important historical earthquakes with magnitudes between Mw 6.2 and Mw 6.6 heavily affected Laodicea and Hierapolis. In 26 BC an earthquake with a magnitude of Mw 6.2 completely destroyed the ancient city. In 17 BC an earthquake with a magnitude of Mw 6.6 resulted in damages in 12 ancient cities, causing surface ruptures and landslides.

Focusing on Denizli’s geologic and tectonic structure (Figure 2), the N-S extension system in Western Turkey and the Africa plate subducting under the Anatolian plate strongly affect the Denizli region and surrounding areas. The most active grabens in the province are the NW-SE extending Gediz graben and the E-W extending Menderes graben (Figure 2). The Denizli basin is situated close to the confluence of these two grabens.

As seen in Figure 3, the magnitude range of most of the earthquakes in this region is from 2 to 4. This indicates that active crustal movements and numerous outflow thermal waters are responsible for high micro-seismic activity. Denizli is one of the most seismic regions in Turkey. Figure 4 shows the number of earthquakes with depth distribution for Denizli. As mentioned, Denizli was destroyed by an earthquake with a magnitude of 4.7 on 19 August 1976. Although it was a moderate earthquake, it resulted in the collapse of 40 houses and heavy damage to 1284 houses [38].

3. Numerical Model

The Frontinus Gate has three openings in squared travertine blocks, with masonry arches, flanked by two round towers. Some of the blocks above the arches collapsed and/or were damaged. One of the round towers of the opening is still in good condition today, standing at a maximum height of 8.17 m. The external and internal diameter of this tower are 10.20 m and 9.09 m, respectively. Almost half of the other tower has collapsed. This damaged tower is 3.18 m at its highest point, with 8.25 m and 7.11 m external and internal diameters, respectively. Two of the three masonry segmental arched doorways have almost the same dimensions. The intrados are circular but less than a semicircle with 4.58 m above ground level [1,2,4,5]. In the literature, there are several archeological, experimental and numerical studies investigating the material characteristics of the stone masonry structures in Hierapolis [8,9,39,40,41,42]. Most of these studies emphasize the material properties of travertine stone block masonry and ashlar (regular) masonry based on destructive and nondestructive testing procedures. The Frontinus Gate has traces of damage to the stones at round towers and arches from earthquakes or from other natural disasters.

The average mechanical properties of the Frontinus Gate were determined using information inferred from previous studies mentioned above and considered obvious damage. In this study, elasticity and shear modulus were assumed as 2800 MPa and 860 MPa, respectively. The numerical model was created with the discrete element method code 3DEC [43], widely used in masonry models [20,21,22,23,24,25,26,27,28,29,30,31,32,33]. 3DEC represents the discontinuous medium as an assemblage of discrete blocks, with the discontinuities treated as interfaces between the blocks. 3DEC allows large displacements and rotations of blocks, thus being able to simulate structural collapse. In this study, rigid blocks have been employed, as often done in the analysis of structures made of strong stone units, for which deformation and failure modes are governed by sliding or separation along discontinuities [19].

The solution is based on the integration of the equations of motion of the blocks, using an explicit time-stepping algorithm. The 3 translational equations of motion of a rigid block center of mass may be expressed as:

$$\mathrm{m}{\ddot{\mathrm{u}}}_{\mathrm{i}}+\mathsf{\alpha }\mathrm{m}{\dot{\mathrm{u}}}_{\mathrm{i}}={\mathrm{f}}_{\mathrm{i}}^{}$$

(1)

where u_i denotes the displacement vector of the block center; m, the block mass; and α, the mass-proportional viscous damping parameter, which reproduces the energy losses in the system beyond frictional dissipation. The force vector f_i is the sum of the applied forces, including self-weight, and the contact forces, which are a function of the relative block displacements, and therefore includes the elastic forces. The 3 rotational degrees of freedom are governed by Euler’s equations, expressed in the principal axes of inertia of the block as:

$$\begin{array}{l}{\mathrm{I}}_1{\dot{\mathsf{\omega }}}_1+\alpha {\mathrm{I}}_1{\mathsf{\omega }}_1+\left({\mathrm{I}}_3-{\mathrm{I}}_2\right){\mathsf{\omega }}_3{\mathsf{\omega }}_2^{}={\mathrm{m}}_1^{}\\ {\mathrm{I}}_2{\dot{\mathsf{\omega }}}_2+\alpha {\mathrm{I}}_2{\mathsf{\omega }}_2+\left({\mathrm{I}}_1-{\mathrm{I}}_3\right){\mathsf{\omega }}_1{\mathsf{\omega }}_3^{}={\mathrm{m}}_2^{}\\ {\mathrm{I}}_3{\dot{\mathsf{\omega }}}_3+\alpha {\mathrm{I}}_3{\mathsf{\omega }}_3+\left({\mathrm{I}}_2-{\mathrm{I}}_1\right){\mathsf{\omega }}_2{\mathsf{\omega }}_1^{}={\mathrm{m}}_3^{}\end{array}$$

(2)

where ω_i denotes the rotational velocity vector; I_i, the principal mass moments of inertia. The moment m_i is the sum of moments produced by the contact forces and the applied forces.

The 3DEC model is shown in Figure 5a. The geometry was based on drawings of the structure. Some simplifications were adopted, for example the smaller blocks around the tower openings were not considered. In order to reproduce the correct structural deformability, the joint stiffnesses are based on the dimensions of the blocks adjacent to the joint. The normal joint stiffness is given by k_n = E/d, where E is the Young’s modulus and d the average of the dimensions normal to the joint of the 2 blocks in contact; the shear stiffness is given by k_s = G/d, where G is the shear modulus. These formulas were applied to assign the stiffness to each group of joints (e.g., bed joints and vertical joints of each layer of the towers, joints of the pillars, and arched gates). A Coulomb friction law was adopted as the joint constitutive model, assuming a friction angle of 30°, zero cohesion, and zero tensile strength. This friction angle is a conservative estimate, in the lower range of the available data. The material properties are summarized in

Table 1. Material properties.

Unit weight	24 kN/m3
Young’s modulus	2800 MPa
Shear modulus	860 MPa
Joint friction angle	30⁰
Joint cohesion	0
Joint tensile strength	0

.

In 3DEC, static solutions under gravity loading are first obtained by a relaxation solution algorithm. Then, the dynamic analysis is performed by prescribing the seismic input motion at the base block. The time step required for stability of the explicit algorithm was in the order of 10^-5 s, which allows an accurate representation of the block motion and contact updates. Histories of velocity, displacement, and normal and shear stresses were recorded at the locations shown in Figure 5b, where larger movements are to be expected.

4. Nonlinear Dynamic Analysis

Seismic behavior of the Frontinus Gate was simulated through non-linear dynamic analyses. The objective is to discuss the dynamic global response of the masonry gate to an earthquake ground motion which is consistent with the earthquake hazard levels, in accordance with the Turkish Building Seismic Code 2019. The first seismic input was given by a recorded earthquake, the velocity records of the event that occurred in Dinar (Mw 6.0) on Oct 10, 1995 downloaded from ground PEER NGA strong motion database. In addition to the real ground motion, four synthetic acceleration time series were generated. First, for the coordinates of the Frontinus Gate, target response spectrums were established for two ground motion levels. The first earthquake ground motion level has a 2% probability to be exceeded in 50 years, with a return period of 2475 years (level 1). The second earthquake ground motion level has a 10% probability to be exceeded in 50 years, with a return period of 475 years (level 2). For each target response spectrum, two acceleration time series were generated based on inter plate regimes and linear site effect.

Selected and created earthquakes were best fit with the adopted target spectrums. Comparison of the response spectrum of the Dinar earthquake (Mw 6.0) on Oct 10, 1995, with the target spectrum for ground motion level 2, is shown in Figure 6. Comparisons of the response spectrum of synthetic loadings with corresponding target spectrums are plotted in Figure 7.

These real and synthetic ground motions were applied as base excitation to the numerical model. It is assumed that the structure is subjected to earthquakes in both horizontal directions. Due to the fact that the seismic effects in two orthogonal directions are unlikely to reach their maximum value at the same time, the 30 percent orthogonal loading rule is applied. Since the two time series were generated for each ground motion level, 100% of the synthetic earthquake is combined with 30% of the other synthetic ground motion in the orthogonal direction. For ground motion level 1, two dynamic analyses were performed, as for 100% of the first synthetic earthquake combined with 30% of the second synthetic ground motion, and 30% of the first synthetic earthquake combined with 100% of the second synthetic ground motion. The same procedure was applied for the dynamic analysis of the second ground motion level.

Identifying the dynamic characteristics of the masonry structures under the assumption of elastic behavior is an important first step in the study, as it provides a good understanding of the dynamic response under low intensity ground shaking, until nonlinear phenomena become dominant. The first ten natural frequencies and eigenmodes of vibration were calculated for the rigid block model. The kinematic variables of the system of rigid blocks are the six degrees of freedom of each block as three translations and three rotations. The global stiffness matrix for the rigid block system is assembled, based on the system deformability governed by the joint stiffnesses. The mass matrix is obtained from block masses and rotational inertia terms [36]. Vibration mode shapes as well as their corresponding frequencies were determined. Frequencies for five dynamic eigenmodes are given in

Table 2. Frequencies for 5 dynamic eigenmodes.

Mode	Frequency (Hz)
1	11.60
2	12.00
3	18.32
4	19.98
5	22.64

. The shapes of modes 1, 3, and 5 are plotted in Figure 8, Figure 9 and Figure 10.

The first period is approximately 0.09 seconds, corresponding to a bending mode shape on the horizontal plane. Mode 3 involves shear and torsional movements, and mode 5 mostly the larger tower. It can be seen from the displacement magnitude contours for mode 3 and mode 5 that displacements increase along the height of the round tower.

In the analysis of the response under seismic action, the nonlinear properties of the joints were applied. First, static analysis was carried out under self-weight. Then, nonlinear dynamic analyses were performed under seismic loading. For the Dinar, Turkey 1995 earthquake with a duration of 29 s and increments of 0.005 s, the masonry gate model experienced the seismic excitations in horizontal and vertical directions simultaneously. The structural response of the numerical model subjected to the Dinar earthquake in terms of contour diagram of the final displacement magnitude is presented in Figure 11. As is seen in this figure, the displacements increase along the height of the round tower and the maximum displacement of 7.1 cm is obtained, particularly in areas where the normal stresses are low. The connection of the round tower with the arched gate is also a critical location.

In the following runs, for each ground motion level the numerical model subjected to synthetic earthquakes in two horizontal directions with a duration of 11 s. A total of four dynamic analyses were performed under synthetic loading combination. According to the results obtained from ground motion level 1, under 30% of the first synthetic earthquake combined with 100% of the second synthetic ground motion, the highest round tower collapsed, while the three masonry gate arches were able to sustain such a large action (PGA about 1g). Representative responses of this analysis in terms of the contour diagram of final displacement magnitude are depicted in Figure 12. For the second combination as 100% of the first synthetic earthquake combined with 30% of the second synthetic ground motion, some small blocks fell and the highest round tower was very close to collapse, a result consistent with the previous run. For ground motion level 2, a lower level of seismic input, under the first case as 100% of the first synthetic earthquake combined with 30% of the second synthetic ground motion, the maximum displacement was 0.16 m (Figure 13), whereas under the second combination the tower experienced maximum displacement of 0.13 m. Again, the top of the higher tower showed more displacements, but no collapse was observed for this lower seismic level. In addition to the displacements, the normal and shear stresses observed from dynamic analyses were evaluated. Considering the stresses, a much better seismic response was observed under the Dinar earthquake than for the synthetic seismic excitations. Figure 14 shows some damage patterns as openings and dislocation of stones under the Dinar, Turkey 1995 earthquake and the synthetic earthquakes. The results predicted by the model for the Dinar earthquake are compatible with the available data which indicate that no major damage was detected in the structure after this seismic event. It can be concluded that when the numerical model is subjected to real or synthetic seismic excitations, the main damage was observed in the higher elevations of the round towers, around the door openings, and also extending to the arch intrados. These results show that the expected forms of structural damage to this gate due to seismic excitations are sliding, dislocation and falling of stones mostly along the towers when the horizontal inertial forces exceeded certain limits. Under both combination of synthetic earthquake for ground motion level 1, the highest tensile stresses were observed along the round tower, which would have resulted in a flexural failure mechanism. As shown in representative figures (Figure 14), openings and dislocation of stones took place due to the sliding on horizontal bed joints, which are associated with shear behavior. It is seen that the observed openings have a pattern which would cause typical diagonal cracking.

Under real and synthetic earthquakes, time histories of displacements relative to the base of the model were evaluated at five points, shown in Figure 15. It is obvious from these that under the Dinar earthquake and the synthetic motions, the round towers and the gate show quite different movements, with the arched gate displaying lower values. The maximum values of displacements relative to the model base during dynamic analyses for real and synthetic earthquakes for the five points shown in Figure 15 are given in

Table 3. Maximum displacement relative to base.

Maximum Relative Displacement (m)
Dinar, Turkey 1995 Earthquake	Point 1	0.0416
	Point 2	0.0826
	Point 3	0.1024
	Point 4	0.0226
	Point 5	0.0440
Ground motion level 1 Case 1: 100% of the first synthetic earthquake combined with 30% of the second synthetic earthquake	Point 1	0.3472
	Point 2	0.7614
	Point 3	0.9380
	Point 4	0.2090
	Point 5	0.2853
Ground motion level 1 Case 2: 30% of the first synthetic earthquake combined with 100% of the second synthetic earthquake	Point 1	0.6616
	Point 2	collapsed
	Point 3	collapsed
	Point 4	0.1890
	Point 5	0.6338
Ground motion level 2 Case 1: 100% of the first synthetic earthquake combined with 30% of the second synthetic earthquake	Point 1	0.0621
	Point 2	0.1444
	Point 3	0.1298
	Point 4	0.0440
	Point 5	0.0709
Ground motion level 2 Case 2: 30% of the first synthetic earthquake combined with 100% of the second synthetic earthquake	Point 1	0.0906
	Point 2	0.1841
	Point 3	0.1981
	Point 4	0.0488
	Point 5	0.0438

. We can see that point 4, located at the masonry arch, experienced relatively lower displacements under real and synthetic loadings. Although point 1 is at the highest point, in some runs, the maximum displacement occurred at point 3 (near the gate-tower junction) under the real earthquake and three of four synthetic loadings, showing that this area is subject to intense motion.

5. Discussion

The Frontinus Gate in Hierapolis in Denizli, Turkey has a very long history. This monumental masonry gate is composed of three openings in squared travertine blocks, with masonry arches, flanked by two round towers. It is difficult to investigate the structural health of existing masonry monuments using seismic code requirements that are prepared for new masonry constructions. The work carried out was aimed at getting a better insight on the seismic response of the Frontinus Gate. To accomplish this objective, a discrete element model was created using the present geometry of the structure. Dynamic nonlinear analyses were performed under real earthquakes and synthetic ground motions. Interpretation of the numerical model results lead to the following conclusions:

• Modal analysis, assuming elastic behavior, is a very useful first step to obtain insight into the dynamic characteristics of the structure.
• The nonlinear dynamic analyses indicate that the Frontinus Gate, in the current state, under an earthquake with a magnitude of 6, would be expected to suffer moderate damage. Under the synthetic earthquake with a return period of 475 years, the observed structural response and damage patterns are serious although it does not mean the collapse of the monument.
• In the case of occurrence of an earthquake with a return period of 2475 years, the models show that the monument would not be able to withstand it, resulting in a large number of blocks collapsing completely.
• Considering the seismic sensitivity of the Frontinus Gate, the results obtained are consistent with the present state of the structure, with the round towers having more seismic damage potential than the masonry arches.
• Although the Frontinus Gate is still standing at the site, showing an outstanding endurance, under seismic excitation with a magnitude of 6 or greater, a high vulnerability of large portions of the structure can be expected.