Seismic Vulnerability of Roman Aqueduct Bridge Structures: Numerical Modelling and Application to the Antioch-on-the-Orontes Aqueduct (Antakya, Türkiye)

1. Introduction

Roman aqueducts are among the most emblematic achievements of ancient hydraulic engineering. Built throughout the Mediterranean world, they supplied water for domestic, agricultural, industrial, and recreational purposes and played a central role in the functioning of Roman cities. Many of these structures are still partly preserved, testifying both to the sophistication of Roman engineering and to the remarkable long-term durability of their construction techniques. Beyond their archaeological and historical value, these monuments also provide an opportunity to investigate how masonry structures behaved under repeated environmental forcing over nearly two millennia.

One of the major hazards that may have affected these structures is seismic shaking. The Mediterranean region is characterized by high seismic activity, with frequent moderate to large earthquakes and several tectonically active belts, particularly in Italy, Greece, and Türkiye. This regional overlap between the distribution of Roman aqueducts and Mediterranean seismicity is illustrated in Figure 1. While some aqueduct damages located near active faults have long been suspected to be earthquake-related (e.g., [1,2]), such an interpretation is not always straightforward. Structural degradation may also result from weathering, poor maintenance, foundation problems, human disturbance, or other natural processes. Except in rare cases where an aqueduct is directly displaced by fault rupture [3,4], the seismic origin of observed damage remains a working hypothesis that must be tested against alternative explanations.

This question has acquired renewed contemporary relevance. On 6 February 2023, the Kahramanmaraş earthquake sequence included a Mw 7.8 mainshock followed about nine hours later by a second major event of Mw 7.5–7.6. Together, these earthquakes produced catastrophic destruction across southern Türkiye and northern Syria and re-emphasized the seismic threat affecting the broader Antioch/Antakya region [7,8,9,10]. Although the numerical modelling presented here predates that earthquake sequence and does not incorporate post-2023 observations from the aqueduct itself, the 2023 disaster underscores the continuing relevance of assessing the seismic vulnerability of historical masonry structures in this part of the eastern Mediterranean.

Numerical modelling provides a useful way to evaluate whether the geometry and material properties of a structure are compatible with earthquake-induced damage. Previous studies have demonstrated the value of such approaches for individual monuments, including masonry arches and aqueducts, through both archaeoseismological and numerical investigations (e.g., [11,12,13,14,15]). Comparable integrated assessment frameworks have also been applied to other historical masonry monuments in Türkiye, combining damage survey, regional seismicity, geotechnical investigations, ambient vibration testing, and finite-element analysis prior to strengthening interventions (e.g., [16]). However, most existing models are tailored to a single structure and therefore remain difficult to generalize. Roman aqueducts exhibit substantial architectural variability, including differences in arch span, pier height, deck thickness, reinforcement systems, and building materials. A broader comparative framework is therefore needed to identify the structural configurations that are most vulnerable to seismic shaking and to distinguish generic earthquake-related damage patterns from site-specific effects.

The present study addresses this issue by combining a generic finite-element assessment of arch-and-pier aqueduct bridge structures with a case study from the Antioch-on-the-Orontes aqueduct at Harbiye, near present-day Antakya in southern Türkiye. The aim is twofold: first, to investigate how basic architectural parameters and material properties influence the seismic response of Roman aqueduct bridges; and second, to test whether the damage distribution observed in the Antioch aqueduct is consistent with earthquake-induced stress concentrations, and whether the identified repairs may have improved the structural stability of the monument. More broadly, this work aims to support the archaeoseismological interpretation of damaged aqueduct bridges in seismically active regions.

2. Geological and Archaeological Setting

2.1. Antioch-on-the-Orontes/Antakya Within the Regional Tectonic Framework

Ancient Antioch-on-the-Orontes, corresponding to present-day Antakya in Hatay Province, is located within one of the most tectonically complex sectors of the eastern Mediterranean. The city lies close to the junction between the northern Dead Sea Fault system, the southwestern termination of the East Anatolian Fault Zone, and the broader plate-boundary domain that also includes the Cyprus Arc [7,10,17,18]. This setting makes the Antakya–Amik region a structurally complex fault-interaction zone rather than an area controlled by a single fault trace (Figure 2a).

At the local scale, the Antioch aqueducts are situated near the northern Dead Sea Fault system. In this region, the fault system includes the Hacıpaşa fault strand to the south and the Karasu fault strand extending northward along the Karasu Valley and the southwestern East Anatolian Fault system (Figure 2a; Akyüz et al., 2006; Karabacak et al., 2010; Karabacak and Altunel, 2013). The Hacıpaşa Segment is particularly relevant for the Harbiye aqueduct because it lies less than about 25 km east of the monument. Its proximity makes it one of the most important local seismic sources to consider when discussing possible earthquake damage to the aqueduct, although the broader Antakya region may also be affected by ruptures on neighbouring fault strands.

The long-term seismic importance of the Antioch/Antakya region is well established. Historical catalogues mention at least sixteen earthquakes that caused severe damage in Antioch-on-the-Orontes during the last two millennia [20,23,24]. During the Roman period, particularly destructive earthquakes occurred in 64 BC, 37 AD, and 115 AD, each of which caused heavy damage to the city and its infrastructure. More recent large events include the 1822 and 1872 earthquakes [20,21,22,25].

