Performance-based seismic assessment of a historical masonry arch bridge: Effect of pulse-like excitations

Seismic analysis of historical masonry bridges is important for authorities in all countries hosting such cultural heritage assets. The masonry arch bridge investigated in this study was built during the Roman period and is on the island of Rhodes, in Greece. Fifteen seismic records were considered and categorized as far-field, pulse-like near-field, and non-pulse-like near-field. The earthquake excitations were scaled to a target spectrum, and nonlinear time-history analyses were performed in the transverse direction. The performance levels were introduced based on the pushover curve, and the post-earthquake damage state of the bridge was examined. According to the results, pulse-like near-field events are more damaging than non-pulse-like near-field ground motions. Additionally the bridge is more vulnerable to far-field excitations than near-field events. Furthermore, the structure will suffer extensive post-earthquake damage and must be retrofitted.


Introduction
Masonry arches are one of the earliest structural forms that have been used for thousands of years as parts of structures such as bridges and cathedrals [1]. It is widely accepted that the use of the arch as a structural form was developed independently in China and the Middle East more than 5000 years ago [2]. Masonry arch bridges still play a crucial role in transportation infrastructure in European countries because they form a significant part of the road and railway bridge inventory [3]. Historical masonry arch bridges must be preserved because of their importance as cultural heritage assets that carry lessons from past generations, in addition to their functional value as important infrastructure. Historical masonry arch bridges have been exposed to material degradation, scouring effects on boundary conditions, and changes in applied loads during their lifetime. In addition, they have been exposed to exceptional natural hazards such as earthquakes and floods. Furthermore, they were not originally built for the heavy loads they are often exposed to currently, and there are serious concerns about their safety for our society today [4].
Numerical modeling is an important part of seismic vulnerability assessment methodology [5,6]. The seismic vulnerability of masonry arch bridges must be assessed, because of the vulnerability of masonry subjected to earthquakes. Several different modeling strategies have been proposed for the simulation of the behavior of masonry arch bridges exposed to seismic events. Finite element (FE) modeling has been widely used to model masonry arch structures [3]. Using one-dimensional frame elements is the simplest approach for modeling masonry arch structures [7,8]. However, the effect of spandrel side walls and backfill soil material cannot be considered, and to tackle this limitation, nonlinear springs or nonlinear truss elements need to be used to simulate the interaction between the backfill, spandrel walls, and adjacent arches [9]. Two-dimensional (2D) FE models have also been used to explicitly model the fill soil and spandrel walls with higher computational effort compared with the one-dimensional modeling approach [10,11].
Finally, as the most detailed and complex method, threedimensional (3D) FE modeling has been proposed, with detailed 3D modeling of all different parts of a bridge and the assignment of nonlinear constitutive relationships to both masonry and fill soil [12,13].
However, the considerable uncertainties related to the fill-structure interactions are possible limitations of the 3D FE modeling approach for masonry arch bridges. Contact interface elements (CIEs) were used comprising zero-thickness plane interface elements with normal and tangential stiffnesses derived based on the equations proposed for discrete element modeling of masonry structures [14,15]. A comprehensive study to investigate the effect of CIEs was performed, which emphasized the significant effect of CIEs on the modal properties as well as on the seismic behavior of a bridge [14]. A strategy was proposed to model a plane CIE with high normal stiffness to avoid the interpenetration of masonry and soil media such that no tension stiffness is considered [4,16,17]. Therefore, the CIE transfers only the compression and shear forces. Normal stiffness can be derived based on the soil−structure interaction equations, and the tangential stiffness is considered as 0.01-0.1 times the calculated normal stiffness [18,19]. However, the Coulomb friction model can also be used for a CIE with zero cohesion and friction coefficients [4,16,17].
Several studies have shown that structures subjected to near-field (NF) seismic events typically exhibit higher levels of damage owing to the strong velocity pulses of the excitations [20][21][22]. The results of linear time-history analysis of a two-span masonry bridge by applying one far-field (FF) and one NF event revealed that, in general, NF events cause more destructive damage than FF events [23]. However, by comparing the seismic response of a single-span masonry bridge subjected to FF and NF seismic events, it was concluded that the bridge is more susceptible to FF events [24]. This conclusion was further confirmed by performing linear seismic analyses on a single-span masonry bridge considering the soil−structure interaction [25]. Furthermore, similar results were obtained by performing a comparative study on a threespan masonry bridge by applying a set of eight FF and NF seismic events [26]. This finding is consistent with the conclusions drawn for nuclear structures [27].
Pulse-like near-field (PL-NF) ground motions are often caused by forward-directivity effects, and several different quantitative methods have been proposed for identifying pulse-like ground motions [28][29][30][31]. Studies have revealed that pulse-like events impose more significant demands on structures than non-pulse-like near-field (NPL-NF) seismic excitations [32][33][34]. However, these records cannot be labeled as aggressive without reference to structural characteristics [35]. Therefore, there is still a gap regarding the need to investigate the influence of PL-NF earthquake excitations on the seismic behavior of masonry arch bridges and compare the effect of FF and NF seismic excitations by applying a large number of events with several different characteristics.

