Safety Evaluation of a Water-Immersed Bridge Against Multiple Hazards via Machine Learning

A safety analysis process for an immersed bridge is proposed, in which material nonlinear properties, added mass, soil spring, plastic hinge length and property, nonlinear pushover analysis, and time history analysis are included. Outcomes of interest in a bridge analysis, such as displacement ductility, is not an easy task if the authentic model is adopted. The least-squares support-vector machine (LSSVM) is therefore proposed for replacing the analysis model. Water-immersed bridges under varied water depths, scouring depth, and seismic intensity are investigated. The feasibility of using machine learning as a surrogate model is discussed. Results indicate that the trends of fragility curves derived from LSSVM are very similar to those of the authentic model. Therefore, LSSVM is a suitable tool to reduce the computational burden. Analyses also show that stream velocity and pier length are two important factors for the safety of an immersed bridge, and the immersed effect should not be ignored if pier length is long.


Introduction
In Taiwan, interweaving rivers and creeks are spread over the entire island. To achieve smooth southbound and northbound traffic flow, river-crossing bridges are prerequisite infrastructure facilities. Because these infrastructure facilities are usually designed with a very long lifespan, it is also required to consider the effects from all aspects when designing these structures. As Taiwan is situated in a subtropical monsoon area, river-related hazards should not be treated lightly because of the abundant rainfall and the striking of typhoons during June-August each year. In the meantime, Taiwan is also located at the interfacing area between the Eurasian and Philippine Sea Plates. Therefore, it sits on the seismic zone. When these plates drift and squeeze with each other, Taiwan will be shaken by earthquakes. As such, impacts to structures as caused by earthquakes should not be ignored, and this is the key point that will be dealt with by each engineer.
Many researchers have focused on building the relationship between earthquake intensities and bridge/building damages. For example, Bazos et al. [1] suggested that fragility curves are a practical tool. Based on the observation at bridge sites, Hsu and Fu [2] found bridges were damaged in different ways in the Chi-Chi earthquake such as unseating span failure, abutment failure, joint failure, substructure damage, footing settlement, and so on. Elnashai et al. [3] developed site-specific ground motions to analyze several typical failure bridge modes and concluded that displacement is a good indicator to measure bridge safety. In addition to the seismic hazard, 44 bridges were damaged from Hurricane Katrina [4]. Andric' and Lu [5] reported that flood-induced damage was the primary damage for bridges in the United States. Similarly, Taiwan faces the same challenge. Liao et al. [6] proposed a bridge safety evaluation process against seismic and flood hazards, in which a scour prediction equation was constructed to measure the flood risk, and the probabilistic seismic-hazard