The 6 February 2023 Kahramanmaraş earthquake sequence also affected this region. The first mainshock (Mw 7.8) initiated on or near the Narlı fault, a secondary fault south of the main EAF strand, and then propagated bilaterally along the East Anatolian system over about 310 km, involving several fault segments, including the Karasu/Amanos sector (Figure 2a); roughly nine hours later, a second major earthquake (Mw 7.5–7.6) ruptured the Çardak–Sürgü fault system [7,8,9,10]. The sequence caused catastrophic damage across southern Türkiye and northern Syria and was followed by intense aftershock activity, including the Mw 6.4 Hatay earthquake of 20 February 2023 near the Samandağ–Antakya fault zone [9,26]. Although the numerical modelling presented in this paper predates the 2023 sequence and does not incorporate direct post-earthquake observations from the aqueduct itself, these events clearly underline the continuing seismic exposure of the Antioch/Antakya region.

2.2. The Antioch Water-Supply System

The structure analyzed in this study formed part of the broader water-supply system of Antioch-on-the-Orontes, which conveyed water from the springs of Daphne, near present-day Harbiye, to the ancient city. The general layout of this hydraulic system and the mapped aqueduct remains are shown in Figure 2c. This hydraulic network comprised bridges, channels, tunnels, and conduits converging toward Antioch and, during the Roman period, was complemented by wells and cisterns [19]. According to these authors, the first water-supply works were initiated shortly after the foundation of the city, in the 3rd century BC. Later historical sources, especially the Byzantine chronicler Malalas (490-578 AD), indicate that additional hydraulic constructions were undertaken under Julius Caesar around 50 BC. The same sources further suggest that parts of the aqueduct system were repaired, and probably enhanced, after the earthquake of 37 AD under Caligula, and again after the earthquake of 115 AD under Trajan and Hadrian [19,27].

The Daphne branch of the system drew water from a source located on a calcareous plateau about 7 km from the city. In order to maintain gravity-fed flow from this upland source to Antioch, the hydraulic line had to cross a rugged topography deeply incised by valleys, requiring the construction of major bridge-aqueduct structures. In the Harbiye sector, two arched aqueduct bridges crossed the Orontes valley. The present study focuses on the larger of these structures, which constituted one of the major engineering elements of the Antioch water-supply network.

2.3. The Harbiye Bridge-Aqueduct

Among the two arched aqueduct bridges that crossed the Orontes valley in the Harbiye sector, the structure investigated here is the larger double-level bridge-aqueduct. It is located at 36.151826°N, 36.154263°E and is attributed to the Roman period [27]. The location and present-day remains of the monument are shown in Figure 2b. The modelled monument is approximately 61.5 m long and reaches up to 30 m in height. It consists of two superposed structural units: a lower structure built with large piers and arches, and an upper structure formed by a more regular series of smaller arches. These two structural units and the main preserved architectural features are illustrated in Figure 3.

The geometry of the lower structure is strongly controlled by local topography and foundation conditions. Its piers rest on ophiolitic bedrock, which is highly deformed and sheared, especially on the southern side of the valley [28]. Because the valley is asymmetrical, with steep north-facing slopes and stepped south-facing slopes marked by small terraces, the bridge had to be anchored on four piers of different heights. From south to north, Pier 1 is the largest, reaching about 30 m in height and 15 m in width; it is founded partly on the riverbed and partly on the ophiolitic substratum. Pier 2 lies 21 m farther north and is founded about 6 m higher, on ophiolitic bedrock and sediments. Farther north, the structure continued on three narrower and more closely spaced piers, one of which was later destroyed during road construction. Piers 1 and 2 were originally connected by a very large arch spanning the river. This major arch has collapsed, and remnants of it are reportedly still visible in the Orontes riverbed. Because its exact original geometry is unknown, the missing arch is reconstructed in the model using a simplified inferred shape.

The upper part of the monument appears to have consisted of a regular series of arches separated by piers of similar height. However, because the central part of the aqueduct has collapsed, this upper structure is preserved only on both sides of the river. The reconstruction used for modelling therefore assumes a continuous and regular upper-level arrangement based on the surviving remains. The standing remains also show localized travertine cover on some preserved piers and arches, as indicated in the original field documentation and figure descriptions.

The monument preserves a marked material and structural dichotomy between its lower and upper parts. The lower structure was built of large cut stones quarried from nearby calcareous outcrops. This stone-built part was subsequently modified by repairs, including the infilling of arches and the addition of buttresses, especially around the main lower piers. Such reinforcement of aqueduct piers by means of buttresses and walls is not exceptional in Roman aqueduct architecture and has been documented elsewhere, for example in the Hadrianic aqueduct of Corinth [29]. Most of these reinforcements were constructed using cemented rubble stones. Several visible damage features affect the lower structure, including diagonal cracks, shearing of stone blocks, and displaced or rotated blocks near the base of Pier 1. These features may be interpreted as likely earthquake-induced.

By contrast, the upper structure is built of brick facings with a core made of cemented rubble stones, a method called opus caementicium developed after 2nd century BC [30]. The contrast in building technique between the two levels strongly suggests that the upper part belongs to a later construction or reconstruction phase. This interpretation is supported by the material analyses of Benjelloun et al. [27], who sampled and characterized bricks and mortars from the monument and showed that the clay used for the bricks was sourced from the Amik Lake area, located a few kilometres from Antioch. Radiocarbon and paleomagnetic dating carried out on charcoals preserved in the bricks and mortars of the upper structure indicate that this phase can be attributed to the reign of Caligula (37–41 AD), consistent with the historical account of the Byzantine chronicler Malalas (490-578 AD) that the aqueduct was rebuilt or substantially repaired after the 37 AD earthquake.

