Testing the earthquake damage and vulnerability of the Cherichira aqueduct bridge, Kairouan (Tunisia) with discrete element modeling

The Cherichira aqueduct, originating from Roman times, supplied the city of Kairouan, Tunisia, with water and has had alternating phases of damage and repair after the Roman and during the Aghlabid and Fatimid era. A crucial section of the lifeline is the Cherichira aqueduct bridge (CAB), and scholars have discussed the possibility that earthquake ground motions caused damage which disrupted the water supply of Kairouan. However, little was known about the dynamic behavior of the bridge and its vulnerability to earthquake ground motions. Computer-aided design based on a detailed laser scan model of the remains of the bridge and published data were used to reconstruct the CAB of the Aghlabid period. Subsequently converted into a discrete element model, the digital version of the CAB was subjected to analytic ground motion signals and full 3D simulations of local earthquakes. The CAB model shows a fundamental eigenfrequency close to 1 Hz in the direction transverse to its trend, and single-component ground motions in this direction with peak particle velocities above 1.0 m/s cause damage to the top of the CAB. Among the earthquake scenarios with full 3D ground motions applied, only the activation of a nearby thrust fault caused distinct damage. While fractures in the ruins of the CAB cutting through the upper part of the bridge which includes the water canal are a likely cause for disrupting the water flow and are similar to the damage pattern produced in the model calculations, a solely seismogenic cause of the total collapse of some parts of the CAB cannot be verified by the simulations.


