Modelling the sequential earthquake-tsunami response of coastal urban transport infrastructure

. Transport networks in coastal, urban areas are extremely vulnerable to seismic events, with damage 10 likely due to both ground motions and tsunami loading. Most existing models analyse the performance of structures under either earthquakes or tsunamis, as isolated events. This paper presents a numerical approach that captures the sequential earthquake-tsunami effects on transport infrastructure in a coastal area, taking into consideration the combined strains of the two events. Firstly, the dynamic cyclical loading is modelled, applied to the soil-structure system using a finite difference approximation to determine the differential settlement, 15 lateral displacement, and liquefaction potential of the foundation. Next, using a finite volume method approach, tsunami wave propagation and flooding potential are modelled. Finally, the hydrostatic and hydrodynamic loads corresponding to the wave elevation are applied to the post-earthquake state of the structure, to obtain a second state of deformation. The sequential model is applied to an embankment in Manzanillo, Mexico, which is part of a main urban road, the response is analysed using ground motion records of the 1995 Manzanillo earthquake- 20 tsunami event. the stage at the end of the earthquake and the stage at the end of tsunami-induced flooding, is seen on the beach slope. In the case study, the sequential model shows that 385 seismic loading causes the greatest effects in pore pressure increase and in ground failure due to vertical and horizontal displacements in the embankment.


