Different Approaches to Assess the Seismic Capacity of Masonry Bridges by Non-linear Static Analysis

A large portion of the existing masonry arch bridges in Italy are still in service in the infrastructure system and are located in a geographical area of high seismic risk. Most of them were built more than 100 years ago taking into account only gravitational loads during the design phase without any seismic analysis. For this reason, a seismic vulnerability assessment has been appointed as mandatory regular activity from the Italian government in order to define the priority of seismic retrofit interventions. In this study, a multi-span slender masonry bridge considered as the most vulnerable typology of masonry bridges to seismic action will be assessed. Then, the results obtained from three different seismic assessment approaches will be discussed and compared. In particular, two approaches based on different FE modeling and the last one built on rigid blocks analysis are considered. Finally, a detailed 3D finite element analysis allowed representing all the collapse mechanisms (global and local) of the bridges are presented. Simplified approaches, even though cannot describe all the collapse mechanisms of the bridges due to seismic action can lead to reliable results.


INTRODUCTION
Masonry arch bridges represent a significant percentage of cultural heritage in Italy and Europe (Melbourne et al., 2007). Most of them were built mainly between the end of the nineteenth century (Brencich and Sabia, 2008) and the early twentieth century, a period characterized by the development of important transportation means such as roads and railways. The increase in traffic loads and recent significant seismic events, combined with the material deterioration, have led to a more in-depth study of the behavior of these structures through the use of different analysis approaches (D'Altri et al., 2019). In many geographical regions with high seismic hazard, the assessment of the seismic vulnerability of existing bridges is required since they have been designed to withstand only gravitational loads. Although seismic actions were not considered in the design, most of the masonry bridges anyhow display good resistance (Da Porto et al., 2016) to horizontal loads due to their strong structural robustness. However, multi-span bridges with slender piers and almost semi-circular arches are still vulnerable to seismic action. The "Claro River Bridge" (Da Porto et al., 2016) is a perfect example of this type of structure, which collapsed after an earthquake occurred off the coast of central Chile on Saturday, February 27th 2010 with a magnitude of 8.8 on the moment magnitude scale.
Besides, many researchers have addressed in recent decades fatigue problems in masonry arch bridges. Amongst the existing methodologies, the stress-life approach and numerical analysis are commonly used by scientists. In particular, the stress-life approach to evaluate the fatigue strength of a specific existing bridge has been used in Laterza et al. (2017) and Casamassima and D'amato (2019). The authors highlighted the effectiveness of this methodology and praised for further studies to derive new curves. On the other hand, the fatigue assessment of arch bridges has been also studied using probabilistic tools. Few authors investigated the remaining service life using the Weibull distribution to calibrate their model (Clark, 1994;Casas, 2009). The authors proposed a new model for the reliability-based assessment of the fatigue performance of existing bridges.
For the purpose of vulnerability assessment, the static and dynamic linear analyses may lead to conservative results, since it is not possible to consider the non-linearity of the masonry material. Non-linear dynamic analysis (Pelà et al., 2013) is certainly the most refined method of analysis, which at the same time requires particularly a high computational cost and provides results that may be unfamiliar to structural engineers. Non-linear dynamic analysis can be a quick tool for analyses of simple structural models, which is actually the case of macro-elements (Pantò et al., 2016;Cannizzaro et al., 2018) or fiber-beam element approaches (De Santis and de Felice, 2014).
Another method of analysis used for seismic design of masonry structures is the non-linear static analysis. Indeed, non-linear static analyses can be developed with different analysis strategies: 3D finite element models Scozzese et al., 2019), macro-element (Caliò et al., 2012;Cannizzaro et al., 2018), fiber-beam models, FEM/DEM models (Milani and Lourenço, 2012;Baraldi et al., 2018) etc. Other authors stressed out the necessity of experimental campaigns to calibrate numerical analyses and the relevance of the choice of the control node before the non-linear static analysis of the bridge (Pelà et al., 2005;Rota et al., 2010;Zampieri et al., 2015b). In addition, a nonlinear static analysis procedure called "non-linear kinematic analysis" is also possible if the bridge is discretized into rigid blocks (Galassi et al., 1839(Galassi et al., , 2018aCavalagli et al., 2017;Galassi and Tempesta, 2019a,b) and the capacity curve of the collapse mechanism is obtained.
In this paper, the authors reported the critical review of the results from the non-linear static analysis using two modeling strategies: FEM with 3D brick elements on the one hand and then with fiber-beam elements. Besides, the results obtained from the non-linear kinematic analysis will be also discussed.

THEORETICAL BACKGROUND
In this section, the models used for the analysis of the structure under study are briefly described. Two different approaches are considered: the first consists of a 3D finite element model (3D-element approach) and the second is made of fiber-beam elements (fiber-beam element approach). The results obtained from the two approaches will then be compared with the results from non-linear kinematic analysis (N-L-A).

3D Element Approach
In this approach, a three-dimensional model is created using MIDAS FEA (Midas Fea v1.1., 2016). The model includes all the structural parts participating to the seismic resistance of the bridge. Prismatic finite elements with a triangular base were used for the mesh and each finite element is characterized by 18 degrees of freedom in space with an average size of 0.35 m. The numerical analysis of the structure considers the masonry as a continuous medium. "Total Strain Crack Model" (TSCM) 1 refers to the constitutive law relationship describing the mechanical behavior of the masonry. It displays an elastic behavior in traction with linear softening after the ultimate tensile strength is reached (Figure 1). It is worth noting that the ultimate tensile strength (f t ) in this study corresponds to 1/100 of the compressive strength (f c ). On the other hand, the constitutive relationship for compressive stresses is linear elastic-perfectly plastic was considered to describe the compressive behavior of the masonry (Figure 1). Finally, the filling material is modeled using a Mohr-Coulomb failure criterion.

Fiber-Beam Element Approach
An adequate number of experimental tests referred in Brencich and Gambarotta (2005) and De Felice and De Santis (2010) confirm that the plane section hypothesis of the beam element is also valid when it experiences a non-linear deformation. For these reasons, the solution considered to discretize the arch by means of straight beams (composed of fibers) is a valid trade-off between the model itself and the accuracy of the solution obtained.
The non-linear constitutive law of the beam is illustrated as follow.
In Figures 1B,C, the longitudinal strain along the x-axis of the generic section A(x) of the fiber-beam element is defined as , where χ is defined as the curvature in the z and y directions and ε is the axial strain. According to the plane section hypothesis, the strain of the section is defined as ε x, y, z = a y, z d (x) , where a y, z = −y z 1 . Therefore, considering the elastic modulus E x, y, z , the distribution of stresses σ x, y, z = E x, y, z ε x, y, z is immediately defined. From these relationships, it is possible to obtain the stiffness matrix k (x) of the section and the vector defining the resisting actions D (x) = M z (x) M y (x) N x in the case of continuous bridge structure throughout Equations (1) and (2): On the other hand, when the fiber beam model is used, the generic section of the beam element is divided into n(x) fibers. Equations (1) and (2) become Equations (3) and (4) and are calculated as follows: A i (x) is the area of the i-th fiber of the reference section x. Finally, the stiffness matrix for a single element is obtained as: where F represents the flexibility matrix, defined as follows: L is the length of the element, k(x) the section stiffness matrix and b(x) the matrix function of the force.

Non-linear Kinematic Analysis
Non-linear kinematic seismic analyses of the masonry bridge allow the design of the seismic capacity curve of the bridge collapse mechanism (taking as a control point the keystone centroid of the central arch). For the longitudinal non-linear kinematic analysis, the structure is considered as made of n rigid blocks (Figure 2A).  In order to consider the resisting effect of the backing, a macroblock between arches and top piers is considered as shown in Figure 2B.
The capacity curve of the bridge is calculated by applying the principle of virtual work (PVW) to the deformed shape of the structure: Where: P i is the resultant of the forces directly applied in the i-th block (e.g., the weight of the block).
Q j is the resultant of the forces not directly applied in the block but transmitted by the structure and which generate a horizontal seismic force.
θ is the generalized displacement (e.g., a rotation) taken as a reference; θ is the rotation of the arch segment H1H2 ( Figure 3A). d x,i and d x,j are the horizontal virtual displacements of the point of application of the i-inertia force of the arch and the j-inertia force of the infill ( Figure 2B). d y,i and d y,j are the vertical virtual displacements of the point of application of the force P i and the point of application of the force Q j (Figure 2B).
Finite compressive strength of the masonry is considered moving the hinge from the edge of the arch to the neutral axis depth u c in the ultimate condition as represented in Figure 2C.
With the application of PVW for finite incremental values of the rotation θ , the collapse multiplier α(θ ) and the associated displacement d i (θ ) for each i-block can be obtained. The procedure ends with the derivation of a displacement for the corresponding collapse multiplier α(θ ) = 0. In this case, the equation of PVW becomes.
The first step prior to the transversal kinematic analysis is to subdivide the structure in macro-elements as shown in Figure 3. If n is considered as the number of bridge spans, then the structure is subdivided into n+1 macro-blocks. The two abutments as well as half of the two outer spans, are considered as two macro-blocks whereas the other n-1 blocks are centered bridge spans. Each arch crown is identified as the starting and ending of the inner blocks of the bridge. This assessment methodology is efficient when dealing with the kinematic approach. Further information can be found in Zampieri et al. (2015a).

Non-linear Static Assessment
Once the curves describing the structural behavior (capacity curve) have been obtained following a non-linear static analysis, the structure is assessed. To do this, the force-displacement curve of the MDOF system is converted into a bilinear curve of an equivalent SDOF system.
Where, F * e d * are respectively, the force and displacement of the equivalent SDOF F b e d c are respectively, the force and displacement of the actual system MDOF Γ is the modal participation factor, defined as: Where, τ is the dragging vector relative to the direction of the investigated bridge ϕ is the vector relative to the fundamental vibration mode of the real system normalized by setting d c = 1 M is the matrix of the mass of the real system The assessment of the performance of the structure is carried out considering the following formulas.
The period of the bilinear equivalent system SDOF (T * ) is given by: FIGURE 4 | Geometrical properties of the case study.
Frontiers in Built Environment | www.frontiersin.org Where, m * = Φ · M · τ is the mass of the equivalent SDOF system k * is the stiffness of the elastic part of the bilinear curve If T * ≥ T C virgule the demand in displacement for the inelastic system (S Dp (T * )) is assumed to be the same as that of an elastic system of the same period (S De (T * )): If T * < T C virgule the demand in displacement for the inelastic system is greater than that of an elastic system of the same period, obtained as follows: Where q * = S e (T * ) · m * /F * y is the ratio between the elastic response force and the yield strength of the equivalent system. If q * ≤ 1 si ha che S Dp T * = S De T * .

Non-linear Kinematic Assessment
In the non-linear kinematic seismic analysis, the capacity curve of the structure must be transformed into equivalent SDOF system. The spectral acceleration a * (θ ) and the spectral displacement d * (θ ) of equivalent SDOF system are computed as follows.
d x,k is the finite horizontal displacement of the generic point P of the system (assumed as representative to plot the pushover curve) and δ x,k is its virtual horizontal displacement; M * is the mass participation of the structure to create the mechanism computed as: With non-linear kinematic analysis, the seismic assessment is expressed in terms of displacement and the verification is satisfied when: Where d * , u d * , s T * S and S De (T * S ) represent respectively, the ultimate displacement, elastic displacement, secant period and spectrum elastic displacement for the period T * S . Frontiers in Built Environment | www.frontiersin.org

CASE STUDY: MULTI-SPAN BRIDGE
The bridge typology investigated in the case study is quite similar to the Rio Claro Bridge before its collapse in 2010. The bridge is made of seven masonry arch spans and the geometrical characteristics have been chosen by the authors to reflect current existing arch bridges. Due to the lack of resources, no experimental campaign has been done in this work. It is worth noting that the mechanical characteristics of the masonry brick elements found in Table 1 derived partly from the experimental campaign discussed in Reccia et al. (2014) while the other parameters were suggested in Barbieri (2019). As shown in Figure 4, the bridge consists of semi-circular arches with a span of 12.5 m and a rise of 6.1 m. The vaults are 1.0 m thick. The piers, also in masonry, have a variable height of 18.1, 15.0, and 9.6 m. The masonry abutments are 7.5 m high. A 4.5 m thick of the backing is considered starting from the arch level. The filling material consists of loose material, for a height of 0.8 m from the extrados in keystone to the arches and, above this, a layer of ballast equal to 0.8 m; both layers are confined by a 0.8 m thicker spandrel walls. The transversal dimension of the bridge is 5.0 m.  To match as close as possible the deterioration observed in existing bridge, vertical cracks between the arch and the pier as illustrated in Figure 4 have been considered in all the arches during the numerical modeling.
As already described in the previous paragraph, two types of finite element models are considered. The first concerns a 3D finite element model presented in Figure 5, in which the resisting elements of the bridge are discretized by means of tetrahedral elements. The mechanical parameters of these elements are listed in Table 1.
The next model considered is illustrated in Figure 6 and concerns a fiber-beam element model in which the resisting elements of the bridge such as piers and arches are discretized by means of beam elements. On the other hand, the backing and filling materials are discretized by means of cut-off bar elements without mass, to which resistance characteristics have been assigned for compression only, to work as connecting rod elements. Figure 6B shows the three fiber-beam element models analyzed, in which different elements of the structure participating in the seismic resistance of the bridge are considered from time to time. In particular: in the Fiberbeam element model 1 virgule only piers and arches are considered to be resistant; in fiber-beam element model 2 virgule the addition of the backing and in fiber-beam element model 3 virgule the additional filling is considered to be resistant.
The weights of the remaining elements, which are not considered as resisting elements from time to time, are applied as forces to the arches.

Modal Analysis Results
The two FEM modeling strategies of the bridge under study are now compared. Firstly, the comparison between the vibration modes obtained through a modal analysis is done. Table 2 contains the results of the period (T) and From the modal analysis although the deformed modes of vibration are substantially similar, the periods differ by 44.5% in the longitudinal direction (X) and by 21.0% in reference to the transversal direction (Y). Less evidence is observed for  Frontiers in Built Environment | www.frontiersin.org mass participation. Indeed, the difference is up to 5.1% in the longitudinal direction (X) and 3.5% points for the transversal direction (Y).

Non-linear Static Analysis Results
The next comparison between the two types of finite element modeling concerns a non-linear static analysis. In this regard, the analysis was carried out by applying a force proportional to the masses of the structure in both longitudinal and transversal directions of the bridge (which corresponds respectively, to the X-axis and the Y-axis of the models). As far as the fiberbeam element model is concerned, we consider other cases: The contribution from the filling material is neglected (namely "Fiber-beam element model 2") and the backfilling and filling are simultaneously considered (called "Fiber-beam element model 3"). The curves resulting from the above analyses are illustrated in Figure 8, which shows the normalized force against the displacement of the control point. It is worth noticing that the control point in this research is taken at the arch intrados of the bridge mid-span.
Considering Figure 8, it is possible to draw the following conclusions. For fiber-beam element models, as the number of elements considered to participate in the seismic resistance of the bridge increases, the value of ultimate force increases. In fact, the resistance changes from 11057.6 kN for the Fiber-beam element model 1 to 24769.1 kN for the Fiber-beam element model 3. On the other hand, the value of ultimate displacement decreases, from 0.05 m for the Fiber-beam element model 1 to 0.025 m for the Fiber-beam element model 3.
Similar observations can be made by comparing the 3D FEM model and the corresponding Fiber-beam element model.
With regard to the non-linear static analysis in the transversal direction, the behavior between the different modeling strategies is quite similar. Indeed, the normalized pushover curves for the three fiber-beam models overlay each other ( Figure 8B). Furthermore, the ultimate resistance varies from 16807.6 kN for the 3D FEM model to 17692.2 kN for the fiber-beam element models. Moreover, the ultimate displacement varies from 0.115 m for the 3D FEM model to 0.103 m for the fiber-beam element models. It should be pointed out that the ultimate resistance and ultimate displacement remain the same for the three fiber-beam element models.
For each step of the pushover curve from the relevant analysis, the verification of the shear-slip should be done throughout the following expression: Where l ′ is the length of the compressed part of the section t is the thickness of the section f vk0 is the characteristic shear strength in the absence of vertical loads σ n is the average normal stress due to the vertical loads acting on the design elements Γ M is the safety coefficient, considered equal to 1.0 in the case of non-linear static analysis.
Considering the 3D-element modeling approach, the structure displays damages at the base of the piers. In the transversal direction, the damage is less significant as shown in Figure 9A. In particular, the observed damage is relevant for pier 6, arches 1, 2, and 7 and for the right abutment ( Figure 9B).
From the analysis in the transversal direction through the 3D-element approach, it is possible to appreciate the collapse of the spandrel that occurs before reaching the capacity in that direction. For both spandrels (namely spandrel A and spandrel B) (see Figure 10), this collapse consists of the formation of cracks in the areas mainly subjected to the traction of the element. These cracks are more pronounced as the applied load increases (force proportional to the masses of the structure).
The damage involving the spandrel, as illustrated in Figure 10, consists of the formation of cracks mainly at arches 1 and 7 for the spandrel A and arches 3, 4, and 5, for the spandrel B.

Non-linear Kinematic Analysis Results
From the non-linear kinematic analysis, it was possible to determine the collapse mechanism of the bridge in longitudinal and transversal direction (Figure 11). It is worthy to record that the analyses performed and presented in this sub-section refer to the use of the rigid block model already introduced. On the first hand, the longitudinal collapse mechanism (Figure 11A) is characterized by the formation of hinges at the base of the piers and on the arches. The collapse mechanism is defined through an iterative process that provides a minimum kinematically admissible collapse multiplier. a) Collapse mechanism in the longitudinal direction b) Collapse mechanism in the transversal direction.   On the second hand, the transversal collapse mechanism ( Figure 11B) is defined without using an optimization process based on the previous authors' work and reported in Zampieri et al. (2015a) The capacity curves obtained from non-linear kinematic analysis show that a * 0 is equal to 0.162 g and 0.157 g, respectively, for longitudinal and transversal analysis while d * 0 is equal to 0.182 and 2.155 m, respectively.