Methodology and case study overview
In this study, the modeling procedure of a case study, the Roman Bridge located on the island of Rhodes, in Greece, was developed. The 3D geometric documentation was conducted using digital images, laser scanners, and total stations. The dimensions were obtained from the 3D product of geometric documentation. Using these data, an FE model was developed.
In the first part of the study, a dynamic modal analysis was performed, and the model was calibrated using an analytical equation to derive the first natural frequency value of masonry bridges.
As the examined bridge structure is located in a high seismicity zone, the second part of the study focused on investigating the seismic behavior of the bridge. To this end, 15 seismic events were selected that were categorized as FF, PL-NF, and NPL-NF. The records of these seismic events were scaled to a target response spectrum, and nonlinear time-history analyses (NLTHAs) were performed by applying the scaled seismic records in the transverse direction to assess the seismic vulnerability of the bridge. Three performance levels were calculated based on the pushover curve of the bridge, and the values of displacement capacity and demand were compared to assess the bridge performance levels.
The Roman Bridge, illustrated in Fig. 1 is located on the island of Rhodes, in Greece. It was built across the Rhodini stream before its outfall into the Mediterranean Sea at the main exit of the city on the east coast of the island. The bridge is 38.85 m long, 8.4 m wide, with a thickness of 0.6 m for the arch and spandrel components. The arch span is 6.4 m, and the height of the bridge is 5.2 m. Bridge building was a key part of the underlying Roman infrastructure [36,37], and the studied stone masonry bridge dates back to the Roman period. It is one of the few existing ancient bridges that still survive in Greece and remains in continuous use today. Later, repairs to the spandrel walls were reported, and temporary wooden scaffolds were built beneath the arches owing to falling stones. The island of Rhodes is in a high seismicity zone with peak ground acceleration (PGA) from 0.35g to 0.55g with 10% probability of exceedance in 50 years, equivalent to a return period of 475 years [38]. Documentation also reveals that several earthquakes have occurred in the island region, some of which are also associated with tsunamis [39,40]. For this purpose, the seismic vulnerability of the Roman Bridge as a cultural heritage asset and infrastructure must be assessed.

Three-dimensional geometric documentation
Geometric documentation is crucial for developing detailed and accurate 3D simulation models [41]. Geodetic, photogrammetric, and laser scanning data acquisition methods were used to provide an accurate 3D model of the bridge. In total, 2576 aerial images from drones and 271 ground digital images were processed using imagebased modeling software to develop the 3D dense point cloud of the case study. The aerial images (using drones) were processed separately from the ground images because they had a lower resolution. Furthermore, 24 scans were performed using 3D laser scanners to provide point clouds to fill the gaps in 3D dense point clouds derived from the digital images. A local reference coordinate system was established using two total stations, and the coordinates of the required points, such as targets for the point cloud registration and ground control points for the orientation of the images, were determined. The workflow for obtaining the final 3D point cloud is illustrated in Fig. 2. Further details regarding this methodology can be found in the literature [42,43]. After providing the final 3D point clouds, the triangular irregular network method was used to represent continuous surfaces. The light 3D model of the bridge shown in Fig. 3(a) was developed by reducing the number of triangular meshes without compromising the surface detail or color. Furthermore, cross sections were produced with different interval distances, as illustrated in Fig. 3(b). A 3D light model was used for visualization, and the dimensions of the bridge were determined from both products [44].