Bridge Modeling
In this research, nonlinear behaviors of steel bars were established according to the theory proposed by Priestley et al. [12]. As the steel bars selected were 36 mm in diameter (#11), its limit strain value ε su was set at 0.09, and the stress-strain hardening initial value ε sh was set at 0.0115. The corresponding stress-strain curve was drafted as shown in Figure 3. In this research, the bridge pier was configured as a tapered cylindrical column where the top and bottom diameters were 1.8 and 2.4 m, respectively. Because different steel bars were used, the mechanical behaviors of the sectional profile will vary. For this reason, the variable section was divided into 5 sectional profiles in order to calculate the nonlinear mechanics behaviors. Accordingly, 4 tapered cylindrical columns were generated in SAP2000. For example, the topmost portion used the sectional properties of the 1.8 and 1.95 m in diameter, as shown in Figure 4. The interpolated sectional diameters of different height levels were 1.8, 1.95, 2.1, 2.25, and 2.4 m.
The strength of the concrete will be affected by the confining force of stirrups. For this reason, they will be divided into two portions for calculation. One of which is the protection layer, which is the unconfined concrete; and the other is the core concrete, which is a confined layer. For the calculation of confining concrete, the method proposed by Mander et al. [13] was adopted here. Because the steel bars used for this research were 36 mm (#11) in diameter, the limit strains (ε su ) of the steel bar and concrete were 0.09 and 0.004, respectively. Note that the limit strain for concrete is 0.003 if American Concrete Institute (ACI) is adopted. After being calculated, the stress-strain curves of the desired sectional profile were obtained ( Figures 5 and 6) for confined and unconfined concretes.

Bridge Modeling
In this research, nonlinear behaviors of steel bars were established according to the theory proposed by Priestley et al. [12]. As the steel bars selected were 36 mm in diameter (#11), its limit strain value su ε was set at 0.09, and the stress-strain hardening initial value sh ε was set at 0.0115.
The corresponding stress-strain curve was drafted as shown in Figure 3. In this research, the bridge pier was configured as a tapered cylindrical column where the top and bottom diameters were 1.8 and 2.4 m, respectively. Because different steel bars were used, the mechanical behaviors of the sectional profile will vary. For this reason, the variable section was divided into 5 sectional profiles in order to calculate the nonlinear mechanics behaviors. Accordingly, 4 tapered cylindrical columns were generated in SAP2000. For example, the topmost portion used the sectional properties of the 1.8 and 1.95 m in diameter, as shown in Figure 4. The interpolated sectional diameters of different height levels were 1.8, 1.95, 2.1, 2.25, and 2.4 m.
The strength of the concrete will be affected by the confining force of stirrups. For this reason, they will be divided into two portions for calculation. One of which is the protection layer, which is the unconfined concrete; and the other is the core concrete, which is a confined layer. For the calculation of confining concrete, the method proposed by Mander et al. [13] was adopted here.
Because the steel bars used for this research were 36 mm (#11) in diameter, the limit strains ( su ε ) of the steel bar and concrete were 0.09 and 0.004, respectively. Note that the limit strain for concrete is 0.003 if American Concrete Institute (ACI) is adopted. After being calculated, the stress-strain curves of the desired sectional profile were obtained ( Figures 5 and 6) for confined and unconfined concretes.

Bridge Modeling
In this research, nonlinear behaviors of steel bars were established according to the theory proposed by Priestley et al. [12]. As the steel bars selected were 36 mm in diameter (#11), its limit strain value su ε was set at 0.09, and the stress-strain hardening initial value sh ε was set at 0.0115.
The corresponding stress-strain curve was drafted as shown in Figure 3. In this research, the bridge pier was configured as a tapered cylindrical column where the top and bottom diameters were 1.8 and 2.4 m, respectively. Because different steel bars were used, the mechanical behaviors of the sectional profile will vary. For this reason, the variable section was divided into 5 sectional profiles in order to calculate the nonlinear mechanics behaviors. Accordingly, 4 tapered cylindrical columns were generated in SAP2000. For example, the topmost portion used the sectional properties of the 1.8 and 1.95 m in diameter, as shown in Figure 4. The interpolated sectional diameters of different height levels were 1.8, 1.95, 2.1, 2.25, and 2.4 m.
The strength of the concrete will be affected by the confining force of stirrups. For this reason, they will be divided into two portions for calculation. One of which is the protection layer, which is the unconfined concrete; and the other is the core concrete, which is a confined layer. For the calculation of confining concrete, the method proposed by Mander et al. [13] was adopted here.
Because the steel bars used for this research were 36 mm (#11) in diameter, the limit strains ( su ε ) of the steel bar and concrete were 0.09 and 0.004, respectively. Note that the limit strain for concrete is 0.003 if American Concrete Institute (ACI) is adopted. After being calculated, the stress-strain curves of the desired sectional profile were obtained ( Figures 5 and 6) for confined and unconfined concretes.    Depending on the damage mode, the plastic hinge used for the finite element model had the following properties: flexural damage (corresponding to the bending moment of plastic hinge M2 or M3 in SAP2000), shear damage (corresponding to shear plastic hinge V2 or V3 in SAP2000), and flexural-shear damage (corresponding to the axial force interacting with the bending moment of plastic hinge P-M2-M3 in SAP2000). In SAP2000, the aforementioned three types of plastic hinges can be simulated in three different ways-link, sectional hinge, and fiber-which will vary in accuracy and analysis time. Aviram et al. [14] suggested that the fiber hinge possesses the best accuracy; therefore, it was selected in this article. The fiber hinge resembles the bending moment of the plastic hinge featuring an axial force interaction (P-M2-M3). During analysis, the profile of a target section will affect the analysis time according to the mesh size. For the sake of balancing accuracy and time, the trial-and-error method was used for deciding the mesh size. The finalized meshes for pier and caisson are described in Figures 7 and 8.   Depending on the damage mode, the plastic hinge used for the finite element model had the following properties: flexural damage (corresponding to the bending moment of plastic hinge M2 or M3 in SAP2000), shear damage (corresponding to shear plastic hinge V2 or V3 in SAP2000), and flexural-shear damage (corresponding to the axial force interacting with the bending moment of plastic hinge P-M2-M3 in SAP2000). In SAP2000, the aforementioned three types of plastic hinges can be simulated in three different ways-link, sectional hinge, and fiber-which will vary in accuracy and analysis time. Aviram et al. [14] suggested that the fiber hinge possesses the best accuracy; therefore, it was selected in this article. The fiber hinge resembles the bending moment of the plastic hinge featuring an axial force interaction (P-M2-M3). During analysis, the profile of a target section will affect the analysis time according to the mesh size. For the sake of balancing accuracy and time, the trial-and-error method was used for deciding the mesh size. The finalized meshes for pier and caisson are described in Figures 7 and 8. Depending on the damage mode, the plastic hinge used for the finite element model had the following properties: flexural damage (corresponding to the bending moment of plastic hinge M 2 or M 3 in SAP2000), shear damage (corresponding to shear plastic hinge V 2 or V 3 in SAP2000), and flexural-shear damage (corresponding to the axial force interacting with the bending moment of plastic hinge P-M 2-M 3 in SAP2000). In SAP2000, the aforementioned three types of plastic hinges can be simulated in three different ways-link, sectional hinge, and fiber-which will vary in accuracy and analysis time. Aviram et al. [14] suggested that the fiber hinge possesses the best accuracy; therefore, it was selected in this article. The fiber hinge resembles the bending moment of the plastic hinge featuring an axial force interaction (P-M 2-M 3 ). During analysis, the profile of a target section will affect the analysis time according to the mesh size. For the sake of balancing accuracy and time, the trial-and-error method was used for deciding the mesh size. The finalized meshes for pier and caisson are described in   This article used the formula defined by Paulay and Priestley [15] to calculate the length (Lp) of plastic hinge as below: where h means pier length (mm), bl d means diameter of the longitudinal steel bar (mm), and y σ means steel bar yielding stress (MPa). Soil resistance around the foundation was simulated by a spring element, in which the front and bottom sides of the foundation were separately considered. During simulation, the soil spring was assumed as a bilinear spring, in which the maximum value did not exceed the passive soil pressure. The soil springs on the front side were deployed at 1 m intervals, and the soil spring was removed if scour occurred.

Effect of Immersing in Water
As the natural frequency of the river-crossing bridge will vary according to the water level, the influence of immersing the structure in water should be discussed. The immersing effect is considered using the Morison equation [16].
The Morison equation is a semiempirical equation for calculating the force acted on a body under oscillatory flow. There are two parts in the Morison equation, an inertia force that depends on the acceleration of the surrounding flow and a drag force that depends on the velocity of the surrounding flow. Two empirical hydrodynamic coefficients are needed when the Morison equation is adopted: an inertia coefficient (Cm) and a drag coefficient (Cd). These two coefficients depend on the Keulegan-Carpenter number, Reynolds number, and surface roughness [17]. Yang and Li [18] proposed the added mass ratio method (AMRM) according to the Morison equation, by which a rational method was inducted for calculating the added mass and was verified with the finite element method. In this article, AMRM was adopted, and the added mass ( cir m Δ ) calculation equation is described as below:  This article used the formula defined by Paulay and Priestley [15] to calculate the length (Lp) of plastic hinge as below: where h means pier length (mm), bl d means diameter of the longitudinal steel bar (mm), and y σ means steel bar yielding stress (MPa). Soil resistance around the foundation was simulated by a spring element, in which the front and bottom sides of the foundation were separately considered. During simulation, the soil spring was assumed as a bilinear spring, in which the maximum value did not exceed the passive soil pressure. The soil springs on the front side were deployed at 1 m intervals, and the soil spring was removed if scour occurred.

Effect of Immersing in Water
As the natural frequency of the river-crossing bridge will vary according to the water level, the influence of immersing the structure in water should be discussed. The immersing effect is considered using the Morison equation [16].
The Morison equation is a semiempirical equation for calculating the force acted on a body under oscillatory flow. There are two parts in the Morison equation, an inertia force that depends on the acceleration of the surrounding flow and a drag force that depends on the velocity of the surrounding flow. Two empirical hydrodynamic coefficients are needed when the Morison equation is adopted: an inertia coefficient (Cm) and a drag coefficient (Cd). These two coefficients depend on the Keulegan-Carpenter number, Reynolds number, and surface roughness [17]. Yang and Li [18] proposed the added mass ratio method (AMRM) according to the Morison equation, by which a rational method was inducted for calculating the added mass and was verified with the finite element method. In this article, AMRM was adopted, and the added mass ( cir m Δ ) calculation equation is described as below: This article used the formula defined by Paulay and Priestley [15] to calculate the length (L p ) of plastic hinge as below: where h means pier length (mm), d bl means diameter of the longitudinal steel bar (mm), and σ y means steel bar yielding stress (MPa). Soil resistance around the foundation was simulated by a spring element, in which the front and bottom sides of the foundation were separately considered. During simulation, the soil spring was assumed as a bilinear spring, in which the maximum value did not exceed the passive soil pressure. The soil springs on the front side were deployed at 1 m intervals, and the soil spring was removed if scour occurred.

Effect of Immersing in Water
As the natural frequency of the river-crossing bridge will vary according to the water level, the influence of immersing the structure in water should be discussed. The immersing effect is considered using the Morison equation [16].
The Morison equation is a semiempirical equation for calculating the force acted on a body under oscillatory flow. There are two parts in the Morison equation, an inertia force that depends on the acceleration of the surrounding flow and a drag force that depends on the velocity of the surrounding flow. Two empirical hydrodynamic coefficients are needed when the Morison equation is adopted: an inertia coefficient (C m ) and a drag coefficient (C d ). These two coefficients depend on the Keulegan-Carpenter number, Reynolds number, and surface roughness [17]. Yang and Li [18] proposed the added mass ratio method (AMRM) according to the Morison equation, by which a rational method was inducted for calculating the added mass and was verified with the finite element method. In this article, AMRM was adopted, and the added mass (∆m cir ) calculation equation is described as below: where ρ con represents the density of concrete, and the added mass ratio p cir (H, D) of each cylinder section is calculated as below: in which H means the cylinder height, and D means the diameter of the cylinder section. Because the bridge piers of this study were tapered, their added mass was calculated from the pier bottom for every 1 m interval, and then the calculated added mass will be applied to each point, as shown in Figure 9.
Appl. Sci. 2019, 9, where con ρ represents the density of concrete, and the added mass ratio in which H means the cylinder height, and D means the diameter of the cylinder section. Because the bridge piers of this study were tapered, their added mass was calculated from the pier bottom for every 1 m interval, and then the calculated added mass will be applied to each point, as shown in Figure 9.

Simulation of Water Depth, Stream Velocity, and Scouring Depth
When calculating bridge scouring depth, two parameters must be obtained: water depth and stream velocity. Because the water depth and the stream velocity should be observed over time, normally they are calculated through simulation. Presuming that the mean value of water depth and stream velocity in different return years can be regarded as a log-normal distribution with three parameters, the water depth and stream velocity of the same return year can also be regarded as a log-normal distribution. Based on the above assumptions and the information revealed in the design drawing (i.e., the discharge rate for return periods of 50 and 100 years), the mean value of water depth 200 y and stream velocity 200 V for a 200-year return period were obtained. Furthermore, the coefficients of variation of water depth and stream velocity were presumed identical under different return periods. In this way, we may acquire the probability distributions of water depth and stream velocity for 50, 100, and 200 year return periods, as per the detailed calculation method. Provided below is the hydrological frequency equation:

Simulation of Water Depth, Stream Velocity, and Scouring Depth
When calculating bridge scouring depth, two parameters must be obtained: water depth and stream velocity. Because the water depth and the stream velocity should be observed over time, normally they are calculated through simulation. Presuming that the mean value of water depth and stream velocity in different return years can be regarded as a log-normal distribution with three parameters, the water depth and stream velocity of the same return year can also be regarded as a log-normal distribution. Based on the above assumptions and the information revealed in the design drawing (i.e., the discharge rate for return periods of 50 and 100 years), the mean value of water depth y 200 and stream velocity V 200 for a 200-year return period were obtained. Furthermore, the coefficients of variation of water depth and stream velocity were presumed identical under different return periods. In this way, we may acquire the probability distributions of water depth and stream velocity for 50, 100, and 200 year return periods, as per the detailed calculation method. Provided below is the hydrological frequency equation:x where x and S represent mean value and standard deviation of the sample, respectively. K T represents frequency factor and is calculated below; where C s means the skewness. Parameter t in Equation (4) can be acquired by the formula below: where P is the inverse of the return period, C 0 = 2.515517, To calculate scour depth, especially to obtain the uncertainty of scour, a simulation of water level and stream velocity was necessary. Liao et al. [19] indicated that the correlation coefficient between water depth and stream velocity was about 0.92. Based on this, a Copula function was used in this research for generating random samples for water depth and stream velocity. After acquiring the water depths and stream velocities for 50, 100, and 200 year return periods, they were revolved in the formula proposed by Melville and Coleman [20], respectively, in order to acquire the corresponding distribution of scouring depth, as shown below.
In the formula, d s represents the scour depth (m), K s means foundation shape factor, K t means time factor, K G means approach channel geometry factor, K θ means foundation alignment factor, K yb means foundation size factor, K I means flow intensity factor, and K d means sediment size factor. Among these factors, K yb poses the most significant impact. Detailed calculations of K yb in this article is provided below. Based on the classification of caisson width and elevation of river bed, a pier-equivalent width was derived for calculating the bridge scour [21]. Liao [22] utilized an optimization approach to update the parameter weights for the pier-equivalent width formula proposed by Melville and Raudkivi [21], as displayed below.
where b e represents the equivalent width of bridge, b means the width of pier, b pc means the width of caisson, y means the depth from water face to riverbed, and Y means the elevation from caisson top to riverbed. The foundation size factor K yb will vary according to the ratio between bridge equivalent width b e and water depth and its calculation formula is as below: Based on Liao's updated formula, the probability of exceedance and occurrence for the scour distribution for the 100 year return period are shown in Figures 10 and 11, respectively.
Based on Liao's updated formula, the probability of exceedance and occurrence for the scour distribution for the 100 year return period are shown in Figures 10 and 11, respectively.

Seismic Hazard Analysis
Based on the seismic design code in Taiwan, horizontal spectrum acceleration coefficients SDS and SDm were obtained for the investigated bridge site. Based on the short period peak ground motion (PGA) value, the EPA (effective peak ground acceleration) of the design earthquake (about 475 years return period) and the maximum consideration earthquake (about 2500 years return period) can be calculated as below.  Based on Liao's updated formula, the probability of exceedance and occurrence for the scour distribution for the 100 year return period are shown in Figures 10 and 11, respectively.

Seismic Hazard Analysis
Based on the seismic design code in Taiwan, horizontal spectrum acceleration coefficients SDS and SDm were obtained for the investigated bridge site. Based on the short period peak ground motion (PGA) value, the EPA (effective peak ground acceleration) of the design earthquake (about 475 years return period) and the maximum consideration earthquake (about 2500 years return period) can be calculated as below.  Figure 11. Scouring depth probability density function for the 100 year return period.

Seismic Hazard Analysis
Based on the seismic design code in Taiwan, horizontal spectrum acceleration coefficients S DS and S Dm were obtained for the investigated bridge site. Based on the short period peak ground motion (PGA) value, the EPA (effective peak ground acceleration) of the design earthquake (about 475 years return period) and the maximum consideration earthquake (about 2500 years return period) can be calculated as below.
This study presumed that the bridge lifespan was set at 50 years. On this basis, the exceedance probability (EP) of an earthquake, in magnitude greater than or equal to 475 years and with 2500 years of reoccurrence period that has occurred at least once in such duration, shall be calculated as below: Presuming that the probability of magnitude greater than or equal to M, which is indicated by Equations (14) and (15), during the lifespan follows a log-normal distribution, then the exceedance probability distribution and the occurrence probability distribution of the bridge during the lifespan can be acquired as indicated in Figures 12 and 13.
Presuming that the probability of magnitude greater than or equal to M, which is indicated by Equations (14) and (15), during the lifespan follows a log-normal distribution, then the exceedance probability distribution and the occurrence probability distribution of the bridge during the lifespan can be acquired as indicated in Figures 12 and 13.

Fragility Curve Analysis
A fragility curve is used to indicate the failure probability of the structure under different levels of hazard, and the fragility curve is mainly affected by the demand and the capacity of the bridge. The demand means the displacement ductility (RD) analyzed by simulating a scoured bridge under Presuming that the probability of magnitude greater than or equal to M, which is indicated by Equations (14) and (15), during the lifespan follows a log-normal distribution, then the exceedance probability distribution and the occurrence probability distribution of the bridge during the lifespan can be acquired as indicated in Figures 12 and 13.

Fragility Curve Analysis
A fragility curve is used to indicate the failure probability of the structure under different levels of hazard, and the fragility curve is mainly affected by the demand and the capacity of the bridge. The demand means the displacement ductility (RD) analyzed by simulating a scoured bridge under

Fragility Curve Analysis
A fragility curve is used to indicate the failure probability of the structure under different levels of hazard, and the fragility curve is mainly affected by the demand and the capacity of the bridge. The demand means the displacement ductility (R D ) analyzed by simulating a scoured bridge under specific earthquake magnitude. As for the capacity (R C ), it is calculated according to the definition proposed by Federal Emergency Management Agency (FEMA, HAZUS-MH/MR3, [23]), where 1 < R c < 2 means minor damage, 2 < R c < 4 means medium damage, 4 < R c < 7 means major damage, and R c > 7 can be considered that the structure completely collapsed. When the demand of displacement ductility is bigger than the capacity, it can be assumed that the bridge is damaged. Here, the lower bounds of each performance level (e.g., 1, 2, 4, and 7) were used as the threshold values. The displacement ductility demand (R D ) of the bridge is the ratio between the ultimate displacement (∆ u ) and yielding displacement (∆ y ), that is, R D = ∆ u ∆ y . Assuming the pier of the bridge serves as the primary damage mechanism, in this case, the aforementioned term of "displacement" refers to the relative displacement between pier top and pier bottom. In view that it would be impossible to acquire an accurate ultimate displacement of the bridge by a nonlinear static analysis method, such as a pushover analysis, the ultimate displacement was therefore calculated by a nonlinear time history analysis, where the maximum value of the bridge displacement duration excitation is the ultimate displacement.
The yielding displacement (∆ y ) is another key in determining the displacement ductility demand (R D ), where ∆ y is acquired from a nonlinear pushover analysis. Because the pushover curve is irregular in nature, the method of acquiring ∆ y has always been a subject discussed by all sectors. In this study, the equal energy method proposed by Sung et al. [24] was adopted to calculate the yielding displacement, and it will be explained in the following three steps.
First, acquire area (A 1 ) under the pushover curve. Then, draft the tangent lines on the origin and the ending points of the curve, where the corresponding slopes are considered as the initial stiffness (K 1 ) and yielding stiffness (K 2 ). Accordingly, the intersection point (P 1 ) is also obtained, as per Figure 14a below. Next, rotate the initial stiffness (K 1 ) line from its original point clockwise until area (A 2 ), which is enclosed by the rotating and K 2 lines, becomes equal to area (A 1 ). Then, we have intersection point (P 2 ), as per Figure 14b below. Likewise, rotate the yielding stiffness (K 2 ) line from its ending point counter clockwise until area (A 3 ), which is enclosed by the rotating and K 1 lines, becomes equal to area (A 1 ). Then, we have intersection point (P 3 ), as per Figure 14c below. Finally, a line which is orthogonal to the line of P 2 P 3 and passes through the P 3 point can be drafted in which the intersection point is the yielding point (P y ), and the corresponding displacement is the yielding displacement (∆ y ), as per Figure 14d below. Indicated in Figure 15 is an example of the yielding displacement (∆ y ) calculated from a pushover curve with 0 m of water depth and scouring depth.
ultimate displacement ( u Δ ) and yielding displacement ( y Δ ), that is, D y R = Δ . Assuming the pier of the bridge serves as the primary damage mechanism, in this case, the aforementioned term of "displacement" refers to the relative displacement between pier top and pier bottom. In view that it would be impossible to acquire an accurate ultimate displacement of the bridge by a nonlinear static analysis method, such as a pushover analysis, the ultimate displacement was therefore calculated by a nonlinear time history analysis, where the maximum value of the bridge displacement duration excitation is the ultimate displacement. The yielding displacement ( y Δ ) is another key in determining the displacement ductility demand (RD), where y Δ is acquired from a nonlinear pushover analysis. Because the pushover curve is irregular in nature, the method of acquiring y Δ has always been a subject discussed by all sectors. In this study, the equal energy method proposed by Sung et al. [24] was adopted to calculate the yielding displacement, and it will be explained in the following three steps. First, acquire area (A1) under the pushover curve. Then, draft the tangent lines on the origin and the ending points of the curve, where the corresponding slopes are considered as the initial stiffness (K1) and yielding stiffness (K2). Accordingly, the intersection point (P1) is also obtained, as per Figure  14a below. Next, rotate the initial stiffness (K1) line from its original point clockwise until area (A2), which is enclosed by the rotating and K2 lines, becomes equal to area (A1). Then, we have intersection point (P2), as per Figure 14b below. Likewise, rotate the yielding stiffness (K2) line from its ending point counter clockwise until area (A3), which is enclosed by the rotating and K1 lines, becomes equal to area (A1). Then, we have intersection point (P3), as per Figure 14c    After acquiring the displacement ductility demand (RD), the failure probability can be analyzed in order to obtain the fragility curve. When analyzing the failure probability, a limit state needs to be defined, in which After acquiring the displacement ductility demand (R D ), the failure probability can be analyzed in order to obtain the fragility curve. When analyzing the failure probability, a limit state needs to be defined, in which PE = R C R D − 1 was used in this study, indicating that when PE > 0, the displacement ductility capacity (R C ) of a bridge is bigger than displacement ductility demand (R D ), and the considered bridge is under safe conditions. Otherwise, it means that the bridge will be subject to damage. As such, the failure probability, when a bridge is subject to seismic force PGA = x, can be expressed by the following formula: in which R is Assuming R C and D C are log-normal distributions, then the formula used for calculating the failure probability can be revised as: where µ ln(R C ) refers to the mean capacity and σ ln(R C ) refers to standard deviation of the capacity, and µ ln(R D ) refers to the mean demand and σ ln(R D ) refers to standard deviation of the demand. Based on the research report [25] proposed by National Center for Research on Earthquake Engineering (NCREE) of Taiwan, the standard deviations of capacity of the respective performance level (e.g., slight, moderate, extensive, and collapse) are 0.5, 0.45, 0.4, and 0.4 respectively. Φ(·) represents the standard normal cumulative distribution function. The required mean value and standard deviation of demand were obtained from the regression of power law. To do so, 3 sets of data were used: the displacement ductility demands (R D ) for the 30 year return, 475 year return, and 2500 year return periods, respectively. Each set contained 7 structural responses from time history analyses using 7 earthquake excitations. Listed below are the formulas used to calculate the PGA, mean, and standard deviation displacement ductility demand.
where a, b, c, and d represent the regression parameters. Note that at least 7 ground motions should be used if a time history analysis is used to measure the structural performance (LRFD seismic bridge design, AASTHO 2007). In addition, the "response-spectrum-compatible" time histories were used here. A response-spectrum-compatible time history refers to the response spectrum of the selected earthquakes falling between 0.2 and 1.5 T (T is the fundamental period), which may not be less than 90% of the corresponding design spectral acceleration for a damping ratio of 5%. In addition, the average value of the response spectrum within the designated period range may not be less than the average value of the corresponding design spectral accelerations. The 7 ground motions selected here were therefore converted into response-spectrum-compatible data for return periods of 30, 475, and 2500 years.

Least-Squares Support-Vector Machine (LSSVM)
The aforesaid fragility curve was drafted according to the many results of structural analyses, and it is a time-consuming process. For this reason, LSSVM was used to construct a surrogate model and to rapidly establish the fragility curve for evaluating the performance of a bridge. The concept of LSSVM is to consider the key essential parameters as input vectors for a machine, and together with given outputs, the machine can be trained to learn in order to estimate the result through the learning process. The parameters (training input) considered in this article were the PGA, water level, and scouring depth. The displacement ductility demand (R D ) obtained from the structural analysis was used as the output data during the training. The LSSVM is briefly introduced below.
SVM is a classifier that is able to solve a nonlinear problem using convex quadratic programs (QPs), as shown in Equation (20).
in which y k represents the class and (w T K (x i ) + b) indicates the classifier, w is a vector of weights that are orthogonal to the hyper-plane, c is a constant number that is greater than zero, and ξ k is the slack variable. When ξ k is greater than one, this indicates that the k-th inequality is violated. N is the number of data, and K is the kernel function, in which a Gaussian radial basis function (RBF) is one of the common kernels and was adopted here, as shown in Equation (21).
in which vector X is an input, σ represents the kernel function parameter, and X i is the support vector.
The least-square support-vector machine (LSSVM, Suykens and Vandewalle, [26]) does not attempt to solve the QP problem. LSSVM actually tries to solve a system of linear equations after altering the SVM via introducing the error variable (ε), as described in Equation (22).
in which γ is a constant number. It is seen that two modifications, equality constraints and a squared error variable, leads to solving a set of linear equations in LSSVM. When training the LSSVM, the training samples are separated into two portions. In the present case, 525 analysis data points previously obtained were split by an 8:2 ratio, randomly, in which 80% was used as the training samples and 20% was used to compare the estimate result. During training, these samples were separated into 10 groups for cross-validation in order to improve its accuracy and avoid overfitting. To implement the LSSVM in a computer, this study utilized the LSSVMlab toolbox provided by Suykens et al., which contains Matlab/C implementations for a number of LSSVM algorithms. To measure the performance of the trained LSSVM, 4 sets of indices were used: mean absolute error (MAE), mean absolute percentage error (MAPE), root-mean-square error (RMSE), and coefficient of determination (R 2 ). The hyper-parameters (α and β) of LSSVM are determined by particle swarm optimization (PSO).
The joint-failure probability is assumed to be the product of three probabilities [27]: the probability of seismic hazard, the probability of scour hazard, and bridge failure probability for a given limit state. Accordingly, the joint failure probability can be expressed as: where P(SC i ) is the probability of occurrence of a given scour depth (i), acquired from the scouring hazard analysis; P EQ j is the probability of occurrence of a certain PGA (j), acquired from seismic hazard analysis; and P(DS k ) is the exceedance probability of a bridge exceeding a certain performance level (k), acquired from time history analyses. The hyper-parameters (α and β) of LSSVM are determined by PSO. PSO is briefly introduced below. To find the optimum, PSO mimics the flight of a flock of birds and simplifies the social behavior with integrated individuals. PSO is a popular algorithm that has been used in many optimization problems [28]. In PSO, a population/swarm means a set of particles, and their corresponding position is one possible solution. The optimal direction for the next step of each particle is calculated based on the velocity according to best experiences of the particle itself and the swarm. The position, X d i , and the velocity, V d i , of the i-th particle in a population is recalculated as shown below.
where pbest is the best solution up to the current iteration of the particle itself, gbest is the best solution among pbests up to the current iteration, r i is a random number between 0 and 1, c 1 and c 2 are acceleration parameters for cognitive and social behaviors, and ω is the inertia weight. Deb's rules is one of the popular methods to handle the constraints, and it was adopted in the current study.

Results and Discussion
Indicated in Figure 16 is the fragility curve of the respective scouring depth under the same performance level. As shown, the failure probability was not proportional to the scour depth. In addition, the performance level of slightly damaged had the least failure probability. To show this trend in a clearer way, Figure 17 denotes the fragility curve of the respective performance level under the same scouring depth and water level. As shown, the fragility curves indicate that the deeper the scour depth, the higher the failure probability, and the lower the performance level, the higher the failure probability. Such conclusion is consistent with the concept held by general public reasoning that the strength of a bridge foundation will degrade after being scoured and eroded. However, varied water depth tends to affect the failure probability of the bridge. Indicated in Figure 18 are the fragility curves of each water depth under the same scouring depth (0 m) and the same performance level. Taking the slight damage level for example, the failure probability increased from 43% to 56% if the water depth increased from 0 m to 8 m at PGA of 0.2 g. As seen, the influence of water depth was on the nonconservation side, and the influence of water depth should be considered, especially, when longer piers are used.
Indicated in Table 1 is the result of errors and relative coefficients obtained during LSSVM training for the case when the scour depth was 0 m. It is shown that the prediction performance was not perfect, but it was acceptable. To be more specific, Lewis [29] suggested that if the MAPE is less than 10%, the prediction performance is considered as highly accurate. If the MAPE is between 11% to 20%, the prediction performance is good. After completing the LSSVM training, the corresponding result can be obtained by inputting the input parameters, such as water depth and PGA, in the trained machine to compute the fragility curve. Indicated in Figure 19 is the fragility curve under different scouring depths when the water depth was 3.5 m. Note that the purpose of LSSVM is to provide results that are not calculated by an authentic model such as SAP2000. For example, scour depths shown in Figure 19 were 3.25, 5, 6.25, and 7.5 m. By comparing Figures 16 and 19, it is learned that the LSSVM will provide a similar kind of failure probability trend. When viewing from a moderate limit state, the failure probability between 4 and 8 m of scouring depth, as indicated in Figure 16, was 0.50-0.70 (at 0.3 g), and that in Figure 19 was 0.52-0.72 (at 0.3 g) between 3.25 and 7.5 m of scouring depth. This indicated that the failure probability slightly increased when the immersion depth increased from 0 to 3.5 m. Similar conclusions can be inferred under other performance levels indicated in Figures 16 and 19. return period), seven analyses are needed. Such analyses often need several days or weeks to complete. If LSSVM is adopted, one can acquire the displacement within a few minutes. Indicated in Figure 21 is the joint failure probability distribution subjected to 100 year return period scouring when the bridge is immersed underwater at a depth of 4 m.      Through the aforesaid comparisons, we can learn that the fragility curves obtained from the fragility curve established by LSSVM and structural analysis are very close. The failure probability of multiple hazards can be calculated according to the scouring hazard curve (Figure 10), seismic hazard curve (Figure 12), and bridge conditional failure probability obtained from nonlinear pushover and time history analyses, given these specific hazards. A direct comparison of using two different approaches (authentic and LSSVM) are displayed as shown in Figure 20. It is seen that the trends of two approaches were similar. The greater the scour depth, the bigger the difference. Please note that the interpolated data between authentic ones may have a higher/lower exceedance probability (EP) than that of the authentic one, as shown in Figure 20. The proposed surrogate model is applied to Equation (23) only when the authentic model is not available. That is, if EP from the authentic model is unavailable, interpolated data from LSSVM was adopted to build plots in Figure 21. The main purpose of utilizing LSSVM for fragility analysis is to reduce the computational cost and retain the correct trends for engineering preliminary analysis. A failure decision for a given bridge depends not only on the failure probability but also on other factors such as the safety threshold, PGA design, scour depth, and performance level. Taking the bridge in Figure 20 for example, based on the bridge site, the PGA design was 0.24 g, assuming that the safety threshold was 0.00135, which is equivalent to a common reliability index (i.e., β = 3) used for infrastructure. Under these conditions, it was found that among four different scour depths with collapse performance limit states considered, the case of 8 m scour can be considered as a false negative event. Thus, if the predicted fragility is used accompanied with other factors. such as reliability target, scour depth, and PGA, for bridge safety judgement, caution should be taken. If computation time is an issue, LSSVM provides an alternative to shorten the analysis time. The effectiveness of the LSSVM approach depends on the complexity of a given bridge. If the authentic approach is used, many pre-analyses such as material nonlinear properties, added mass, soil spring, plastic hinge length and property, nonlinear pushover analysis, and time history analysis are needed to build the fragility curve. Taking nonlinear time history analysis in this study for example, for each seismic intensity (e.g., the 475 year return period), seven analyses are needed. Such analyses often need several days or weeks to complete. If LSSVM is adopted, one can acquire the displacement within a few minutes. Indicated in Figure 21 is the joint failure probability distribution subjected to 100 year return period scouring when the bridge is immersed underwater at a depth of 4 m.

Conclusions
In view that Taiwan is affected by typhoons and earthquakes all year long, this research proposes a process for analyzing immersed river-crossing bridges when subjecting to multiple hazards. In addition to the immersing effect, scouring depth and the seismic hazard are included in failure probability calculations. According to the research results, the following important conclusions are drawn: When calculating the added mass inferred from Morison equation, the Cm and Cd coefficients are from the semiempirical formula, and both will be affected by the Reynolds number and stream velocity. Because the velocity employed in this research was not so fast, it can be assumed as a constant. II. Based on the variation of different return years, the scouring depth increases along with the

Conclusions
In view that Taiwan is affected by typhoons and earthquakes all year long, this research proposes a process for analyzing immersed river-crossing bridges when subjecting to multiple hazards. In addition to the immersing effect, scouring depth and the seismic hazard are included in failure probability calculations. According to the research results, the following important conclusions are drawn: When calculating the added mass inferred from Morison equation, the C m and C d coefficients are from the semiempirical formula, and both will be affected by the Reynolds number and stream velocity. Because the velocity employed in this research was not so fast, it can be assumed as a constant. II.
Based on the variation of different return years, the scouring depth increases along with the increase of return year. III.
The immersing effect is not so apparent as expected when assessing its impact to the safety of a river-crossing bridge, and if the pier is not long enough, the immersing effect may not be so prominent. However, the influence of immerse effect on the bridge is on the non-conservation side, it is important to know that the immersing effect should not be ignored when the bridge pier is long. IV.
This research uses the displacement ductility (R) as the performance measurement, and such displacement ductility is the ratio between ultimate displacement and yielding displacement. It is seen that the yielding displacement plays a significant role in safety measurements. In light of this, this study proposes a way to calculate the yielding displacement. V.
The error indicator of LSSVM is very ideal, and the trend of fragility curves drafted according to such results is very close to that established with the authentic model. If computation time is an issue, LSSVM provides an alternative to shorten the analysis time for fragility analyses without scarifying the accuracy.

Conflicts of Interest:
The authors declare no conflict of interest.