SEISMIC VULNERABILITY ASSESSMENT
From the force-displacement curves that describe the behavior of the structure under the effect of seismic loads, it was also possible to investigate the seismic vulnerability of the bridges under study. This assessment is carried out by comparing the ultimate capacity in displacement of the structure (obtained by the pushover analysis) with the structural demand.
The Performance Indicator (PI), which best describes the relationship between the structural capacity and the demand are defined as follow respectively, for the non-linear static and the non-linear kinematic analysis.
The seismic vulnerability can also be investigated throughout the response spectrum from Eurocode 8 considering the following parameters: Soil category: B; Design ground acceleration: a g = 0.2 g; Soil factor: S S = 1.2; Topographic coefficient: S T = 1.0; Maximum spectrum amplification factor F 0 = 2.5; Damping coefficient (5.0 %).
The analysis of the 3D-element approach shows the following results (see Figures 12A,B). In particular, the most severe conditions are obtained with reference to the longitudinal direction (X-axis in the model) with PI values equal to 0.62 and a final displacement of 0.011 m, compared to a PI of 0.93 found for the transversal direction (Y-axis in the model) and a final displacement of 0.034 m.
The results from the analysis of the fiber-beam element model are shown in Figures 12C,D. In particular, the most severe conditions are obtained with reference to the longitudinal direction (X-axis in the model) with PI values of 0.52 and a final displacement of 0.013 m, compared to a PI of 0.85 for the transversal direction (Y-axis in the model) and a final displacement of 0.030 m.
From the non-linear kinematic analysis, the results obtained are presented in Figure 13. In particular, in this case, the verifications are both satisfied. From these, it can be concluded non-linear kinematic provides results that are not comparable with the ones obtained with the other methods illustrated. This difference probably depends on the flexural behavior of the slender piers that cannot be considered with a rigid blocks schematization of the bridge. However, the displacement and the period d * u and T * S , depending on a * 0 e d * 0 may not be representative for masonry bridge. These values have been calibrated indeed for masonry buildings. Therefore, the methodology could produce results which are comparable to those obtained from non-linear static analysis provided that the values of d * u and T * S are calibrated specifically for bridges in the scope of future research.
By comparing the PI values obtained from the analyses, it is possible to notice how the fiber-beam element model gives lower PI values than those of the 3D finite element model. These results show that the fiber-beam element model, while satisfactorily approximating the behavior of the structure (also in comparison to the 3D finite element model) is more favorable to safety by lowering PI values in respect to the 3D finite element model.