Taken together, these characteristics make the Harbiye bridge-aqueduct an exceptional case study. It is not only part of a major urban hydraulic system, but also a monument in which structural hierarchy, material contrasts, visible damage, inferred collapses, and later reinforcements can all be examined together. This combination of archaeological observations, historical information, and preserved architectural complexity makes it particularly suitable for testing whether the present damage pattern is mechanically consistent with seismic loading and whether the observed repairs correspond to structurally effective responses to past earthquakes.

3. Materials and Methods

3.1. Numerical Modelling Framework

3.1.1. Modelling Strategy and Rationale

Masonry structures can be modelled using different numerical approaches depending on the objective of the study, the available data, and the architectural complexity of the monument. Two broad classes of approaches are commonly used for arch structures: finite element methods (FEM) and discrete element methods (DEM) [11,12,13,14,31]. DEM approaches are well suited to representing the composite nature of masonry, including the interactions between blocks, mortar, and joints. However, for complex architectural systems they remain difficult to implement at the scale of an entire monument. In FEM approaches, masonry is generally treated as a homogenized material with equivalent mechanical properties. Although this simplification does not explicitly resolve the internal heterogeneity of stones, bricks, mortar, and joints, it provides an efficient framework for comparing the structural response of multiple architectural configurations and material combinations.

In the present study, our objective is not to reproduce every local construction detail of a specific monument, but to identify the first-order controls on the seismic response of aqueduct bridge structures and to locate the zones most susceptible to stress concentration under seismic loading. This approach is appropriate to archaeological contexts, where uncertainties commonly affect the original geometry, the internal structure of the masonry, the exact properties of the materials, and the chronology of repairs and rebuilding phases. We therefore adopted a simplified finite-element strategy in which aqueduct bridges are treated as assemblages of two main substructures, piers and arches, and masonry is represented as a homogenized material. This strategy allows systematic comparison between generic aqueduct geometries and the Antioch-on-the-Orontes case study, while remaining compatible with the level of archaeological and architectural information available. The overall numerical workflow, including the static and dynamic analyses, is summarized in Figure 4.

3.1.2. Numerical Tool

The numerical analyses were performed using FINELG, a finite-element structural analysis software developed at the University of Liège. FINELG allows both static and dynamic analyses under linear or nonlinear formulations. In this study, we used a linear formulation as a first-order approximation of the seismic response of aqueduct structures. The results are therefore interpreted in terms of relative displacement amplitudes and stress distributions, rather than as simulations of nonlinear failure, crack initiation, or progressive collapse.

3.2. Model Construction

3.2.1. Geometrical Representation and Meshing

Because the aim of the study is to investigate the in-plane response of aqueduct bridge structures rather than their full out-of-plane behaviour, we adopted a 2.5D modelling strategy instead of a full 3D volumetric representation. The geometry was discretized using quadratic and triangular thick-shell elements with four and three nodes, respectively. Each node is characterized by six degrees of freedom, corresponding to three translations and three rotations. In practice, the chosen direction of the imposed earthquake motion primarily activates the in-plane degrees of freedom, whereas the out-of-plane behaviour is simplified.

The base of each model was assumed to be fixed in order to represent anchoring to the substratum and to avoid detachment. In this first-order approach, the structures were assumed to be directly founded on bedrock. Site effects and local soil amplification were therefore not considered as explicit parameters of the analysis. This simplification was adopted to keep the focus on the structural response of the aqueduct itself rather than on the coupling between structure and site conditions.

3.2.2. Mechanical Parameters

For the generic models, masonry was treated as homogeneous and isotropic, so that each model was defined by a single set of mechanical parameters: Young’s modulus (E), Poisson’s ratio (ν), and volumetric mass density (ρ). These parameters control the stiffness, elastic response, and inertial behaviour of the structure under seismic loading. Although different mechanical properties could theoretically be assigned to individual structural components, the generic models were designed to test the first-order influence of material properties on aqueduct stability rather than to reproduce the full heterogeneity of specific masonry assemblages.

The selected parameter ranges were defined from published values for the principal material categories used in Roman aqueduct construction. Roman aqueducts were commonly built using locally available raw materials, especially stone, brick, and cemented rubble-stone materials. Stone masonry was generally made from rocks available at or near the construction site, most commonly limestone or sandstone. Intact rocks may display Young’s modulus values ranging from about 1000 to 70,000 MPa and Poisson’s ratios from about 0.05 to 0.33 [32]. These values, however, cannot be directly transferred to masonry assemblages, because cut blocks are separated by joints and mortar, which act as preferential weak zones and reduce the effective stiffness of the structure. As a result, the effective Young’s modulus of stone masonry is expected to be lower than that of the intact rock from which the blocks were quarried.

Ancient bricks are also mechanically variable. Their properties depend on raw material, production process, firing temperature, mineralogical transformations during firing, and subsequent weathering. Vitruvius already emphasized the importance of raw-material selection for brick quality, although local availability was probably a major practical constraint [33]. Modern studies show that firing conditions and mineralogical changes can strongly influence brick strength and elastic properties [34,35]. In addition, the current properties of ancient bricks may differ from their original properties because of long-term weathering and alteration [36]. Published values remain relatively scarce because sampling from historical monuments is often limited, but reported Young’s modulus values for ancient bricks commonly fall within the lower part of the range considered here, with values of about 1000–7000 MPa reported by Stefanidou et al. [37]. As with stone masonry, the behaviour of a brick assemblage is also influenced by mortar and joints and is therefore generally less resistant than that of individual bricks.