Introduction
During Roman times, an aqueduct located west of the present-day city of Kairouan was constructed to supply agricultural fields and nearby villages with fresh water from the Bir el Adine mountains (Fig. 1). A neuralgic section of the aqueduct is a bridge spanning the Oued el Mouta valley with a length of 180 m and a maximum height of 14 m, the Cherichira aqueduct bridge (CAB). The present state of this bridge raises questions about its fate in history as some sections are totally collapsed and/or the material missing, while a substantial part including four pillars and the corresponding arches with the superstructure and water canal are still standing. However, this section also shows damage to a degree which implies a dysfunction of the lifeline. Large vertical cracks cut through the water canal and parts of the superstructure of the pillars and tufa deposits at these locations indicate a prolonged period of heavy leakage. The bridge shows at least two phases of degradation and repair in the Aghlabid and Fatimid period, and earthquakes have been proposed as a potential cause for the damage (e.g., Bahrouni et al. 2014Bahrouni et al. , 2020. In a recent paper Hinzen et al. (2021), in the following referred to HI21, presented the results of two field campaigns at the CAB including laser scanning the ruin of the bridge, dating of mortar samples, and horizontal-to-vertical spectral ratio (HVSR) measurements to quantify local site conditions. In addition to a proposed timeline of events, they presented a seismotectonic model of the area (Fig. 1) and used 18 earthquake scenarios to calculate synthetic ground motions at the CAB site and estimate macroseismic 1 3 intensities. That work suggests using a 3D discrete element model (DEM) of the CAB to examine its dynamic behavior and to test the hypothesis of a seismogenic origin of the observed damage; this is the main goal of this contribution.
Computer models using DEMs (Cundall 1971) which have been used in archaeoseismology in various types of applications have an advantage (over finite element models) in that the motion of building material, modeled as individual blocks with varying degrees of stiffness, can be traced when shifted, rotated or toppled after they broke free from the structure (e.g., Hinzen et al. 2013;Psycharis 2015;Lemos et al. 2015;Schweppe et al. 2017;Hinzen and Montabert 2017). This feature is of particular importance when the pattern of fallen blocks in a ruined building may be the primary evidence of past earthquakes. In most cases, ground motion parameters must be deduced from fallen classical columns or deformed and toppled walls (Hinzen 2009;Ambraseys and Psycharis 2012;Stefanou et al. 2015;Papadopoulos et al. 2019). Favorable conditions for archaeoseismological DEMs are structures where worked ashlars with flat contact surfaces are stacked without any mortar or clamps (Hinzen et al. 2013;Psycharis 2015). In these cases, only the friction between blocks (which are usually from the same material) and classical mechanics (Goldstein 2002) determines the dynamic behavior of the structure. However, such models have also been applied to structures with mortar binding between blocks (e.g., Lemos 2007;Sarhosis et al. 2016Sarhosis et al. , 2018. A DEM (Cundal and Strack 1979) is based on the repeated performance of three steps, with (1) detection of contacts between blocks or particles, (2) calculation of interaction forces, and (3) numerical time integration to solve the interaction of the elements (Bićanić 2004). The technique has been applied previously to different types of structures. For example, Clementi et al. (2020) presented a detailed study of a clocktower in Amatrice, Italy which was damaged during the 2016 earthquake sequence using a DEM and the Non-Smooth Contact Dynamics Method (Clementi et al. 2019). Tóth et al. (2009) and Rafiee et al. (2013) used DEMs to study the effect of the backfill on the mechanical behavior of a multi-span masonry arch and the mechanical behavior of a stone masonry bridge, respectively.  (upper left) shows the location of the Kairouan and CAB area (black box). Black lines show the major fault zones. Red lines are the surface traces of the assumed fault segments. The red crosses show the reference points for the sources and the yellow shaded areas indicate the surface projection of the idealized fault plane geometries. Circles numbered 1, 2, and 3 indicate the suggested epicenters for three rupture processes in each source. Black squares show location of villages possibly damaged by the earthquake scenario. The blue line indicates the trend of the Cherichira aqueduct and the blue diamond the location of the CAB. The digital terrain model is from ALOS World 3D (Takaku et al. 2014) Ancient aqueducts have also been the focus of 3D dynamic modeling: Volant et al. (2008) made an interdisciplinary study of the seismic behavior of the Pont du Gard Roman aqueduct bridge (Fabre et al. 2000) near Nîmes in southern France using FEMs. They discretized the threelevel arches of the aqueduct by 3D linear elastic triangular elements and showed that the fundamental eigenfrequency of the bridge is around 0.4 Hz. This low frequency prevented severe damage from the model earthquake of magnitude 6 at 10 km distance. Alternatively, they showed, that the smaller bridge Pont de la Lône with an eigenfrequency of 7 Hz for the same earthquake suffered from significant tensile stress at the base of the canal where failure was observed.
In this study, we start with the laser scan model of the CAB from HI21 and reconstruct the probable shape of the bridge with the Aghlabid (IX -X century AD) water canal before it was presumably damaged by the 859 earthquake and based on this construct a CAD (computer aided design) model. The latter is transformed into a DEM which is loaded with (1) analytic ground motions to estimate basic dynamic performance, eigenfrequencies and eigenforms and (2) with the synthetic seismograms from HI21 to evaluate the likeliness of earthquake ground motions as the cause of the observed damage.

Setting of the Cherichira aqueduct bridge
The aqueduct of which the CAB is a crucial part brought water from the Bir el Adine mountains west of Kairouan to the city. First built during Roman times, it was repeatedly destroyed and reconstructed (HI21 and references therein). Figure 1 shows the main potentially seismogenic faults which formed the basis for modeling 18 earthquake scenarios. The scenarios comprise activation of the Cherichira (CHE) and Cululis (CUL) strike-slip faults, the Qasar Al Ma (QAM) thrust fault and the El Qarn thrust fault, which presupposes the thrust faults as the main seismogenic sources (HI21). For the latter three scenarios, activation of the SW section (ELQ1), the NE section (ELQ2) and an activation of both sections in one event (ELQT) were included. For each of the six sources, three rupture scenarios with unilateral rupture spreading from each end and bilateral spreading starting in the middle of the fault sections were tested; the rupture scenarios are labeled 01, 02 for the unilateral rupture spreading towards and away from the CAB, respectively, and 03 for the bilateral rupture directions. Table 1 summarizes location and size parameters; more details are given in HI21.

Implementing the discrete element model
In several previous archaeoseismological DEM studies, the models replicated the buildings stone-by-stone with the measures coming from detailed classical surveys and/ or laser scans of the original structures (e.g., Psycharis 2007, Ambrasey and Psycharis 2012, Hinzen et al. 2013, Schweppe et al. 2017. Initially, the outside appearance of some sections of the CAB (particularly for pillars) appeared to be built regularly with worked stones arranged in layers. However, closer inspection revealed irregular construction for most parts of the bridge. The photos of the inner structure of the bridge, including the pillars in Fig. 2 and canals in Fig. 3 show mostly rounded or edgy field stones of mixed size imbedded in a flaky mortar. Thus, a stone-by-stone approach to modeling the CAB was considered infeasible.
Recently, there has been an increasing use of DEM to study cemented materials (e.g., Obermayr et al. 2013;Brown et al. 2014). Radia et al. (2019) showed that DEM can be used to predict the elastic and fracture behavior of brick and mortar materials and capture non-continuous phenomena such as multi-cracking. In this case, a realistic representation of the discrete inhomogeneous structure largely depends on the implementation of an inter-particle bonded contact model.
We implemented a DEM of the CAB in three main steps based on particles representing the stones and a bonding model for the mortar. For ease of calculation, spheres with a diameter of 0.2 m were used as particles. The implemented rolling-friction model in the DEM causes a non-spherical roll motion with a better computing performance than for complex particle forms. As shown in Fig. 4, the laser scan model (HI21) was the basis for an innovative reconstruction of a CAD model of the CAB using Fusion360 (see Data and Resources). First, the scan of the surrounding terrain was used to add a simplified subsurface as shown in 4b. In the reconstruction we followed the draft of the center of the bridge given by Gauckler (1900) and Mahfoudh et al. (2004) where the pillars rest on breaking edges and four arches are flanked by massive abutments of which only the northern one presently exists. As mentioned by HI21, there is still some uncertainty about the exact structure of the missing sections C and D of the CAB (Fig. 5); e.g., whether eventually another pillar existed further north. This uncertainty cannot be solved at this point without new archaeological excavations.
Once the CAD model of the bridge was completed, a box measuring 183 × 17 × 15 m was set around it and the volume of the CAB was subtracted from that block; thus, giving a hollow form of the CAB (Fig. 4d). To transfer the result of the CAD modeling to the DEM software ThreeParticle (3P, see Data and Resources), the form was converted into STL (Standard Triangle Language) format. Once imported to 3P the process of filling the hollow form with particles was started. In order to resemble the random arrangement of the stones in the bridge, systematic packing of the spheres, e.g., cubic closest packing, was not an option. Instead, sections of the model were filled with spheres in a random arrangement and calculations were made that let the spheres freely settle only under the action of gravity. To avoid cavities, the model was shaken part of the CAB with blue particles; g to i filling the rest of the CAB in layers with red particles; k after activating the bonding between particles, the hollow form is removed 1 3 horizontally during the settlement of the spheres. During the fill phase, no bonding between the spheres was active and a low coefficient of friction (µ = 0.1) was used. To allow different material parameters for the Roman and the Aghlabid parts of the bridge, the central (Roman) part with the pillars and arches was filled first (Fig. 4f) and then the rest of the bridge followed in layers with a different particle type until the model was filled all the way up including the sides of the water canal. The total number of spheres in the model is 263,174. After completion of the filling, a total of 587,969 bonding elements were activated and the hollow form inactivated. The implemented bonding model is an extension of the classical Timoshenko beam theory (1922) and particle pairs are connected by so-called bond elements; e.g., Obermayr et al. (2013) used this to model cemented sand particles using a linear Timoshenko beam to realize an elastic bond between the centers of two Particles. This model seems suitable for the mortar bond as it allows the bonds to withstand higher compressive than tensile stresses and uses an accidental failure of individual bonds, which is controlled by a random strength reduction factor.
The trend of the CAB section is N331°E and for the DEM the y-direction is parallel to this trend, and thus, the x-direction transverse to the bridge with the z-axis pointing upward. The model building steps are further explained in an animation in the Online Resource 1 and Online Resources 2-4 show examples of filling the hollow form with particles.
In the study of ancient monuments, ambient vibration measurements and structural health monitoring has become an established technique in the past decades (e.g., El-Borgi et al. 2008;Psycharis 2015;Salonikios and Morfidis 2018;Montabert et al. 2022). Reviews of the technique are given by, e.g., Carden and Fanning (2004), Doebling et al. (1996), and Gandham (2021). The concept is to measure the basic eigenfrequencies and eigenmodes of a structure through ambient vibrations (or controlled excitation), often combined with finite element modeling where the model parameters are adjusted, so that the model characteristics agree with the measurement results. Long-term surveillance can reveal changes to eigenfrequencies and, thus, indicate environmental influence and sudden or progressive deterioration (Hinzen et al. 2012;Combey et al. 2022). Another application is to explore the effect of retrofitting devices to strengthen the structures (e.g., Chrysostomou et al. 2008). The CAB is essentially the remains of a monument which is (1) partially collapsed, (2) severely fractured, and (3) showing traces of continued (natural and/or anthropogenic) deterioration. Thus, monitoring the current status would reveal eigenmotions significantly different from those of the intact Sensors are consecutively numbered from left to right (N to S) and bottom to top. The insert shows a detail of the sensor placement at the center pillar building (Prabhu 2011) and would therefore be of only limited use in the modeling process of the complete CAB. Another practical aspect is that there are severe safety concerns to access the top of the CAB for monitoring purposes. The HVSR measurement taken by our previous field campaign (HI21) at the base of the northern most pillar showed a clear peak at 0.87 Hz, this might be related to the basic eigenfrequency of the ruin in its current state in this section.
In view of the current state of preservation, we subdivided the structure of the CAB into nine sections labeled A to I (Fig. 5) which are described in HI21 in detail, and which will be adopted in the following. Figure 5 also shows the location of virtual sensors, each measuring 0.5 × 0.5 m in the y-and z-direction and covering the whole width of the CAB model in the x-direction. These sensors are arranged in five vertical profiles through the pillars and one horizontal profile along the top of the CAB just below the water canal (Fig. 5). Parameters (like the average velocity of all particles within the sensor) can be extracted during the postprocessing for all stored time steps.

Material parameters
Laboratory testing of building material from the CAB could not yet be done and standard values were applied for the limestones and the bonding elements. The density, shear modulus, and Poisson's ratio for the particles were set to 2.7 Mg/m 3 , 24 GPa, and 0.28, respectively. Static and rolling friction between particles is 0.7 and 0.6, respectively. Young's modulus and Poisson's ratio for the bonding was set to 10 GPa and 0.25, respectively. The Mises Stress of the bonding elements, representing the tensile strength of the mortar was set to 150 Pa for the particles in the lower Roman part of the CAB and due to the reduced building quality to 125 Pa for the Aghlabid part of the CAB (e.g., Maroušková and Kubát 2017). For a better representation of the accidental failure of individual bonds the strength of the bonding is calculated as a random number in an interval between 80 and 110% of the nominal values.

Calculations
The displacement of selected particles was output for each time step during the calculations through a user defined Application Programming Interface (API), so that the movements of interesting parts of the CAB could be monitored. In addition, the 3P software offers several tools for post-calculation analysis of the results. For all calculated cases, a video of the CAB behavior (Online Resources 05-19) and a series of still photos for further analysis was produced. The color of the particles of the DEM can be bound to instantaneous values of numerous parameters. Of these, we used the instantaneous velocity as a standard inspection tool for all calculations sometimes bonding or Mises stress. Another tool is the placement of the sensors within the model to select monitoring locations for various parameter (Fig. 5).
The critical time step for the integration process for the CAB model is 1.0083 10 -4 s. All calculations were done with a fixed time step of 2.166 10 -5 s (20% of the critical time step). For a typical 30-s time window of a seismogram, this results in 1,385,130 iterations for the 263,000 particles. On a workstation with a 16 core AMD Ryzen 9 processor running at 4.15 GHz, a 30-s real-time calculation consumes about 170 h CPU time. Calculation results were stored at 0.025-s intervals with a file size of 470 Mbyte and a total of 535 GB for a 30-s calculation.

Application of analytic signals
Two types of analytic ground motion signals were applied to the CAB model to test for eigenfrequencies, eigenmodes, damping, and amplitude thresholds above which structural damage occurs. The first signal type is a cycloidal pulse (CP) (Jacobsen and Ayre 1958;Anooshehpoor et al. 1999). This trigonometric one-sine pulse simulates the fault parallel movement of the ground next to a strike-slip fault Roussos 1998, 2000;Hinzen and Montabert 2017;Fig. 6 Basic forms of a one-sign cycloidal pulse and Morlet wavelet as ground displacement, velocity and acceleration Schweppe et al. 2017). The displacement is monodirectional with a final static displacement, and the time from zero to the final displacement determines the maximum velocity and the acceleration of the start and stop phase (Fig. 6).
The second signal is a Morlet wavelet (MW) (Goupillaud et al. 1984), a complex exponential carrier multiplied by a Gaussian envelope. The basic frequency and the number of cycles within a given time window is determined by a constant. The displacement is bidirectional and the rise and fall of the amplitudes of the repeated oscillations has similarities to surface waves of local earthquakes. The basic form of the signals is shown in Fig. 6 and all MWs were applied in a 10-s time window with the peak velocity amplitude occurring at 5 s.

Cycloidal pulse
A first series of CP tests was calculated with pulse frequencies of 4.0, 2.0, 1.0, and 0.6 Hz, all at a peak ground velocity (PGV) of 1.0 m/s with the ground motion directed along the positive x-axis (transverse to the trend of the CAB). Figure 7 shows snapshots of the velocity of the particles along the eastern side of the bridge. Even though the PGV is the same for the three tests, the decrease of frequency, concomitant with an increase of peak ground displacement (PGD), causes a stronger velocity response of the CAB. While the motion within the CAB desists after about 1.0 s in case of the 4 Hz pulse, reverberations last beyond 2.5 s in case of the 1 Hz and 0.6 Hz pulses. In Fig. 7, the blue and red colors indicate instantaneous velocities in opposite directions transverse to the bridge trend. The tests with the lower frequencies (1.0 and 0.6 Hz) show that in the initial phase, the whole structure moves in the positive x-direction indicated by the red colors between 0.3 s and 1.0 s and 0.5 s and 1.3 s along the timeline, respectively. However, at the transition to the first backswing, visible between 1.15 and 1.6 s in Fig. 7, the motion is no longer coherent along the bridge. Due to the U-shaped valley forming the subsurface and thus the changing height of the CAB, the backswing starts in the center, the highest part of the bridge indicated by the blue color at the top around 1.5 s-1.7 s and then spreads in an undulating pattern towards the sides of the structure. The dependency of the CAB response on the pulse frequency is also visible in Fig. 8   As expected, in all cases the center of the bridge produces the largest velocity amplitudes with the highest values in case of the 2 Hz pulse which also produced larger velocities along the northern and southern end in sections F, H and A of the CAB where the height decreases to less than 2 m above ground. While the first motion is coherent along the bridge, the peaks and troughs of the velocity signal during the backswing (particularly in case of the 1.0 Hz CP test) illustrate the loops and nodes of the undulating motion. The above-mentioned extended duration of the CAB response for the lower frequencies is also obvious in this graph. Animations of the velocity response of the CAB in these four CP tests are available as Online Resources 5-8.
While Fig. 8 illustrates the instantaneous velocities with respect to time, Fig. 9 summarizes the maximum velocities from each sensor along the top of the bridge (Fig. 9b) and for the five vertical profiles through the pillars (Fig. 9c). For the CP with 0.6 and 1.0 Hz, the vertical profiles show a smooth and constant increase of the maximum velocity from bottom to top along the pillars, indicating that the motion is dominated by the basic eigenmode. Values start at the bottom with the PGV of the ground motion of 1.0 m/s and reach values at the top of 1.49 to 1.86 m/s and 2.02 to 2.41 m/s for 0.6 Hz and 1.0 Hz, respectively. The massive abutments (V1 and V5) at the ends show smaller values and the center pillars V2 and V3 larger maximum values. At the higher CP frequencies, the velocity increase with height is no longer monotone. At 2 Hz, the velocity decreases from 1.0 m/s at the base to 0.71 m/s at a height of 8.5 m due to stronger inertial forces and then increases to 2.1 to 2.7 m/s at the top. At 4 Hz the initial velocity decrease is even greater (0.64 m/s at a height of 7.5 m) and the values at the top are from 1.8 to 2.9 m/s with a significant change above 12 m height where the top part with the water canal starts. The no longer smooth and parallel curves of the vertical profiles show that at higher frequencies the motion of the pillars is dominated by higher modes.
This trend of a more irregular motion with an increase in frequency is also corroborated by the horizontal sensor profiles (Fig. 9b). With sensors of 0.5 m width placed next to each other, sharp velocity discontinuities along the trend of the bridge indicate inelastic deformation, the loss of bonding between particles, and thus the formation of cracks. At 0.6 Hz, the ends of the bridge experience a maximum velocity close to the PGV. In the center between a 30 m distance from the southern end to 125 m (sections B, C, D, and E Fig. 5) three peaks of maximum velocity occur with values of 1.6 to 1.9 m/s. However, the curve is smooth with no indication of major disruptions. At 1.0 Hz the response is similar but the range of amplitudes larger than the PGV and has broadened at the northern section to 153 m (reaching well into section A) and the maximum velocity in the center of the bridge is 2.4 m/s. Around 95 m profile distance, disruptions in the curve are the result of loss of bonding and thus the formation of cracks cutting the top of the CAB in the transverse direction. The increase of the pulse frequency to 2 Hz changes the trend of the maximum velocity curve along the top in several ways: (1) at the southern end (0-45 m, sections F and G) and north of the pillars (108-155 m, section B and part of A), a velocity of 1.5 m/s is reached, indicating a stronger response of these shallow parts of the bridge; (2) two peaks in the maximum velocity curve occur north and south of the center of the pillars with 2.6-2.9 m/s and (3) in the center there is a local minimum of 1.9 m/s, a consequence of a different wavy motion at the top of the CAB. Several disruptions indicate the formation of cracks. At the highest tested pulse frequency of 4.0 Hz, such disruptions are even more frequent, particularly in the section above the arches; however, due to the frequent cracks, the peak velocity values are smaller than for a 2.0 Hz pulse. In addition to the sensor recorded velocities, the displacement of selected single particles was monitored during the calculations. Figure 10 shows the displacement of an observation point just below the upper (Aghlabid) water canal in the middle of the bridge for the four different CP frequencies. With the 1.0 Hz CP test, the bridge responds with its fundamental eigenfrequency in the transverse direction of 0.94 Hz, and the damping is 9.4%. A further test with the same CP with a PGV of 1.0 m/s in the bridge-parallel direction shows a higher eigenfrequency and stronger damping of 2.0 Hz and 24%, respectively.
In the next set of tests with a CP of 1.0 Hz, the PGV was increased from 1.0 to 1.5, and 2.0 m/s and the development of the instantaneous velocities within the horizontal sensor array is shown in Fig. 11. The undulating of the top of the CAB increases with the increasing strength of the pulses.
The plot of the maximum velocity amplitudes in the vertical sensor arrays (Fig. 12c) shows that at a PGV of 1.0 m/s the basic eigenmode in the transverse direction dominates the response of the CAB. At 1.5 m/s and more pronounced at 2.0 m/s PGV, the influence of higher modes is visible. The maximum velocities observed at the top of the bridge for the three tests are 2.4, 3.6, and 4.7 m/s, thus the ratio between the maximum velocity at the top of the CAB the PGV for the simple CPs is 2.4 in all three cases. In both Figs. 11 and

Morlet wavelets
A series of tests with MWs was calculated to further explore the structural limits of the CBA model. In contrast to the single CP the MWs generate repeated quasi harmonic motions.
The ground motions were again applied in the bridge-transverse direction only. The basic frequency was set to 1.0 Hz and the PGV increased from 1.0 m/s to 2.0 m/s. The instantaneous velocities in the horizontal sensor array for the MW tests are shown in Fig. 13. While the top of the CAB moves coherent during the first 3 to 4 oscillations, with an increase of the vibration amplitudes, the undulation along the bridge increases severely and with an increase in PGV, additional cracks develop. The MW test with 2.0 m/s PGV causes partial collapse in the center of the bridge, affecting the upper Aghlabid section. Two such sections are indicated The plot of the maximum particle velocities measured with the sensors (Fig. 14) shows that the undulation of the motion at the top of CAB is more severe than in case of the CP tests (Fig. 12). While the CP produced a clear increase in the maximum velocities along the whole length of the CAB with an increased PGV, the MW tests do not produce the same results. In the distance range from 40 to 140 m (sections B to E) the maximum velocity curves are close to each other. With the exception of the sections between 77-85 m and 88-94 m where the top collapsed, the maximum particle velocities of the test with 2.0 m/s PGV are even below those for the 1.0 m/s PGV. This is an effect of an increased number of broken bondings and a shift from elastic to inelastic response. The vertical sensor profiles (Fig. 14c) support this interpretation; the maximum particle velocity at the pillar tops is 3.58 m/s in case of the 1.0 m/s PGV and 3.91 in case of 1.5 m/s PGV. For the strongest MW with 2.0 m/s PGV, the maximum particle velocity is reduced to 3.34 m/s (not considering the high velocity above 4.0 m/s caused by parts than fell off). The strong MW with 2.0 m/s PGV causes several cracks in the top and also in the shallower sections of the bridge. Animations of the velocity response of the CAB to the three MW tests are available as Online Resources 11-13. Figure 15a shows velocity snapshots taken at 5.0 and 5.8 s for a test with 1.0 m/s PGV further illustrating the undulating nature of the motion along the length and height of the bridge. A snapshot of the distribution of the kinetic energy within the model in Fig. 15b shows hotspots of energy (black arrows) along the bridge and a pattern which agrees The smaller graphs on the right side show the sensor-measured velocities of the central part of the CAB; white arrows indicate fractures that developed in the top of the bridge and red arrows mark collapsed sections in general with the observed major cracks and spilling sections in the ruin (HI21) as shown in Fig. 15c.

Application of earthquake ground motions
Upon completion of test calculations with analytic signals, synthetic ground motions of the earthquake scenarios described in HI21 were used as boundary conditions in the DEM of the CAB. The complete set of seismogram plots can be found in the electronic supplements of HI21 (https:// static-conte nt. sprin ger. com/ esm/ art% 3A10. 1007% 2Fs42 990-021-00062-9/ Media Objec ts/ 42990_ 2021_ 62_ MOESM1_ ESM. pdf, last accessed 05-2022).
The horizontal components were rotated into bridgeparallel (P) and bridge-transverse (T) directions, and the elastic response spectra for these rotated ground motions were calculated for a damping of 9.4% and 24% for the T, and P component, respectively, based on the results of tests with the CPs (Fig. 16). For all six earthquake sources (Fig. 1), the model with a unilateral rupture directed towards the CAB site (mechanism 01) causes the largest response amplitudes, followed by the bilateral rupture models (03), and a smaller response for the ruptures Among the 18 scenarios, only three caused significant damage to the model: During ELQ1/01 (earthquake/rupture scenario) cracks developed in the upper part of the bridge in the same manner as seen in the tests with analytic signals and some smaller parts of the water channel separated and fell down. The two strongest ground motions of scenarios ELQT/01 and ELQT/03 caused a partial collapse of the upper section of the CAB model. However, for the scenario where the rupture progressed away from the CAB (ELQT/02), the significantly lower velocity response amplitudes above a period of 0.5 s (Fig. 16) did not cause collapses. Figure 17 shows the instantaneous velocities from the 364 sensors of the horizontal array along a time axis of 30 s for the six earthquake scenarios of rupture mechanism 01. The z-axis of the plots, representing the velocity, is the same for all six tests, but the color scale of the 3D field was correlated to the maximum amplitudes to better highlight the individual values. The time axes start at source time except for ELQ2 and ELQT which start at 10 s, due to the larger distance to the epicenter (Fig. 1). Animations of the velocity response of the CAB to the six earthquake scenarios are available as Online Resources 14-19. Figure 18 summarizes the maximum velocities observed with the sensors within the CAB model for the five vertical sensor arrays in the pillars and the horizontal array along the top of the model. For CHE/01, CUL/01, and ELQ2/01 maximum velocities even in the center to of the CAB are well below 1 m/s. QAM/01 reaches 1.1 m/s above the pillars,

Discussion
In engineering seismology, many different ground motion intensity measures (IM) have been introduced often including modeling of collapse structures (e.g., Ambraseys and Psycharis, 2011;Buratti 2012). In addition to peak values of displacement, velocity, and acceleration (PGD, PGV, and PGA), the elastic spectral acceleration, S a , at the fundamental period of a structure is a useful parameter. Bianchini et al. (2009) suggested that the geometric mean of the spectral acceleration over a range of periods, S a,avg , is a better IM than S a or PGA to predict inelastic structural response of buildings during earthquakes. We determined these IMs for the analytic signals as well as all 18 earthquake scenarios for the CAB (Fig. 19) to correlate them with damage observed in the test cases. The cases are categorized as "no damage" with elastic behavior of the CAB model, "cracks" when clear fractures occurred (mainly) in the upper section of the CAB model, and "partial collapse" when parts of the water canal and its substructure and/or arches separated and fell down.
The tests with analytic signals do not indicate a clear damage threshold for either PGA or PGD. But all tests with a PGV above 1.0 m/s show damage with increasing damage level at higher values of the PGV. For S a,avg the threshold is at about 7 m/s 2 for the MW but somewhat lower at 5 m/s 2 for the CPs. The PGV threshold of 1.0 m/s and the S a,avg threshold of 7 m/s 2 also holds for the simulated earthquake signals, of which only three reach the damage range. Of these, ELQ1/01 produced cracks and ELQT/01 and ELQT/03 caused collapse at the upper central part of the CAB (Fig. 17).
The intensity estimates for the 18 earthquake scenarios of HI21 for the site of the CAB ranged from below V to X. The dynamic tests showed damage to the CAB model only for the strongest ground motions (Fig. 16) from these scenarios with estimated site intensities of IX and X. Scenario ELQ1/01 with intensity IX resulted in the cracking of the top of the CAB and downfall of some small parts of the water canal, while scenarios ELQT/01 and ELQT/03 (Mw Fig. 16 Elastic velocity response spectra at 24% damping for the bridge-parallel (P) and 9.4% damping for the bridge-transvers (T) components of 18 earthquake scenarios. Labels in the upper left corner of each plot denotes the earthquake source (see Table 1); line color and the numbers 01, 02, and 03, as explained in the legend, indicate the ground motion component and the earthquake rupture model of Fig. 1 7.2) with estimated intensity X caused partial collapse of the model. Scenario ELQT/02 with the largest Joyner-Boore distance (Joyner and Boore 1979) also resulted in a site intensity of IX; however, here the bridge-parallel component was stronger than the bridge-transverse component and no significant cracking was induced. The intensity estimates in HI21 were based on the PGV and PGA of the simulated ground motions. The eigenfrequency in the more vulnerable bridge-transverse direction at 0.94 Hz is lower than that of residential buildings. Even though the European Macroseismic Scale (EMS, Grünthal 1998) recommends not to use special structures (e.g., lighthouses, bridges, monumental buildings like cathedrals) in the routine determination of intensities, the results of the modeling of the CAB should be discussed in terms of site intensity.
While some special structures in the sense of the EMS tend to be more vulnerable to certain ground motions than ordinary structures, Camelbeeck et al. (2014) showed that during moderate earthquakes in Belgium and Germany churches suffered damage of a higher degree than expected for the local intensity based on ordinary structures. Alternatively, aqueduct of Pont du Gard (Volant et al. 2008) did not suffer damage during local earthquakes mainly due to its low eigenfrequency of 0.4 Hz. Clamenti et al. (2020) showed why the Amatrice clock tower survived the main shock of the 2016 magnitude 6.2 Amatrice earthquake, and Combey et al. (2021) describe a high earthquake resistance of Inca infrastructures.
Most of the CAB building material is rubble and field stones (in a mortar matrix) and would require a low vulnerability class in terms of the EMS, but due to its compactness and dimensions (pillars measure 2.9 × 2.9 m), the intact CAB can be moved to in a higher vulnerability class of C-D. This agrees with the observation that the ground motions of the scenarios on the QAM fault and ELQ1/02 and ELQ1/03 with intensities of VII and VIII, caused no damage in the CAB model larger than grade 2. However, scenario ELQ1/01 with intensity IX did cause damage of grade 3 and the scenarios with intensity X caused damage of grade 4-5. In turn, caution is suggested in the EMS (Grünthal 1998) to when estimating site intensities for special structures. It is highly probable, that simple houses in the surrounding villages were heavily damaged or destroyed during the 859 earthquake as Fig. 17 Results of six tests with 3D earthquake ground motions. Labels at the right of the z-Axis indicate the earthquake scenario (Table 1 and HI21). The 3D plots show the instantaneous velocity measured with an array of sensors along the top of the CAB at 0.5 m intervals (see Fig. 5) with respect to time. White arrows indicate fractures that developed in the top of the bridge and red arrows mark collapsed sections. The individual color scales for the velocity, v, are shown next to each plot described by Bahrouni et al. (2020), while the CAB suffered less damage.
The fractures in the top part of the CAB model comprising damage grades 2-3 (EMS) correlate well with the observed vertical fractures in the CAB and their frequency and distribution along the bridge (Fig. 14). These fractures are the location of previously deposited tufa, clear indicators of continued water spill in cracks. Thus, the model calculations show that earthquakes of magnitude 6.5-6.7 on the ELQ fault (Table 1) could well disrupt the function of the aqueduct. However, these events do not explain a total collapse. The scenario with the strongest ground motions from a magnitude 7.2 earthquake rupturing both sections of the ELQ fault does cause partial collapse of the CAB model, particularly the (Aghlabid) top part in section E including the collapse of the most northern and southern arches, although, it does not explain the total collapse in sections C and D (Fig. 20). Unless the complete CAB in sections C and D was built much less stable than in our current model or there are unknown irregularities in the subsurface, a seismogenic cause of their total collapse seems unlikely.
Alternatively, it is possible that initial fracturing of sections B, F, and also section I by an earthquake was the nucleus for further disintegration by weathering and/or anthropogenic reuse of the building material easily accessible in these sections.  (Table 1 and HI21) for earthquake fault segments (Fig. 1). a Top view and side view (from North) to the model of the CAB, the ground is shown in gray and the bridge in blue; red lines indicate the location of sensors in a horizontal array along the top and five vertical arrays at the center of the CAB. b Maximum velocity recorded with the sensors along the top of the model; frequencies of the cycloidal pulses are given in the legend. c Maximum velocity of the sensors of the five vertical arrays; frequency and PGV of the ground motion is given in the legend of the first graph 1 3

Conclusions
A discrete element model of the Cherichira aqueduct bridge close to Kairouan, Tunisia, further quantifies the decay process of the structure. The eigenfrequencies in the bridge-transverse and parallel directions are in the range of 0.94 Hz and 2.0 Hz, respectively. The model predicts a threshold of about 1.0 m/s of peak ground velocity in the bridge-transverse direction to induce damage in form of  The many uncertainties affecting the modeling of the earthquake response of the CAB, include correct reconstruction of the missing parts, the material parameters, bonding of particles, and details of ground structure interaction. However, the presented model should adequately reproduce the basic dynamic behavior of the bridge. The application of 18 earthquake scenarios based on a seismotectonic model of the area shows that only the strongest simulated events on the El Qarn thrust fault are capable of causing significant damage to the model of the bridge, making this a primary target for future palaeoseismic surveys.
The result of the predicted site intensities of the earthquake scenarios reaching the level of IX-X is not contradictory to the results of the model calculations. A structure like the CAB is different from ordinary houses on which the intensity scales are based, and much less vulnerable. This explains why it is possible that during the same earthquake the CAB was only fractured along its top but residential buildings in the surrounding villages totally collapsed. The damage to the water canal by the posited earthquakes would have been strong enough to interrupt the water supply of Kairouan through the Cherichira aqueduct and might have triggered the subsequent construction of huge cisterns in Kairouan.
The damage patterns in all test cases show first a cracking of the upper section and then a collapse of the water canal with parts of the superstructure and eventually some arches of the CAB. However, the arches and their superstructure are still in place in section E of the CAB and none of the tests caused a total destruction of the sections C and D as it appears in situ. Thus, earthquake ground motions are a plausible explanation of the vertical cracks visible in the upper part of the remains of the CAB, and their pattern agrees well with the simulations. However, as the earthquake scenarios hypothesize a total rupture of the activated fault segments, they constitute a worst-case scenario. The model does not explain the total collapse of the sections C and D while the Roman arches remained largely intact. Based on tests with the DEM of the CAB it is unlikely that earthquake ground motions caused the collapse of these parts of the bridge. To confirm or disprove this hypothesis archaeological excavations of sections C and D are necessary to produce if so new data which might lead to an improved model of the CAB. An alternative hypothesis is (repeated) flash-flooding in the wadi which the bridge crosses as a non-seismic cause, possibly in combination with anthropogenic action. This scenario should be further explored and combined with more physical dating of bridge material to better clarify the sequence of events.