CONCLUSION
In order to study the seismic behavior of the multi-span masonry bridge (chosen as a case study) three different approaches based on non-linear static analysis were investigated. Two of them based on different F.E modeling strategies (3D-element approach and fiber-beam element approach) and the third is a rigid block analysis. The 3D-element approach was able to define all the collapse mechanism of the bridges: Global Longitudinal collapse mechanism, Global Transversal collapse mechanism and the Local Transversal mechanism of the spandrel wall.
• The capacity curve was found from all the different approaches and the seismic vulnerability expressed in term of capacity displacement/demand displacement was investigated. Finally, a comparison of the different results was done leading to the following conclusions: • The fiber-beam element approach, even if it is a simplified approach and limited to describe all the possible seismic collapse mechanisms of the bridges carries out results that are comparable with that of the 3D-element approach. However, this approach underestimates the horizontal stiffness of the bridges; • Kinematic analysis can represent in a clear way the collapse mechanisms of the bridges and its maximum resisting acceleration (a * 0 ); • Non-linear kinematic seismic assessment procedure displays results completely different with respect to the other approaches. This result depends on T * S and d * u parameters and are valid for masonry buildings.
• Future research and development could involve the calibration of T * S and d * u in other to perform reliable non-linear kinematic analysis of masonry bridges.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
PZ coordinated the work, performed some numerical analysis, and edited the document. SP performed most of the numerical analysis and edited the document. CT performed some numerical analysis and edited the document. CP supervised the work.