Cemented rubble-stone materials, including opus caementicium, are even more difficult to parameterize because of their heterogeneous and composite structure. Their mechanical behaviour depends on the properties of the mortar, the aggregate, the degree of cementation, and the state of preservation. Roman opus caementicium consists of mortar mixed with stone aggregates, whose nature was commonly determined by the local geological resources available near the construction site . Mechanical tests on ancient opus caementicium indicate Young’s modulus values ranging from about 800 to 9170 MPa and Poisson’s ratios around 0.18–0.19 [38]. Mortars also show large variability. Depending on composition and state of preservation, reported Young’s modulus values range from several hundred MPa to several thousand MPa. Reproduced Roman mortars studied by Giavarini et al. [38] and Degryse et al. [39] yielded values of about 2500–3600 MPa, whereas Özkaya and Böke [40] reported a lower value of about 630 MPa for a weathered lime-and-aggregate Roman mortar from the Serapis temple at Pergamon.

Taken together, these data justify the broad sensitivity ranges adopted in this study and reinforce the need to treat the assigned properties as effective values for masonry assemblages rather than as intrinsic properties of individual stones, bricks, or mortar components. For the reference generic models, we used E = 7500 MPa, ν = 0.25, and ρ = 2700 kg m⁻³. Sensitivity tests were then carried out by varying Young’s modulus between 1000 and 70,000 MPa, Poisson’s ratio between 0.15 and 0.40, and density between 1600 and 2700 kg m⁻³. These tests were designed to assess how the choice of material parameters influences resonance frequencies, displacement amplitudes, and stress distribution.

For the Antioch model, different properties were assigned to the lower and upper parts of the structure in order to reflect the observed material contrast between the lower stone-built unit and the upper brick-and-rubble superstructure. The lower limestone structure was assigned E = 15,000 MPa, ν = 0.25, and ρ = 2300 kg m⁻³, whereas the upper structure was assigned E = 9000 MPa, ν = 0.25, and ρ = 1800 kg m⁻³.

3.3. Structural Configurations

3.3.1. Generic Aqueduct Models

To explore the influence of architectural design on seismic response, we constructed a series of simplified generic aqueduct bridge models. The corresponding geometries are shown in Figure 5 and Figure 6. These models represent two broad families of aqueduct geometry related to different geomorphological settings. The first family comprises narrow-arched structures with relatively slender piers (Models I). These configurations are representative of aqueducts built across plains, plateaus, or gently undulating relief. Model 1.0 was used as the reference geometry. Model 1.1 differs from the reference model by higher piers, whereas Model 1.2 retains the same pier geometry as Model 1.0 but includes a thicker deck.

The second family comprises wide-arched structures (Models II), designed to represent aqueducts crossing more rugged topography and deeply incised valleys. Model 2.0 is characterized by a single large arch between two comparable piers. Model 2.1 introduces a buttress, whereas Model 2.2 replaces the buttress with a larger pier. In addition to these major geometrical contrasts, the influence of minor arch-shape variations was explored by comparing different arch geometries while keeping the rest of the structure unchanged. The tested arch shapes are presented in Figure 7. These generic models were designed to evaluate how changes in pier height, arch width, deck thickness, and reinforcement by buttresses or massive piers affect the dynamic response of the structure.

3.3.2. Model of the Antioch-on-the-Orontes Aqueduct

A third family of models (Models III) was developed for the Harbiye bridge-aqueduct. The three reconstructed construction stages used for the Antioch model are shown in Figure 8. These models combine the wide lower arch structure and the smaller superposed arches observed in the monument. Their geometry was based on field measurements, observations of the standing remains, and a limited number of necessary architectural inferences. Because some parts of the monument are missing, the reconstructed geometries necessarily remain hypothetical in part. The missing sections were completed assuming a symmetrical and regular organization consistent with the preserved remains.

Model 3.0 represents the inferred initial lower stone structure before the construction of the upper level. Model 3.1 represents the structure before the identified repairs, combining the lower stone-built unit and the upper brick-and-rubble unit but without arch infill or buttresses. Model 3.2 represents the final Roman construction stage, including the visible reinforced zones such as infilled arches and buttresses. This sequence of models was designed to test how successive modifications affected the seismic behaviour of the monument and whether the repaired configuration corresponds to an increase in structural stability.

3.4. Seismic Excitation and Damping

3.4.1. Structure Excitation

The structural response was evaluated using both static and dynamic analyses. In the static approach, an equivalent static load was calculated and applied to the structure. This preliminary analysis is computationally efficient and was used primarily as a first check of the mesh quality and of the consistency of the assigned mechanical properties. In the present study, however, the main results are derived from dynamic analyses, which more explicitly simulate the response of the aqueduct structures to seismic shaking.

In the dynamic analysis, the seismic input is applied in the form of an acceleration time series imposed at the base of the model, at the anchoring points. The excitation can in principle be defined using either a real accelerogram or a synthetic signal; here, we used real strong-motion records. One horizontal component of each accelerogram was applied as basal excitation in the in-plane direction of the structure.