Three-dimensional finite element modeling
The geometry of the bridge was obtained from the 3D models used to develop the 3D FE model with DIANA FEA software [45]. Masonry is a construction material with complex mechanical behavior [46], and in the discrete element modeling approach, masonry units and mortar are modeled separately by means of interface elements to decrease the corresponding uncertainties [47]. For 3D FE modeling of the masonry part of the bridge, the homogenized method was used in which the discretization of the masonry units and mortar was ignored, which is widely used as a simplistic approach for modeling full-scale unreinforced masonry structures [48]. This approach requires less input data and computational effort than the discrete element method [49]. Figure 4(a) illustrates the different parts of the masonry arch bridge comprising spandrel walls, arches, and parapets made of stone masonry and backfill soil. To model the backfillspandrel and backfill-arch interaction, CIEs were used, as shown in Fig. 4(b), as described previously [50]. The effect of soil−structure interaction was not considered, and rigid boundary conditions were set for the models.
The CIE can transfer normal compression and tangential friction with zero tensile strength; thus, the tensile forces are not transferred. Modeling the CIE is recommended not only to represent the actual behavior of the bridge but also to avoid early convergence problems owing to the low stiffness of the backfill material compared with the stiffness of the masonry material [4]. Nevertheless, the higher computational effort and level of input data and the additional need for a skilled analyst to model the CIE can be considered as limitations of modeling masonry arch bridges using the CIE. τ max A high value was considered for the normal and lateral stiffnesses of the CIE to prevent the interpenetration of the two bodies. The CIEs follow the Coulomb friction model in such a way that the two media can carry shear stresses up to a certain magnitude, calculated using Eq. (1), before sliding across one another where c is the cohesion value (considered to be zero), is the normal stress, and is the friction angle of the soil. The friction coefficient (tangent of the friction angle) is recommended to be in the range of 0.4-0.5 [51][52][53].