Introduction
Transportation networks are key elements of the economy and society; therefore, more attention has been paid to the need for sustainable, adaptive and resilient transport systems under natural hazards and other emergency conditions in the last decades. Past experience of the subduction mechanisms of seismic events has shown the 25 vulnerability of transport infrastructure located in coastal urban areas associated with the cascade effect of earthquake-tsunami-flooding-ground motion (Bhattacharya et al., 2017;Koshimura et al., 2014;Rowell and Goodchild, 2017;Sarkis et al., 2018;Williams et al., 2020). Existing methodologies to analyse the response of transport assets during earthquake and tsunami conditions consider these phenomena as independent events.
Probabilistic analyses and fragility curves (Akiyama et al., 2020;Burns Patrick et al., 2021) are examples of 30 multi-hazard approaches that aim to determine the occurrence probability of both phenomena simultaneously and the expected level of damage of structures after them. Recent multi-hazard literature has focused on structures in the urban environment such as buildings and bridges (Attary et al., 2021;Burns Patrick et al., 2021; https://doi.org/10.5194/nhess-2021-381 Preprint. Discussion started: 2 February 2022 c Author(s) 2022. CC BY 4.0 License. Ishibashi et al., 2021;Karafagka et al., 2018), while little attention has been paid to the response of road networks to sequential earthquake-tsunami processes (Iai, 2019).

35
The impacts of a combined earthquake-tsunami on urban transport networks have important economic and social effects in the short, medium and long terms, in both human and material losses. These impacts are fundamentally associated with the loss of connectivity of regions or cities, which primarily disturb the continuity of rescue efforts, reconstruction, urban logistics, health and safety activities (Dizhur et al., 2019;Sarkis et al., 2018). Most of the negative effects are associated with partial or total damage to urban and 40 interurban roads that disrupt land communication between key sites (Cubrinovski, 2013). In coastal areas, one significant factor causing damage in earthquakes is the increase in the pressure within the soil mass induced by the seismic loads, and by the accumulation of hydrodynamic pressure in the non-cohesive saturated soils common in these areas. From observations of past events in coastal areas it is seen that moderate seismic intensity can cause liquefaction, resulting in reduced stiffness and a loss of shear strength in liquefied soils.

45
This can lead to soil settlement, increased lateral ground pressure, and loss of strength (Kakderi and Pitilakis, 2014), bringing about damage to structures, sometimes even causing their collapse, and in the case of embankments, potential slope-stability problems, or even landslides.
Regarding tsunami damage, one of the most critical aspects to consider is the hydrodynamic behaviour of the wave as it approaches the coast. However, information on the physical parameters that characterize these waves 50 is often limited, due to difficulties in achieving accurate measurements during the event.
The impact of a tsunami on the coast is governed by non-linear physics, such as turbulence (Klapp et al., 2020).
The hydrodynamic behaviour during a tsunami is complex and can significantly affect structures. Forces induced by the hydrodynamics associated with a tsunami are governed by various fluid parameters, such as density, velocity, and depth, as well as the geometry of the structure. These forces induced by the pressures and 55 velocities of tsunamis are particularly important in the stability of coastal structures, affecting assets such as piers, embankments and bridge piles (Chinnarasri et al., 2013). Numerical approaches based on finite element methods (Argyroudis and Kaynia, 2015;Mckenna et al., 2021) and finite differences (Mayoral et al., 2016;Mayoral et al., 2017) have been applied to assess the seismic vulnerability and potential damage of transport infrastructure. Likewise, the effects of hydrodynamic loads on structures have been estimated through 60 modelling methods, such as finite volume method (Jose et al., 2017) and smooth particle hydrodynamic method (Altomare et al., 2015;Klapp et al., 2020).
Besides the simulation models available to assess the individual effects of earthquake and tsunami loads on transport infrastructure, a few, limited, modelling approaches to predict sequential damage and effects are available. From a review of the literature on fragility and vulnerability models it was seen that most studies 65 focus on individual transport components or networks, usually considering only one hazard at a time (Argyroudis and Kaynia, 2015;Argyroudis and Mitoulis, 2021;Briaud and Maddah, 2016;Iai, 2019;Maruyama et al., 2010;Nibs, 2004). The studies mentioned focus mainly on the vulnerability of bridges and tunnels, and the main emphasis is on ground movement due to seismic excitation. Research based on simplified models has examined other structures such as embankments, slopes, retaining walls, and abutments. Existing models for 70 https://doi.org/10.5194/nhess-2021-381 Preprint. Discussion started: 2 February 2022 c Author(s) 2022. CC BY 4.0 License. these elements are based on empirical data or expert judgment, mainly for seismic analysis, and focus on a universe of limited typologies.
There is, at present, the opportunity to improve the modelling approaches that consider more structure types, as well as multiple hazards, including earthquakes and tsunamis simultaneously. To extend these studies, the risk analysis can be extended by using one of the methodologies described by (Escudero Castillo et al., 2012) or 75 better yet, by implementing an analytical framework similar to the Drivers-Exchanges-State of the environment-Consequences-Responses (DESCR) framework proposed by (Escudero Castillo et al., 2012;Silva et al., 2020) in which ecosystem-based adaptation and community-based adaptation is considered (Silva et al., 2019). This paper proposes a numerical approach to capture the sequential earthquake-tsunami effects on urban road embankments located in coastal areas. The approach allows the effective stress paths and corresponding soil 80 strain evolution during the event to be considered, taking into account soil-embankment response, tsunami propagation, and hydrodynamic behaviour. The first step consists of computing the initial static stress and deformation conditions, and modelling the dynamic cyclic loading applied to the soil-embankment system in an earthquake. A finite difference approximation is used to determine the differential settlements and lateral displacements, associated with excess pore pressure generation during cyclic loading. This allows the 85 assessment of the liquefaction potential of the soil foundation. Next, a finite volume approach and smooth particle hydrodynamics are combined, to model tsunami wave propagation at oceanic level and flooding potential. Finally, the resulting sea-surface elevation and hydrodynamic loads are applied to the post-earthquake state of the structure to obtain the final deformation state.

2.
Sequential earthquake-tsunami response model on transport infrastructure 90 The sequential model proposed is comprised of three steps, Figure 1. The first step is the numerical analysis of the seismic response of the soil-structure system. The next stage corresponds to the wave generation and propagation model from deep, offshore waters up to the coast. Finally, there is the analysis of the soil-structure system response to hydrodynamic loadings, given the state of deformation induced by the seismic excitation.
This methodology is applied to an embankment of urban road, in the coastal city of Manzanillo, Mexico.

Case study
Manzanillo, on the Mexican Pacific Coast, is the most important commercial port in the country. The zone is seismically complex, due to the triple junction of the Cocos and Rivera subduction plates that are moving 100 towards and below the continent, on the North American Plate, Figure 2. The city has been affected by many major earthquakes, and several in the last century (Tonatiuh Dominguez et al., 2017;Eqe International, 1996;Ovando-Shelley and Romo, 2004). In 1995, an earthquake Mw = 8 caused significant damage, mainly in the coastal area of the city, where collapsed buildings were reported and liquefaction of the saturated sandy soils was observed at several locations in the port area of Manzanillo. This beach has attractive natural characteristics, perfect for practicing sun and beach recreation activities (Cervantes et al., 2015); the road has become a waterfront promenade, with several restaurants and access points to the beach along the road. This road also connects the communities of El Naranjo, in the west, and Colinas de 115 Santiago, in the east. It is also an important public transport route, serving tourists and the labour force going to and from the restaurants, hotels and recreation facilities. Further west is Manzanillo International Airport.
This road is therefore a vital part of the infrastructure of the city, facilitating many social-economic activities.

120
The data on the topography and topology of the section of embankment analysed here were obtained in field work carried out in 2021. A total of 24 topographic profiles along a 2 km long section were surveyed, Figure 4.
The surveying equipment eSurvey E300 PRO + P8II GNSS was positioned at a control point to later link the collected data with the official datum. With a centimeter precision, post-processing of the data gave the coordinates of a geodesic point which served as a reference to carry out the topographic survey. Based on the 125 geometric characteristics of these cross-sections and their relative importance, the cross-section shown in Figure   5 was selected for the sequential earthquake-tsunami analysis.

Seismic and geotechnical characterization of the study area
Manzanillo is in the northern part of the Sierra Madre del Sur mountain range. Geological studies by (Tonatiuh Dominguez et al., 2017) state that granite intrusions, extrusive igneous outcrops and limestone deposits are found in the region. Deposits of clay and organic materials are also found in the two saltwater lagoons nearby.

135
Sand deposits are found along the coast and alluvial deposits are found around the igneous deposits, Figure 6.
Based on geophysical explorations at points S1 to S5, Figure 6, (Tonatiuh Dominguez et al., 2017) shows an in situ characterization of some dynamic properties, such as shear wave velocities. Figure 7 depicts the shear wave velocities obtained from exploration at Site 4, which is considered representative for characterizing sites with soil type B, where the section of the road embankment studied is located.

Tsunamigenic seismic scenario
The seismic features were characterized using records from the seismic station MZ01, Figure 8, located on sandy deposits. It is the only station in free-field conditions in Manzanillo that was operating during the 1995 event (Mw = 8.0). Table 1 summarizes the characteristics of the event. Figure

180
To simulate seismic wave propagation, during the earthquake, a site response analysis was performed using the N90E acceleration history. The record was deconvolved and propagated through the soil profile of the S4 site, performing a one-dimensional deterministic analysis. In this way, a first approximation of the expected dynamic amplification in the soil deposit was obtained, which is useful prior to the non-linear analysis in the time domain.
This analysis also allows the maximum dynamic shear stresses in the soil to be obtained, which are used to estimate the liquefaction potential in the simplified methods. Since dynamic laboratory tests were not performed on samples obtained near the study area, it was considered appropriate to assume the damping curves and shear stiffness modulus ratio presented by (Seed and Idriss, 1970) for sands, Figure 12.   The liquefaction potential of the soil was evaluated using the cyclic stress criterion (Kramer, 1996;Seed and 205 Idriss, 1970). One of the most common methods to evaluate liquefaction resistance is based on standard penetration tests (SPT). This criterion uses the cyclical resistance ratio CRR from blow counts from the SPT test, corrected to 60% of the energy N60. Thus, according to the shear wave velocity profile of the study site, N60 was estimated using the expressions and parameters shown in Tables 2 and 3, obtaining N60=27, for the first 10 m of depth.

215
The safety factor against liquefaction for the embankment SFl, was calculated based on the critical strength ratio CSR, which was derived from the site response analysis. A correction factor for the magnitude of the earthquake was applied. The CRR was obtained using the (Robertson and Wride, 1998) The preliminary result shows that in the first 10 m of the profile there is the potential for soil liquefaction. The seismic response of the embankment was evaluated using a three-dimensional finite difference model, Figure 15, considering a depth of 110 m in the soil profile, as well as the parameters of site S4. Tables 5 and 6 show the soil profile and the embankment characteristics, respectively.   For calibration purposes, one-dimensional models were developed with the program FLAC3d, to solve the 240 equation of motion in the time domain, considering both equivalent linear and non-linear properties, and the results were compared with those obtained in the frequency domain, with the program SHAKE (Fig. 16). The sig3 model available in FLAC was used to approximately represent soil non-linearity. Thus, the sig3 hysteretic model was used to address the variation of the stiffness modulus and the damping ratio during the seismic event.
This model considers an ideal soil, in which the stress depends only on the deformation, and not on the number 245 of cycles loads. With these assumptions an incremental constitutive relationship of the degradation curve can be described by τn / γ = G / Gmax, where τn is the normalized shear stress, γ is the shear strain, and G/Gmax is the normalized secant modulus. The sig3 model is defined using the following expression: where L is the logarithmic deformation, defined as L = log10 (γ), and the parameters a, b and x0, used by the 250 sig3 model, were obtained through an iterative process, in which the modulus degradation curves were fitted with the model. The corresponding damping is given directly by the hysteresis loop during cyclic loading. For the case study, the parameters a, b and x0 take the values 1.014, −0.50 and −1.25, respectively. Figure 16 shows a comparison between the curves used in the deterministic unidimensional model (Seed and Idriss, 1970), and those obtained with the sig3 model. Figure 17 shows a comparison between the response spectra calculated on

260
For the three-dimensional model, the Finn-Byrne model was used to evaluate the liquefaction potential. This model is defined as: where Δεvd is the decrease in volume per half cycle of deformation; γ is the shear strain, and C1 and C2 are constants that can be calculated as: These parameters, presented in Table 7, were calculated using the average values of (N1)60, for the first 10 m of soil and the embankment. Five scenarios were considered to assess the seismic performance of the road embankment, considering the variations in sea level registered in the study area, as well as the case of high tide, and soil saturation due to rain, Table 8 and Figure 18.   Next, an analysis of the seismic response was performed applying the N90E acceleration history, Figure 11. It 285 was observed that the largest vertical displacement, = 73 cm, was recorded in scenario 5, on the crest of the embankment. Figure 19 presents the results of vertical displacement for scenarios 4 and 5 at the end of the earthquake. Figure 20 gives the results of the horizontal displacement for the same cases. The largest horizontal displacement was observed in scenario 5, = 51 cm. Important vertical displacements of the soil were observed, as well as a tendency of lateral displacement of the body of the slope, increasing with depth.

290
At the control points, Figure 21 shows the pore pressure ratio, ru, obtained for each scenario. A rise in pore pressure was observed at points G and F, in all scenarios. However, with the increase of the sea level, the higher pore pressure ratio at point D, located at the toe of the slope, continues to increase in scenarios 3 and 4, and reaches values above 0.7 in scenario 5, where an additional increase in this ratio was observed at point E, located on the failure surface of the slope.

3.4.
Step 2: simulation of tsunami and wave propagation The tsunami simulation was carried out using the model implemented in the GeoClaw code (Berger et al., 2011),

300
which is based on solving the non-linear shallow water equations through the numerical method of finite volumes. The bathymetric and topographic information used was obtained from the GEBCO database generated by the NOAA, with a resolution of 15 arc seconds. Considering the characteristics of the fault mechanism of the 1995 Manzanillo earthquake, Table 9, the (Okada, 1995) fault model was used to estimate the vertical displacement on the seabed caused by the seismic event. Based on the calculated deformations and the characteristics of the earthquake, a tsunami-wave propagation model was run for a simulation period of 1 hour, beginning 15 minutes after the start of the earthquake. Figure   22 shows the simulation results for time 1 min.  Wave speeds and elevations were obtained at 5 virtual gauges in Santiago Bay, near the road embankment studied. These are shown in Table 10 and Figure 23. Figure 24 shows the change in bathymetry at the gauge locations, as well as the free sea surface elevation, respectively. Figure 25 shows the free sea surface elevation,

315
and u and v the wave velocity components at the gauge nearest to the coast, as well as the wave speed s.   It is worth mentioning that the maximum flood levels obtained coincide with the visual data reported by Avila et al., (2005).

3.5.
Step 3: earthquake-tsunami response Applying a smooth particle hydrodynamics approach, the wave-induced hydrodynamic pressures associated 330 with tsunami flooding at the study site were determined from the speeds and free sea surface elevation from   Table 11. The occupied Kernel function is the Wendland function, the time-step algorithm used is Symplectic, the viscosity scheme used is the artificial viscosity (α = 0.01) and the inter-particle spacing was 0.1 m. The simulations were run on NVIDIA 1650 GPU.

345
The wave-induced behaviour at the embankment was analysed, applying the same time-step numerical approach of Step 1, and the conditions of Scenario 5, with an initial maximum tide of 1.28 m and saturated soil. From the results obtained in Steps 1 and 2, the corresponding wave elevation of 1.5 m was modelled. In the model, the hydrodynamic pressures, pore pressures and incremental hydrostatic vertical loading were applied until the steady state in each of the loading points was reached.

370
A sequential methodology was implemented to analyse the response of transportation infrastructure in an earthquake-tsunami. The methodology was applied for Manzanillo, Mexico, where the behaviour of a section of a road embankment was analysed, considering the accumulated effects of the earthquake in terms of horizontal and vertical displacements, and those induced by an increase in sea level after the earthquake. In the case study, different scenarios were analysed depending on the initial elevation of the sea. The seismic response 375 of the embankment showed that the most critical condition occurred in the scenarios of maximum high tide, and maximum high tide with saturated soil. In these scenarios the vertical displacements were 0.57 m and 0.73 m, respectively, both of which occurred at the embankment crest. The vertical and horizontal displacements can be interpreted as a failure of the embankment, in the sense that they may imply loss of functionality, since these gives a 20% reduction in the height of the structure.

380
The simulation of wave propagation of a tsunami showed that the level of the free sea surface could increase 1.5 m in the Manzanillo Bay area (Avila-Armella et al., 2005). Further analysis of the earthquake-tsunami behaviour showed that for the scenario of maximum high tide and saturated soil, the highest variation in the horizontal displacements in the ground between the stage at the end of the earthquake and the stage at the end of tsunami-induced flooding, is seen on the beach slope. In the case study, the sequential model shows that 385 seismic loading causes the greatest effects in pore pressure increase and in ground failure due to vertical and horizontal displacements in the embankment. From this study it is clear that, as in other parts of the world, the road network in Manzanillo has been built to improve communications, but ignoring many coastal processes, such as the potential presence of a tsunami.
Moreover, this infrastructure, was built on coastal dune ridges, inducing coastal squeeze, isolating ecosystems 390 and increasing risk. The planning, design and construction of infrastructure in coastal areas susceptible to earthquakes and tsunamis should be rethought in order to make it compatible with natural processes and thus truly contribute to a better quality of life Silva et al., 2021). In the case of Manzanillo, a pile-driven road is probably a better option.