An additional parameter required for the dynamic analysis is the damping, which represents attenuation within the structure. Damping was modelled using the linear Rayleigh approximation, in which the damping matrix is expressed as a combination of the stiffness and mass matrices: C=αK+βM where C is the damping matrix of the physical system, M is the mass matrix, and K is the stiffness matrix. In the model, the mass depends on the geometry, virtual thickness, and volumetric mass density, whereas the stiffness depends on both the geometry and the mechanical properties assigned to the building materials.

The Rayleigh coefficients α and β were calibrated using the two characteristic eigenfrequencies of the structure and a target damping ratio, following the standard two-frequency Rayleigh damping formulation. In this formulation, the damping ratio is imposed for the two principal vibration modes retained in the analysis, and the resulting α and β coefficients are then used for the dynamic simulations.

3.4.2. Earthquake Input

To simulate strong regional seismic shaking, we used recorded accelerograms from the European Strong-Motion Database [41]. The records were selected to represent strong potentially damaging shaking, using criteria based on relatively large magnitude, short epicentral distance, high peak ground acceleration, and recording conditions representative of bedrock in order to limit site-specific amplification effects. Two ground-motion records with different focal mechanisms were selected in order to test whether the structural response of the aqueduct models depended strongly on source type. The first record corresponds to the 1995 Kozani earthquake, Greece (Mw 6.5), recorded at an epicentral distance of 17 km, with a peak horizontal acceleration of 2.039 m s−2 and a dominant frequency content around 3–3.5 Hz. The second corresponds to the selected E–W component of the 1999 Düzce earthquake, Türkiye (Mw 7.2), recorded at 23 km from the epicentre, with a peak horizontal acceleration of 4.865 m s−2 and a dominant frequency content around 4.5 Hz. The corresponding acceleration time series are shown in Figure 9. In both cases, the horizontal acceleration time series was applied at the base of the model, and the first 40 s of the record were used in the simulations. Because the Düzce record corresponds to a stronger and relatively proximal strike-slip event in a tectonic setting broadly comparable to that of southern Türkiye, it was used as the main excitation for comparing the stress distribution and displacement response of the different models.

Although the 6 February 2023 Kahramanmaraş earthquake sequence would now constitute an obvious candidate for site-specific seismic input, it was not used as the main excitation record in this study. This sequence was an exceptional doublet involving two Mw > 7.5 events on distinct fault branches, with complex rupture propagation and strong source-specific effects. Because our objective was to compare the relative seismic response of different aqueduct configurations under controlled loading conditions, we retained pre-2023 strong-motion records representative of large regional earthquakes. The 2023 sequence is therefore considered here primarily as a major reminder of the present-day seismic hazard affecting the Antioch/Antakya region, whereas its explicit use as input motion would require a dedicated site-specific follow-up study.

3.4.3. Modal Analysis and Damping Ratio

A modal analysis was performed for each model in order to determine the principal vibration modes and their associated resonance frequencies. Although the software can theoretically compute as many modal shapes as there are degrees of freedom, we retained only the two modes with the highest modal mass in the plane of the structure, as these are the most representative of its seismic behaviour. The first mode mainly corresponds to horizontal translation, whereas the second involves deformation in the vertical plane of the structure. These two modal frequencies were used to calibrate the Rayleigh damping coefficients. They also provide a first indication of the susceptibility of the structure to resonance when compared with the dominant frequency content of the input motion. This comparison is relevant because earthquake frequency content depends on magnitude, source-to-site distance, rupture directivity, and site conditions [42]. For bedrock conditions, the frequency range of strongest seismic excitation is commonly considered to lie approximately between 2.5 and 6.67 Hz in Eurocode 8 [43].

The influence of the damping ratio was tested on the reference Model 1.0. The effect of damping on maximum displacement and maximum Von Mises stress is shown in Figure 10. Values between 1% and 5% were explored, showing that when attenuation falls below 5%, both displacement and stress increase markedly. A damping ratio of 5%, consistent with standard seismic engineering practice and with Eurocode 8 [43], was therefore retained for all simulations.

3.5. Output Variables and Damage Criterion

FINELG computes the displacement of each node and the resulting internal stresses at a chosen time step of 0.01 s. In this study, we focused primarily on the Von Mises stress, which provides a scalar measure integrating both normal and shear stress components. Although principal tensile stresses could also have been examined, the Von Mises criterion was retained here as a practical first-order indicator of the overall stress concentration within the structure.

To assess the likelihood of damage, we estimated the global compressive and tensile strength of the masonry using empirical relationships derived from Eurocode 6 [44]. The characteristic compressive strength was estimated using the empirical relation f_k = K f_b^α f_m^β, where f_b is the masonry-unit strength, f_m is the mortar strength, and K, α, and β are parameters defined according to the masonry type. For natural-stone masonry, we used K = 0.45, α = 0.7, and β = 0.3. This yielded an average global compressive strength of about 15 MPa. Because earthquake loading generates both compressive and tensile stresses, and because masonry is particularly vulnerable in tension, a tensile rupture threshold equal to 10% of the compressive strength was adopted. A value of 1.5 MPa was therefore used as a first-order threshold for damage initiation in the models. The purpose of this threshold is not to simulate crack propagation explicitly, but to identify zones where the modelled stress state is consistent with likely structural damage under seismic loading.

4. Results

4.1. Sensitivity of the Models to the Selected Earthquake Records