Material properties
The total strain fixed crack model was applied for the nonlinear modeling of the masonry parts to define the tensile and compressive behaviors of the homogenized stone masonry with a single stress−strain relationship [54]. The area under the exponential softening curve is calculated based on the tension fracture energy ( divided by the definition of the crack bandwidth (h) of an element. For the compression part, the area under the parabolic curve is calculated based on the compression fracture energy ( divided by h to define the stress− strain curve of the masonry, as shown in Fig. 5. Note that the crack bandwidth model was used in this study, and parameter h is related to the volume of the element [55].
The damage-based shear retention model was selected for the masonry shear behavior. The shear retention was simulated based on the material damage caused by cracking, and the shear capacity was reduced to zero when the masonry was sufficiently damaged. The bridge is made of "sfouggaria stone", which was a common construction material on the island of Rhodes in the past, with a compressive strength ( estimated at 9 MPa [56]. Moreover, the compressive strength of mortar ( is considered to be 1 MPa by choosing a soft mortar type [57], owing to the presence of mildly leached mortar that can be raked out. The recommended equations and corresponding values for the material properties of the bridge masonry used in this study are listed in Table 1.
Moreover, the Mohr−Coulomb material model was selected for the simulation of the backfill soil. The corresponding material properties are presented in Table 2 [58]. The friction coefficient value of the CIE was considered to be 0.4.

Finite element mesh
A preliminary analysis to assess the optimal mesh dimensions was performed on the FE model, and an adaptive mesh size of 0.5 m was used. Moreover, to simulate deformation across the arches accurately enough, a mesh size of 0.3 m was used for the arches. The 8-node quadratic element was selected for the FE modeling, which is the dominant type of element for such simulations. Tetrahedral volume elements were also used in the model to fill parts of the geometry. The FE mesh of the bridge is shown in Fig. 6. The model with the CIE is composed of 12699 mesh elements, of which 1936 are CIE elements.

Modal analysis and validation
Modal analysis was performed to derive the first five natural frequencies and corresponding mode shapes of the structure. The direction-dependent participation factors and modal mass percentages were also calculated for each natural mode and are presented in Table 3.
The first five modes with relatively high participation factor values were selected to perform the sensitivity analyses. Figure 7 shows the natural frequency values and corresponding mode shapes with a contour map of the normalized displacement values of the mode vectors of the first five modes of the model. Based on the results in    Table 3 and Fig. 7 it can be reported that the first, third and fourth modes are in the transverse (y) direction, the second mode is in the longitudinal (x) direction, and the fifth mode is in the vertical (z) direction of the bridge. A comprehensive experimental study was conducted to investigate the dynamic characteristics of historical masonry arch bridges by performing an operational modal analysis of eight bridges [59]. Equation (2) was proposed to calculate the first frequency value of masonry arch bridges (y) (Hz) based on the maximum arch span (x) (m). y = −3.935lnx + 16.824. (2) The calculated first natural frequency of the Roman Bridge based on Eq. (2) is 9.46, which is very close to the first natural frequency derived from the modal analysis of the FE models (< 2%). Although Eq. (2) is only based on the span of the bridge by neglecting the effect of material properties and other geometrical characteristics, it shows a good estimation of the first natural frequency value of masonry bridges, as reported in Ref. [60]. However, operational modal analysis based on in situ tests is required to calibrate the model based on the natural frequency values and corresponding mode shapes of the real structure.

Performance-based seismic assessment methodology
In this section, earthquake records are selected and scaled to a target response spectrum. NLTHAs are then performed, and the performance level of the bridge is assessed based on a damage index. Finally, a comparative study is performed to investigate the influence of the type of seismic event on the seismic behavior of the structure.

Earthquake record selection and scaling
In total, 15 seismic events were selected from the PEER NGA-West2 ground-motion database [61]. The flowchart used for selecting seismic events according to different criteria for the three main groups is shown in Fig. 8. Note that the ground type is considered as type B with a shear wave velocity of 360-800 m·s −1 [62,63], and the average of the Campbell and Joyner−Boore fault distances was used as the source-to-site distance (SSD) [64]. PL-NF event records were chosen based on the method presented in the literature [28] with pulse indicator values > 0.85, early pulses in the time history, and a peak ground velocity (PGV) > 30 cm·s −1 . The characteristics of the selected seismic records, including the record sequence numbers (RSNs) in the PEER database, are presented in Table 4.
All of the selected 15 seismic records were scaled to the target response spectrum [62]. For deriving the target response spectrum, the PGA of the island of Rhodes is considered to be 0.47g for soil class B based on the latest seismic hazard map of Europe in the context of the SHARE project [65,66]. The bridge is more susceptible in the transverse direction than in the longitudinal direction because the mode shape of the fundamental frequency is in the transverse direction. Moreover, the Northridge seismic record with an intensity of 6.7 M was scaled to a PGA of 0.8g and applied to the structure in two directions for preliminary seismic analysis. The results showed that the displacement in the transverse direction of the crown point was ten times greater than that in the longitudinal direction, and the crack width (CW) of the model with longitudinal loading was less than that of the model subjected to transverse loading. Therefore, the seismic vulnerability of the bridge was assessed in the transverse direction.
The earthquake records were scaled and matched to the target response spectrum using the SeismoMatch software package [67] in periods from 0.2T 1 -2T 1 , where T 1 is the fundamental period of the bridge in the transverse direction [62]. An improved scaling method was used to adjust the recorded ground motions without baseline correction [68]. Figure 9(a) shows the single-record spectrum after scaling, and Fig. 9(b) illustrates the meanmatched spectrum.
Five indicators were selected for the scaled seismic records, including PGA, PGV, Arias intensity (Ai), significant duration (SD), and specific energy density (SED), which are illustrated in Fig. 10, with average values for each group of records. Note that the Ai and SED are calculated based on: where and are the acceleration and velocity values, respectively, for each time interval, and is the length of the accelerogram. The significant duration is considered as the interval of time in which 5%-95% of the total Ai is accumulated [67]. The average values of the PGA are close to each other, and the average values of the PGV for PL-NF events are closer to the values for FF events and greater than those for the other two types of records. However, the Ai, SED, and SD values of the NF seismic records are less than those of the FF records after scaling.
Scaled seismic event records were applied to the transverse direction of the structure to assess the seismic vulnerability of the bridge, and NLTHAs were performed. Rayleigh damping is used for the models with factors of 3.641 and 0.0006 applied to the mass and stiffness matrices, respectively. Note that the damping factors were calculated by considering a 5% damping ratio for the first and third modes of the bridge models [14,69].

Seismic performance criteria
Because there are no performance limit states for  masonry arch bridges, comprehensive studies are required to define performance criteria. A methodology for determining the quantitative damage criteria for masonry arch bridges was presented in Ref. [70] based on the pushover curve. Three performance levels were considered, and the displacement limit states were determined based on the criteria listed in Table 5. The relationship between the performance levels and the damage states is shown in Fig. 11. Modal pushover analysis was performed based on the load pattern proportional to the first mode shape [71]. The pushover curves are presented in Fig. 12, by defining the displacement limit states based on the criteria listed in Table 5. The displacement values for the F, LS, and NC limit states are 2.76, 5.78, and 18.5 mm, respectively, as highlighted in Fig. 12. The limit states are derived for the bridge in this study, and a comprehensive study, by developing masonry arch bridges with different geometries, morphologies, and material properties is required to generalize the limit states.

Results and discussion
The results for the maximum displacement of the crown point (Disp) of the bridge and the three aforementioned limit states are shown in Fig. 13(a). The percentage of the cracked volume of the structure relative to the total volume of the masonry is calculated as a cracked volume damage index (CVDI) for each record at the final step of the analysis and illustrated in Fig. 13(b). The CWs in the final step of each seismic analysis are shown in Fig. 13(c). Furthermore, the average (Ave) values of the three aforementioned indices are calculated for FF, PL- NF, NPL-NF, and NF events and all event records, as presented in Fig. 13.
By focusing on the response of the bridge to single records, it can be concluded that different damage indices reflect different structural behaviors. The volume of the cracked parts is relatively high for RSN 828, but the other two indices are not significantly large. By comparing the seismic behavior of the bridge subjected to RSN 3750 and RSN 1083, it can be concluded that a higher value of the displacement corresponds to the point on the pushover curve with 7% of the initial (elastic) stiffness displacement corresponds to 90% of the maximum displacement attained on the pushover curve qualitative description structure is mostly elastic with little or no damage; traffic is not interrupted, and damage can be repaired in a couple of days plasticity starts increasing before and after this performance level; bridge is expected to suffer medium to significant damage; it should still be feasible to repair but cannot be used for a short duration damage is heavy and distributed to the extent that the bridge is near to collapse state; bridge may even be out-of-service or replaced completely Fig. 11 Relationships between performance levels and damage states. Fig. 12 Pushover curve of the bridge with performance limits. displacement does not necessarily result in a higher value of the cracked volume damage index and CW. Furthermore, the displacement damage index cannot necessarily reflect the seismic damage state of the bridge, and performance level criteria should be defined for the cracked volume damage index and verified based on the experimental tests for stone masonry bridges. As can be seen, some events, such as RSN 802 or RSN 879 as PL-NF events and RSN 741 as an NPL-NF event, have destructive effects on the bridge, similar to FF events. Therefore, it should be noted that it is not sufficiently reliable to perform a comparative study on the effect of earthquake excitation characteristics on the seismic behavior of structures by applying one or two pairs of event records.
The results show that the PL-NF events are more destructive than the NPL-NF excitations. Moreover, more extensive damage occurs to the structure subjected to FF events than NF events. The cracked volume damage index and displacement values of the bridge subjected to PL-NF events were 9.6% and 76.6% higher than those of the bridge subjected to NPL-NF events. Therefore, PL-NF events are more destructive in terms of the cracked volume damage index.

Regression and correlation analysis
To investigate the correlation between the excitation properties and seismic behavior of the structure in terms of the three damage indices, a linear regression analysis was performed, and the results, including the linear trendline and R 2 values, are presented in Fig. 14. The positive slope of the trend lines reveals that, generally, all parameters have a positive correlation with the damage indices, however, exceptions can also be seen by comparing the characteristics of some seismic records and the corresponding structural demands.
The R 2 values are presented in Fig. 15 on a scale of 0-0.5 to facilitate the comparison. Higher values of R 2 indicate that the change in parameter can explain the change in structural demand in a better way. Ai can be a good seismic intensity indicator for predicting the seismic behavior of a bridge in terms of all damage indices.

Seismic performance evaluation of the bridge
The bridge subjected to seismic events scaled to the target response spectrum with a 10% probability of exceeding it in 50 years with a 475-year return period passed the acceptance criteria for the LS performance level and reported extensive post-earthquake damage. Therefore, the bridge should be repaired after an earthquake; however, this may not be practical for economic reasons. Furthermore, the bridge should be strengthened to prevent future seismic losses.
The applied accelerograms of the RSN 1083, RSN 1013, and RSN 126 events are illustrated in Fig. 16(a) as samples of the FF, PL-NF, and NPL-NF excitations. Figure 16 selected events. The residual displacements that are presented in Fig. 16(b) are due to the fact that most of the bridge elements were in their plastic phase after applying the selected records. Figure 16(c) shows the magnified deformed shape of the bridge at the end of the analysis time of the RSN 1083 event, which is a common shape in almost all simulations. The results revealed that the relative displacement of the crown points of the two spandrel walls increased in the transverse direction by applying the seismic records. Figure 17 shows the crack pattern and principal crack strain values (E knn ) of the model at the final step of analyzing the bridge subjected to the RSN 1083, RSN 1013, and RSN 126 events. The results show that arches are the most susceptible areas, with higher crack strain values than spandrel walls. Cracks can be found in all the models in the arches and close to the spandrel walls, as shown in Fig. 17 for three-sample earthquake records. The area between the arches is the most vulnerable section of the spandrel wall. Therefore, a strategy to strengthen the arches (especially the area close to the spandrel walls) is needed so that the spandrel wall between the two arches can be sufficiently efficient to reduce the seismic vulnerability of the structure.

Conclusions
In this study, geometric documentation of the Roman Bridge on the island of Rhodes, in Greece, was performed using advanced digital imaging, 3D laser scanning, and total stations, and the methodology for providing an accurate 3D model was described. Because the bridge is located in a high-seismicity zone, the seismic vulnerability of the bridge also needs to be assessed. Fifteen seismic records were therefore selected based on site specifications and categorized as FF, Pl-NF, or NPL-NF.  The event records were scaled to the target spectrum and applied to the FE model in the transverse direction. The results of the NLTHA revealed that the structure is more susceptible to PL events than to NPL events by comparing the three damage indices. Furthermore, FF events are more destructive than NF events. Linear regression analysis revealed that the Ai is the most accurate indicator as one of the seismic records' properties for prediction of the structural demands to different types of excitations. The seismic performance of the bridge was evaluated by defining three performance levels based on the pushover curve. The time-history analysis indicates that the structure passes the acceptability check of the LS performance level and undergoes extensive damage. Temporary wooden scaffolds were built beneath the arches owing to the falling of the stones. Moreover, based on the damage patterns concluded from the numerical analyses, the arches are the most vulnerable structural components, which is consistent with the weakness of the real structure. The observed crack patterns indicate that the structure must be retrofitted by strengthening the arches, focusing on the area close to the spandrel walls and the spandrel wall between the arches. The model should be calibrated based on the operational modal analysis results of the on-site ambient vibration tests, and advanced seismic analyses should be conducted considering soil− structure interaction effects for future studies. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/ by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.