The two selected earthquake records produce broadly similar response patterns in the models. In particular, they yield the same resonance frequencies and identify comparable zones of elevated stress. The main difference between them lies in the amplitude of the structural response: the stronger Düzce-Izmit record generates larger displacements and higher Von Mises stresses than the Kozani record. For this reason, the results presented below are based on the Düzce-Izmit accelerogram, which provides the more discriminating dynamic input for comparing model behaviour.

4.2. Influence of Structural Geometry in the Generic Aqueduct Models

The modal analyses show that small variations in arch shape have little influence on the first-order dynamic behaviour of the reference narrow-arch structure (Table 1). For the three arch geometries tested in Model 1.0, mode-1 resonance frequencies range from 12.4 to 13.6 Hz and mode-2 frequencies from 13.8 to 14.0 Hz. By contrast, larger geometrical changes significantly modify the modal response. In Models I, the reference structure (Model 1.0) has resonance frequencies of 13.0 and 13.9 Hz, whereas increasing pier height in Model 1.1 lowers them to 8.0 and 8.9 Hz. Increasing deck thickness in Model 1.2 also lowers them, although less strongly, to 10.6 and 11.2 Hz. In Models II, the wide-arch structure (Model 2.0) yields frequencies of 9.5 and 10.4 Hz. The addition of a buttress in Model 2.1 increases these values to 11.3 and 12.5 Hz, while the presence of a larger pier in Model 2.2 gives 11.1 and 15.0 Hz. The modal frequencies of Models I and II are summarized in Table 2. Overall, increasing structural height lowers the main resonance frequencies, whereas reinforcement by buttresses or larger piers raises them.

These modal results are mirrored by the dynamic analyses. The corresponding maximum displacements are summarized in Figure 11. In both model families, the largest displacements are observed in the configurations with the lowest resonance frequencies, while the smallest displacements occur in the stiffest structures. In Models II, Model 2.2 shows the lowest displacements, and its response in this respect is comparable to that of the narrow-arch reference structure. Increasing pier height also shifts the resonance frequencies toward the usual frequency band of earthquake shaking on bedrock, which helps explain the greater sensitivity of the taller configurations.

The Von Mises stress distributions reveal a robust and recurrent spatial pattern in all generic models. These stress distributions are shown in Figure 12. Stress concentrations are systematically localized in two main areas: at the arch springings, that is, the arch-pier junctions, and at the bases of the piers. Among all generic configurations, Model 1.0 shows the lowest overall stress level. Increasing pier height in Model 1.1 produces a marked increase in stresses over a broad zone near the arch-pier junction and also at the pier bases. Increasing deck thickness in Model 1.2 also slightly enhances stresses at the arch-pier junction, whereas stresses remain relatively low in the deck itself, which appears comparatively stabilized by its own weight. The wide-arch Model 2.0 exhibits a broader stressed zone along the arch and higher stress concentrations in both piers than the narrow-arch reference model. By contrast, the reinforced Models 2.1 and 2.2 show a strong reduction in the extent and amplitude of the stressed zones, with the large-pier solution producing the strongest reduction.

4.3. Influence of Material Properties

The resonance frequencies of the reference model are primarily controlled by Young’s modulus and volumetric mass density. The sensitivity of the resonance frequencies to these parameters is shown in Figure 13. Increasing Young’s modulus raises the resonance frequencies, whereas increasing density lowers them. The influence of Poisson’s ratio is comparatively weak. Across the tested range, mode 1 decreases by only 0.23 Hz and mode 2 increases by only 0.34 Hz, showing that Poisson’s ratio has only a limited effect on the modal behaviour of the structure.

The dynamic analyses show that material properties also exert a strong control on structural response. The corresponding effects on maximum displacement and Von Mises stress are shown in Figure 14. Maximum displacement decreases as Young’s modulus increases, whereas it increases with volumetric mass density and, more slightly, with Poisson’s ratio. The Von Mises stresses follow the same first-order tendencies: they increase with Poisson’s ratio, increase with density, and rise sharply as Young’s modulus decreases. These results indicate that stiffness and density are the dominant material controls on seismic response. For comparable stiffness values, lighter materials such as brick tend to reduce both maximum displacements and Von Mises stresses relative to heavier stone constructions.

4.4. Application to the Antioch-on-the-Orontes Aqueduct

The three Antioch models display lower resonance frequencies than the generic reference structures, consistent with their larger scale and more complex architecture; the corresponding values are summarized in Table 3. Model 3.0 has the highest resonance frequencies, with values of 3.3 Hz for mode 1 and 9.4 Hz for mode 2. In Models 3.1 and 3.2, the increased overall height of the monument lowers the modal frequencies. Model 3.1 yields values of 1.9 and 6.5 Hz, whereas Model 3.2 gives 2.3 and 7.7 Hz. The addition of the secondary arches in the lower structure of Model 3.1 contributes to this frequency decrease. More generally, the Antioch models place at least some of the principal modes within the frequency range of earthquake shaking, especially for mode 1 of Model 3.0 and mode 2 of Model 3.1.

Despite uncertainties in the exact succession of building phases, the Antioch models consistently identify several mechanically sensitive zones. The computed Von Mises stress distributions for the three Antioch construction stages are shown in Figure 15. In all three models, high Von Mises stresses occur at the base of Pier 1, which corresponds to one of the principal damage zones observed in the standing remains. More generally, the lower structure of the bridge contains the main stress concentration zones, especially near the large river-crossing arch and at the bases of the main piers.

The three structural stages differ markedly in their computed response (Figure 15). Model 3.0, representing the inferred initial limestone structure, remains comparatively stable, although cracking is still predicted at the base of the southern Pier 1 and the northern Pier 2. Model 3.1 is the most vulnerable configuration. In this case, Von Mises stresses reach values up to 9 MPa, and the three arches of the lower structure are the most stressed parts of the monument, probably because of the architecture assigned to Pier 2, represented in the model by a thick deck. The destruction of the small northern arch is still visible in the field, and this arch was later filled with opus caementicium and reinforced by a buttress. Model 3.2, which includes the reinforced state of the monument, shows a clear reduction in stress concentrations relative to Model 3.1 and a more stable overall response, particularly in connection with the addition of the buttress and the infilled repair zones.

A final result concerns the spatial correspondence between the computed stress field and the present state of preservation of the monument. The sectors identified as damaged or repaired in the standing remains coincide with the areas of highest computed stress in the unrepaired model. In addition, the piers subjected to the lowest Von Mises stresses are the ones that remain standing today, whereas the strongest stress concentrations occur in the central part of the wide arch, which is the part of the monument that has collapsed.

5. Discussion

5.1. First-Order Seismic Vulnerability of Arch-and-Pier Aqueduct Structures

The modelling results indicate that Roman aqueduct bridge structures are not equally sensitive to seismic loading across their geometry. Instead, the structural response consistently concentrates in two recurrent zones: the arch-pier junctions and the basal part of the piers. These two areas therefore emerge as the most likely locations where damage would preferentially develop under strong ground motion. In this respect, the generic models provide a useful first-order framework for recognizing structurally vulnerable sectors in aqueduct bridges built on arch-and-pier schemes.

Among the geometrical parameters tested, pier height appears to be a primary control on vulnerability. The models with higher piers combine lower resonance frequencies, larger displacements, and wider stress concentration zones than the other configurations. In practical terms, this means that tall and slender structures are more sensitive to dynamic loading than lower and more compact ones. Wide arches also prove mechanically less favourable than narrow arches when they are supported by ordinary piers, because they generate broad stressed zones both at the pier bases and at the arch-pier connections. In the unreinforced wide-arch configuration, the stressed zone may also extend obliquely from the pier base toward the arch, suggesting a possible pathway for large deformation or pier failure under strong shaking.

The simulations also show that structural weight is an important secondary control. A thick deck increases the load carried by both the arches and the piers and tends to amplify stresses in the already identified weak zones. This result helps explain why many Roman aqueducts avoided unnecessarily massive decks and instead relied on repeated superposed arches, which reduce structural mass while preserving the hydraulic gradient. In high structures, superposed arches may also have contributed to stabilizing tall piers by reducing their effective free height and limiting lateral oscillation. In that sense, the model results are consistent with an architectural logic in which material economy and structural stability converge.

The wide-arch models further show that local reinforcement can substantially improve structural behaviour. Buttresses markedly reduce both the extent of the damage-prone zones and the amplitude of critical stresses, particularly around the reinforced pier. Even more effective is the presence of a large pier, which provides a stronger reduction in stress and displacement than the buttressed solution. At the scale of the simplified models, this suggests that widening or strengthening support elements is one of the most effective ways of improving the seismic stability of large-span aqueduct bridges.

5.2. Implications of Material Choice

The material tests indicate that seismic response is favoured by combinations of high stiffness, low Poisson’s ratio, and especially low volumetric mass density. In comparative terms, density plays a particularly important role because it directly increases inertial loading during shaking. This helps explain why lighter construction materials may provide a structural advantage even when their elastic stiffness does not differ dramatically from that of heavier stone assemblages.

Roman aqueducts were built using a limited range of major material categories, principally stone masonry, brick masonry, and cemented rubble-stone materials such as opus caementicium. Stone was commonly selected according to local geological availability, so that the nature of the masonry depended strongly on the construction site and its surrounding resources [45]. In mechanical terms, however, cut-stone masonry generally implies higher volumetric mass density than brick-based construction, while its effective stiffness is reduced by the presence of joints and mortar. The models therefore suggest that, for comparable effective stiffness values, replacing heavier stone assemblages with lighter brick-based materials can reduce both displacement amplitudes and stress concentrations.

Brick construction may also have offered practical architectural advantages. Bricks provide more regular units than quarried stone blocks, can be produced in standardized shapes and sizes, and are particularly suitable for the construction of arches and repeated openings. These constructive advantages probably contributed to the increasing use of brick in Roman architecture, especially from the early Imperial period onward, with large-scale brick production reaching a high level during the Hadrianic period [45,46]. In the context of aqueduct bridges, the regularity and lower density of brick masonry would have been favourable not only for construction efficiency, but also for reducing inertial loads during seismic shaking.

The same reasoning applies to cemented rubble-stone materials such as opus caementicium. Because they combine mortar and stone aggregates into a relatively light and adaptable mass, these materials appear particularly suitable for large support elements, arch infills, and reinforcement zones. Experimental work on imperial Roman mortars further suggests that Roman cementitious materials could combine relatively low stiffness with appreciable fracture resistance and crack-arresting capacity, properties that may have contributed to the long-term mechanical resilience of large Roman masonry and concrete structures [47,48]. These results do not imply that Roman builders selected materials primarily for seismic reasons, but they do suggest that some widely adopted Roman construction choices were mechanically favourable in earthquake-prone regions.

5.3. Archaeoseismological Significance of the Antioch Aqueduct Case

The Antioch-on-the-Orontes aqueduct provides a particularly informative case because it combines visible damage, multiple construction phases, and clear evidence of repairs. In the unrepaired configuration represented by Model 3.1, the highest Von Mises stresses occur systematically in the sectors that are today identified as damaged or repaired, especially in the lower arches and in the reinforced zones adjacent to the main piers. This spatial correspondence does not, by itself, prove a seismic origin for the observed damage, but it strengthens the interpretation that strong ground motion is a plausible contributor to the structural transformations recorded in the monument.

The transition from the lower initial configuration represented by Model 3.0 to the higher two-level configuration represented by Model 3.1 should not necessarily be interpreted as a direct response to earthquake damage. It may also reflect hydraulic or urban requirements, such as the need to increase or maintain water supply to a growing city. By contrast, the later modifications represented in Model 3.2, including arch infills and buttress-like reinforcements, are more directly compatible with structural repair and stabilization.

The comparison between Models 3.0, 3.1, and 3.2 also shows that the repaired state of the monument is mechanically more stable than the unrepaired one. In the modelling, the addition of arch infills and buttress-like reinforcements reduces stress concentrations and improves the overall behaviour of the structure. This result is important because it suggests that the preserved repairs were not only architectural modifications, but also interventions compatible with an increase in structural stability. The correspondence between repaired zones and mechanically sensitive sectors supports the idea that the builders or restorers were responding to real structural weaknesses, whatever their original triggering cause may have been.

More broadly, the Antioch case suggests that the present state of preservation records successive adjustments to zones of recurrent mechanical weakness. The model results do not establish a seismic trigger by themselves, but they provide a mechanically consistent framework linking collapse, repair, and reinforcement to stress localization within the monument. This supports the interpretation of the aqueduct as a structure whose architectural evolution was shaped, at least in part, by repeated mechanical constraints.

5.4. Architectural Implications and Limits of the Approach

Taken together, the generic models and the Antioch case study suggest that some recurring Roman architectural solutions were mechanically sound in seismic environments. Each aqueduct nevertheless reflects a site-specific compromise between topography, source elevation, available raw materials, and the hydraulic constraints imposed by the route. Narrow arches, repeated superposed openings, reinforced piers, and the use of relatively light materials all tend to reduce the severity of the structural response in the models. Even if Roman builders did not formalize these choices in modern engineering terms, the preserved structures indicate that empirical construction knowledge could lead to stable solutions consistent with the criteria highlighted here. More broadly, archaeological evidence from the ancient Mediterranean suggests that empirical antiseismic construction practices did exist in antiquity, although their identification in individual monuments remains context-dependent [49]. In that sense, the Antioch aqueduct may be viewed as a particularly good illustration of Roman engineering practice under strong environmental constraints.

At the same time, the present modelling remains a first-order approach and should be interpreted accordingly. The masonry is homogenized, the analyses are linear, the site response is not modelled explicitly, and some aspects of the Antioch geometry remain reconstructed rather than directly observed. The results are therefore best understood as identifying mechanically plausible weak zones and relative stability trends, rather than as reproducing the full complexity of damage initiation and propagation in real masonry monuments. Used in this way, the models can serve as a useful diagnostic tool for archaeoseismology, especially when combined with field observations, construction history, and historical evidence. Recent post-earthquake surveys of monumental masonry buildings in Hatay and Osmaniye similarly show that wall geometry, construction quality, and mortar properties strongly condition damage patterns, underlining the need to integrate such parameters more explicitly in future monument-specific analyses [50].

6. Conclusions

This study provides a first-order assessment of the seismic vulnerability of Roman aqueduct bridge structures by combining simplified finite-element modelling of generic arch-and-pier configurations with application to the Antioch-on-the-Orontes aqueduct at Harbiye. Across the generic models, the results consistently identify two recurrent weak zones under seismic loading: the arch springings and the bases of the piers. Structural response is strongly influenced by geometry and material properties. Higher piers, wider arches, thicker decks, lower stiffness, and greater mass all increase displacement amplitudes and stress concentrations, whereas buttresses, larger support piers, and lighter construction materials improve structural stability. These results highlight a set of first-order architectural parameters that control the sensitivity of aqueduct bridges to earthquake shaking.

The application to the Antioch aqueduct shows that the modelled stress distributions are spatially consistent with several observed damage and repair zones in the standing monument. In particular, the unrepaired configuration produces the highest stress concentrations in sectors that are now damaged, collapsed, or later reinforced, whereas the repaired configuration shows a clear reduction in stress concentrations and a more stable overall response. Although such correspondence does not in itself prove a seismic origin for all observed damage, it supports the interpretation that seismic shaking plausibly contributed to the long-term structural evolution of the monument and that the preserved repairs were mechanically effective responses to structural weakness.

More broadly, this work shows that numerical modelling can provide a useful diagnostic framework for evaluating the seismic vulnerability of historical masonry infrastructures in tectonically active regions. In the case of Antioch/Antakya, this question has acquired renewed relevance in light of the destructive 2023 earthquake sequence, which once again highlighted the extreme seismic exposure of the region. The approach remains intentionally simplified, but it demonstrates that first-order modelling can help identify mechanically plausible weak zones, compare alternative structural states, and support archaeoseismological interpretation. Future work should extend this framework through nonlinear analyses, refined site-specific geometry, and scenario-based ground motions tailored to the local tectonic setting.