Modelling the sequential earthquake–tsunami response of coastal road embankment infrastructure

Transport networks in coastal, urban areas are extremely vulnerable to seismic events, with damage 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 cyclic loading is modelled, applied to the soil-structure system using a finite-difference approximation to determine the differential settlement, 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–tsunami event.

1 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 vulnerability of transport infrastructure located in coastal urban areas associated with the cascade effect of earthquake, tsunami, flooding and ground motion (Goda, 2021; Goda et al., 2017, 2021; Williams et al., 2019, 2020; Sarkis et al., 2018; Bhattacharya et al., 2017; Rowell and Goodchild, 2017; Koshimura et al., 2014). 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 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; 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).

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, and 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 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. 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 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 velocities of tsunamis are particularly important for 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, 2017a) 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 modelling methods, such as finite-volume method (Jose et al., 2017) and smooth particle hydrodynamic method (Altomare et al., 2015; Klapp et al., 2020; Hasanpour et al., 2021).

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 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 these elements are based on empirical data or expert judgement, 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, better yet, by implementing an analytical framework similar to the drivers, exchanges, state of the environment, consequences and responses (DESCR) framework proposed by Escudero Castillo et al. (2012) and 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 for the effective stress paths and corresponding soil 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 for the 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 the 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

The sequential model proposed is comprised of three steps (Fig. 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 an urban road, in the coastal city of Manzanillo, Mexico.

Figure 1 Evaluation of the sequential earthquake–tsunami effects on transportation infrastructure.

3 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 towards and below the continent on the North American Plate (Fig. 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 M_w=8 earthquake 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.

Figure 2 Seismic environment of the study area. Red stars indicate earthquakes M≥8.0; blue stars indicate earthquakes $\mathrm{7.0}\le M\le \mathrm{8.0}$.

The urban area of Manzanillo and its commercial port are connected to the rest of Mexico via the Colima–Manzanillo, Manzanillo–Puerto Vallarta and Cuyutlán–Manzanillo highways. Part of the embankment on one of these main avenues was selected for the analyses performed in this research. The Boulevard Miguel de la Madrid runs parallel to Miramar Beach, which is on Santiago Bay, northwest of the city of Manzanillo (Fig. 3). 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 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 Playa de Oro International Airport. This road is therefore a vital part of the infrastructure of the city, facilitating many socioeconomic activities.

Figure 3 Aerial views of Boulevard Miguel de la Madrid.

The data on the topography and topology of the section of embankment analysed here were obtained in fieldwork carried out in 2021. A total of 24 topographic profiles along a 2 km long section were surveyed (Fig. 4). The surveying equipment (eSurvey E300 PRO + P8II GNSS; global navigation satellite system) was positioned at a control point to later link the collected data with the official datum. With a centimetre 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 geometric characteristics of these cross-sections and their relative importance, the cross-section shown in Fig. 5 was selected for the sequential earthquake–tsunami analysis. The soil-embankment system was modelled using a three-dimensional approach where the 20 m segment of the coastal profile selected was considered uniform. This was the section considered to have the most critical condition of the road due to the height and the slope of the embankment.

Figure 4 Location of the road embankment analysed (© Google Earth).

Figure 5 (a) Cross-sections taken along the road and (b) the cross-section used in the analysis.

3.1 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. Sand deposits are found along the coast, and alluvial deposits are found around the igneous deposits (Fig. 6). Based on geophysical explorations at points S1 to S5 (Fig. 6), Tonatiuh Dominguez et al. (2017) show an in situ characterization of some dynamic properties, such as shear wave velocities and the shear wave velocities obtained from exploration at site 4, which is considered representative of characterizing sites with soil type B, where the section of the road embankment studied is located.

Figure 6 Survey sites S1–S5. Geotechnical characterization of Manzanillo and Santiago bays obtained from Tonatiuh Dominguez et al. (2017) and shear wave velocity at site S4.

Manzanillo is located in the subduction zone between the Rivera and Cocos tectonic plates with the North American Plate. The rate of displacement at the contact surfaces is about 50 to 70 mm yr⁻¹ (DeMets and Wilson, 1997), so the probability of severe earthquakes is high. Furthermore, in subduction zones, there is a high likelihood of tsunami occurrence, given the vertical displacement generated on the seabed, which can be several metres, extending over tens of thousands of square kilometres. Over time, numerous serious earthquakes have been recorded in this area, causing high-intensity ground movements, both in the continental region and along the coastline, as well as tsunami damage on the coast. According to historical records, at specific sites in this area the wave heights associated with tsunamis have reached 10 m. Tsunamigenic events occurring in this region are of M_w>8 magnitudes, among which are the 1787, 1932, 1985 and 1995 earthquakes in the study area.

3.2 Tsunamigenic seismic scenario

The seismic features were characterized using records from the seismic station MZ01 (Fig. 7) located on sandy deposits. It is the only station in free-field conditions in Manzanillo that was operating during the 1995 event (M_w=8.0). Table 1 summarizes the characteristics of the event. Figure 7 shows the location of station MZ01 and the epicentre of the 1995 event. Figures 8 and 9 show the acceleration time histories and their corresponding response spectra for each parameter measured at station MZ01.

Table 1 Characteristics of the 1995 Manzanillo earthquake and the MZ01 station record.

PGA: peak ground acceleration.

Figure 7 Location of the seismic station MZ01 and the epicentre of the 1995 Manzanillo earthquake (M_w=8) event (© Google Earth).

Figure 8 Acceleration time histories at station MZ01 during the 1995 Manzanillo earthquake (M_w=8.0).

Figure 9 Response spectra from the time histories recorded at station MZ01 during the 1995 Manzanillo earthquake (M_w=8.0).

Figure 10 Acceleration, velocity and displacement time histories (a) without baseline correction and (b) with baseline correction.

To perform the analysis, the record of the N90E component showing the larger accelerations and pseudo accelerations was selected. The acceleration time history was baseline-corrected to avoid residual velocities or displacements at the end of the movement, which could be confused with permanent displacements. Figure 10a and b show the acceleration, velocity and displacement histories before and after correction, respectively.

To simulate wave propagation during the earthquake, the selected acceleration history (N90E) was deconvolved using the soil profile of site S4 (Fig. 6), with the code SHAKE (Schnabel, 1972), which determines the vertical propagation of shear waves in a horizontally stratified semi-infinite soil deposit. Subsequently, the deconvolved earthquake was propagated again upward to verify that the propagation of the deconvolved and the recorded history are consistent. Figure 11 presents the comparison between the recorded and calculated response spectrum.

Figure 11 Comparison between the recorded and calculated spectrum after the propagation of the deconvolved earthquake.

3.3 Step 1: seismic response

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 for 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 (Fig. 12).

Figure 12 (a) Normalized shear stiffness curves $G/G_{\max }$ and (b) damping curves λ.

The stress history was determined and transformed into an equivalent number of cycles of uniform shear stress. In this way, the intensity and duration of ground motions are considered, as well as the variations in shear stress with depth. For practical purposes, the equivalent number of cycles of shear stress τ_av can be estimated as 65 % of the maximum stress τ_max (Seed and Idriss, 1970) (Fig. 13).

Figure 13 (a) Maximum shear stress profile for the study site and (b) cyclic stress ratio profile, using the equation proposed by Seed and Idriss (1970) and the site response analysis.

The number of equivalent stress cycles N_eq depends on the magnitude and duration of the earthquake. Seed et al. (1975) applied a weighting procedure to establish a uniform number of cycles of shear stress N_eq (with an amplitude of 65 % of the maximum cyclic shear stress τ_cyc=0.65τ_max) that would produce an equivalent pore pressure, where the number of cycles of uniform shear stress increases with increasing magnitude of the earthquake. For the case study, an earthquake magnitude M_w=8.0 corresponds to an equivalent number of cycles of uniform shear stress equal to 21 (Fig. 14).

Figure 14 Equivalent number of uniform cyclic loading N_eq for different seismic magnitudes (Martin et al., 1975).

The liquefaction potential of the soil was evaluated using the cyclic stress criterion (Kramer, 1996; Seed and Idriss, 1970). One of the most common methods to evaluate liquefaction resistance is based on standard penetration tests (SPTs). This criterion uses the cyclic resistance ratio (CRR) from blow counts from the SPT, corrected to 60 % of the energy N₆₀. Thus, according to the shear wave velocity profile of the study site, N₆₀ was estimated using the expressions and parameters shown in Tables 2 and 3, obtaining N₆₀=27, for the first 10 m of depth.

Table 2 Correlations to estimate shear wave velocity.

N: SPT blow counts, N₆₀: blow counts 60 % energy corrected.

Table 3 Parameters for estimating N₆₀.

The safety factor against liquefaction for the embankment SF_l was calculated based on the critical strength ratio (CSR), which was derived from the site response analysis (Table 4). A correction factor for the magnitude of the earthquake was applied. The CRR was obtained using the Robertson and Wride (1998) approximation, and the magnitude correction factor (MSF), proposed by Andrus and Stokoe (1999), was considered through the following equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\mathrm{CRR}}_{\mathrm{7.5}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{34}-{\left(N_{\mathrm{1}}\right)}_{\mathrm{60}}}}+{\displaystyle \frac{{\left(N_{\mathrm{1}}\right)}_{\mathrm{60}}}{\mathrm{135}}}+{\displaystyle \frac{\mathrm{50}}{{\left[\mathrm{10}\cdot {\left(N_{\mathrm{1}}\right)}_{\mathrm{60}}+\mathrm{45}\right]}^{\mathrm{2}}}}\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{200}}},\end{array}\text{(2)}& {\displaystyle }\mathrm{MSF}={\left(M_{\mathrm{w}}/\mathrm{7.5}\right)}^{-\mathrm{2.56}},\text{(3)}& {\displaystyle }{\mathrm{SF}}_l=\left({\mathrm{CRR}}_{\mathrm{7.5}}/\mathrm{CSR}\right)\mathrm{MSF}.\end{array}$$

The preliminary result shows that in the first 10 m of the profile there is the potential for soil liquefaction (Tables 5 and 6).

Table 4 Safety factor against liquefaction.

CSR: critical strength ratio, CRR_7.5: cyclic resistance ratio, MSF: magnitude correction factor, SF_l: safety factor against liquefaction.

The seismic response of the embankment was evaluated using a three-dimensional finite-difference simulation. This was carried out using the Finn et al. (1977) constitutive model for liquefiable soils, available in the software programme FLAC^3D (Fast Lagrangian Analysis of Continua in 3 Dimensions) for the site response analyses of the slope and to account for soil non-linearities and loss of strength associated with pore pressure generation during cyclic loading (Fig. 15). For monotonic loading, an elastoplastic Mohr–Coulomb failure criterion was used that enables horizontal and vertical accelerations, displacements, shear forces, and pore pressure to be obtained. The model was composed of 13 632 elements and 16 083 nodes. The dimension of the element was selected based on the geometry and thicknesses of the soil layers. However, numerical distortion of the propagating wave can occur in a dynamic analysis as a result of the modelling conditions (Itasca Consulting Group, 2009). To represent wave transmission accurately, Kuhlemeyer and Lysmer (1973) suggest keeping the size of the spatial element Δ_l less than $\mathrm{1}/\mathrm{5}$ of the wavelength associated with the highest-frequency component of the input wave that contains significant energy f_max (i.e. ${\mathrm{\Delta }}_l\le \mathit{\lambda }/\mathrm{5}$), where the shortest wavelength λ is obtained from $\mathit{\lambda }=V_{\mathrm{s}}/f_{\max }$. In the case study, the smallest average shear wave velocity V_s of the site corresponds to that of the first stiff soil layer (i.e. shear wave velocity is around 260 m s⁻¹ in the first 10 m layer), while the highest excitation at which the energy is concentrated is about 1–5 Hz. Then, λ varies between 260 and 52 m approximately. The maximum spectral responses of the excitation occur even at higher frequencies (i.e. 0.5 and 10 Hz) as shown in Fig. 16. Hence, a Δ_l of 2 m was deemed appropriate. In previous research, using equivalent linear properties (Mayoral et al., 2015, 2016, 2017a, b), meshes with element sizes of 2 m showed good agreement between finite-difference models developed with FLAC^3D and SHAKE.

Figure 15 Soil profile and road embankment in the three-dimensional finite-difference model. Each colour represents soil layers according to depths indicated in Table 5.

Figure 16 Frequency content of the input ground motion.

Table 5 Soil characteristics.

C: soil shear strength, ϕ: friction angle, V_s: shear wave velocity, E: Young's modulus, G: shear stiffness, K: constraint modulus.

Table 6 Embankment characteristics.

C: soil shear strength, ϕ: friction angle, V_s: shear wave velocity, E: Young's modulus, G: shear stiffness, K: constraint modulus.

For calibration purposes, one-dimensional models were developed with the programme FLAC^3D to solve the 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 programme SHAKE (Fig. 17). 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 in 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 of cycle loads. With these assumptions an incremental constitutive relationship of the degradation curve can be described by ${\mathit{\tau }}_{\mathrm{n}}/\mathit{\gamma }=G/G_{\max }$, where τ_n is the normalized shear stress, γ is the shear strain and $G/G_{\max }$ is the normalized secant modulus. The sig3 model is defined using the following expression:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{G}{G_{\max }}}={\displaystyle \frac{a}{\mathrm{1}+\exp \left(-\frac{L-x_{\mathrm{0}}}{b}\right)}},\end{array}$$

where L is the logarithmic deformation, defined as L=log ₁₀(γ), and the parameters a, b and x₀, used by the 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 x₀ take the values 1.014, −0.50 and −1.25, respectively. Figure 17 shows a comparison between the curves used in the deterministic unidimensional model (Seed and Idriss, 1970) and those obtained with the sig3 model. Figure 18 shows a comparison between the response spectra calculated on the surface with the deterministic one-dimensional model, the equivalent linear FLAC^3D and the non-linear FLAC^3D. Good agreement between the results was observed.

Figure 17 Results obtained from Seed and Idriss (1970) and the sig3 model.

Figure 18 Response spectra at the surface with the different types of analyses.

For the three-dimensional model, the Finn–Byrne model was used to evaluate the liquefaction potential. This model is defined as

$$\begin{array}{}\text{(5)}& {\left(\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{vd}}\right)}_{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{ciclo}}=C_{\mathrm{1}}\exp \left(-C_{\mathrm{2}}\left({\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{vd}}}{\mathit{\gamma }}}\right)\right),\end{array}$$

where Δε_vd is the decrease in volume per half cycle of deformation, γ is the shear strain, and C₁ and C₂ are constants that can be calculated as

$$\begin{array}{}\text{(6)}& C_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{1}}^{\mathrm{c}};C_{\mathrm{2}}={\displaystyle \frac{\mathrm{0.4}}{C_{\mathrm{1}}^{\mathrm{c}};C_{\mathrm{1}}^{\mathrm{c}}=\mathrm{87}}}{\left(N_{\mathrm{1}}\right)}_{\mathrm{60}}^{-\mathrm{1.25}}.\end{array}$$

These parameters, presented in Table 7, were calculated using the average values of (N₁)₆₀, for the first 10 m of soil and the embankment.

Table 7 Coefficients used for the Finn–Byrne model.

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 Fig. 19).

Table 8 Sea level for each evaluated scenario.

Figure 19 Model monitoring points.

Seven control points were established in the model to obtain the soil behaviour parameters. Control points were located along the fault surface, at the toe of the embankment, at the highest point of the embankment and within the soil deposit (Fig. 19). Using the concept of static safety factor SF_s, prior the application of the seismic loading, the general stability of the embankment under static conditions was evaluated for each scenario. SF_s was above 2 in all the scenarios considered. The static safety factor was computed based on the concept of the strength over demand concept.

Next, an analysis of the seismic response was performed applying the N90E acceleration history (Fig. 11). It was observed that the largest vertical displacement, =73 cm, was recorded in scenario 5, on the crest of the embankment. Figure 20 presents the results of vertical displacement for scenarios 4 and 5 at the end of the earthquake. Figure 21 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.

At the control points, Fig. 22 shows the pore pressure ratio r_u obtained for each scenario. A rise in pore pressure was observed at points G and F, in all scenarios. However, with the increase in 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), which is based on solving the non-linear shallow-water equations through the numerical method of finite volumes, using adaptive mesh refinement to model small-scale features of the bathymetry as well as structures and coastal elements on the scale of metres (LeVeque et al., 2011). The shallow-water equations are the standard model used for transoceanic tsunami propagation as well as for local inundation (e.g. Yeh et al., 1994; Titov and Synolakis, 1995, 1998). In one space dimension these are

$$\begin{array}{}\text{(7)}& h_t+(hu{)}_x=\mathrm{0}\end{array}$$

$$\begin{array}{}\text{(8)}& (hu{)}_t+{\left(h{u}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}}\right)}_x=-gh{B}_x,\end{array}$$

where g is the gravitational constant, h(x,t) is the fluid depth and u(x,t) is the vertically averaged horizontal fluid velocity. The function B(x) is the bottom surface elevation relative to mean sea level. Where B<0, this corresponds to submarine bathymetry, and where B>0, it corresponds to topography. GeoClaw code implementation allows for the bathymetry and topography to be time-dependent by solving the two-dimensional shallow-water equations (LeVeque et al., 2011):

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{h}_t+(hu{)}_x+(hv{)}_y=\mathrm{0},\text{(10)}& {\displaystyle }& {\displaystyle }(hu{)}_t+{\left(h{u}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}}\right)}_x+(huv{)}_y=-gh{B}_x,\text{(11)}& {\displaystyle }& {\displaystyle }(hv{)}_t+(huv{)}_x+{\left(h{v}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}}\right)}_y=-gh{B}_y,\end{array}$$

where $u(x,y,t)$ and $v(x,y,t)$ are the depth-averaged velocities in the two horizontal directions and $B(x,y,t)$ is the topography.

Table 9 Fault characteristics of the 1995 Manzanillo earthquake.

Figure 20 Permanent vertical displacements for scenarios (a) 4 and (b) 5.

Figure 21 Permanent horizontal displacements for scenarios (a) 4 and (b) 5.

Figure 22 Variations in pore pressure ratios in scenarios (a) 1, (b) 2, (c) 3, (d) 4 and (e) 5.

The bathymetric and topographic information used was obtained from the GEBCO (General Bathymetric Chart of the Oceans) database, with a resolution of 15 arcsec. A mesh of 129 600 cells was used, applying three levels for mesh refinement, with the finest grids used near the embankment segment, where the grid resolution was 210 m. 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 h, beginning 15 min after the start of the earthquake. Figure 23 shows the simulation results for 1 min. A 1 h simulation period was found for the case study analysed according to records of the duration of the event regarding the wave arriving times (García et al., 1997; Borrero et al., 1997). However, for other cases, longer simulation times could be considered, such as those recommended by ASCE (Robertson, 2017). The authors acknowledge that the grid resolution of the propagation model is a possible research topic for the future. However, a higher resolution was not possible at the time the model was developed. The improvement in the grid spacing would help to reduce uncertainties in the expected flood elevations.

Figure 23 Sea surface elevation obtained after a simulation time of 1 min.

Wave speeds and elevations were obtained at five virtual gauges in Santiago Bay, near the road embankment studied. These are shown in Table 10 and Fig. 24. Figure 25 shows the change in bathymetry at the gauge locations and the free sea surface elevation, respectively. Figure 26 shows the free sea surface elevation, the wave velocity components u and v at the gauge nearest to the coast, and the wave speed s.

Table 10 Location of virtual gauges in the simulation.

Figure 24 Location of gauges used to obtain the free sea surface elevations and speeds in the simulation (© Google Earth).

Figure 25 (a) Ocean depths at the gauge locations (top) and (b) wave height distribution obtained during the simulation (bottom).

Figure 26 For the gauge nearest to the coast (gauge 5): (a) free sea surface elevation, (b) horizontal velocity u, (c) wave speed s and (d) vertical velocity v.

It is worth mentioning that the maximum flood levels obtained coincide with the visual data reported by Avila-Armella et al. (2005).

3.5 Step 3: earthquake–tsunami response

Applying a smooth particle hydrodynamics approach, wave-induced hydrodynamic pressures associated with tsunami inundation at the study site were determined over the cross-section of Fig. 5. The hydrodynamic conditions of the simulation reproduce the flow velocities and elevation of the free sea surface of indicator 5, the closest to the coast (Fig. 26). The smooth particle hydrodynamics (SPH) model used is DualSPHysics, which has been widely used for hydrodynamic forces and to model complex fluid flows (González-Cao et al., 2019; Ye et al., 2019).

3.5.1 Validation

To validate the numerical approach, the experimental and numerical model of St-Germain et al. (2012) was reproduced, comparing water elevations of the incident wave in two control gauges. A tank 13.17 m long, 2.7 m wide, divided into two sections, and 1.4 m in height (Fig. 27) was considered. The model is made up of a dam break, consisting of a block of water of 0.85 m height that is released at time 0 through a swinging gate to travel towards the outflow.

Figure 27 Experimental model of St-Germain et al. (2012).

Figure 28 Water heights at points W1 and W2 (St-Germain et al., 2012).

The validation SPH simulation parameters are shown in Table 11.

Table 11 Parameters and main characteristics of the SPH simulation.

The simulation results are water heights at points W1 and W2 (Fig. 28). The standard deviation of the DualSPHysics model with respect to the experimental model is 5.4 cm for W1 and 5.2 cm for W2, which are 6.38 % and 6.16 % with reference to the initial height of the water. The standard deviation of DualSPHysics with respect to the numerical model of St-Germain et al. (2012) is 2.9 cm for W1 and 3.5 cm for W2, 3.4 % and 4.1 % concerning the initial height of the water.

3.5.2 Analysis of inter-particle spacing sensitivity

A simulation resolution analysis was performed for the SPH model. The resolution used was the initial size of the lattice nodes of the fluid particles and the fixed particles, defined as inter-particle spacing (dp). The simulation was made through tests with five different inter-particle spacings to obtain the resolution that allows for a convergence in the results with a lower computational cost (runtime).

The hydrodynamic pressures at the position x=10 m, z=2.1 m of the different resolutions (dp = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 and 0.7 m) were compared. The agreement of the results obtained with different resolutions is quantified by taking as reference the data series of the finest resolution of dp = 0.1 m, by means of the normalized standard deviation (σ_n), the centred root-mean-square deviation (RMSD) and the correlation (R), as shown in the Taylor diagram of Fig. 29 (as applied in González-Cao et al., 2019, and Klapp et al., 2020).

Figure 29 Taylor diagram of hydrodynamic pressures at x=10 m and z=2.1 m, taking as reference the data series of the finest resolution of dp = 0.1 m.

Details of the resolution, runtime and number of simulated particles are shown in Fig. 30. The parameters and main characteristics of the simulations are shown in Table 12. The simulations were run on an Nvidia 1650 graphics processing unit (GPU).

Figure 30 Number of particles and runtime for the different resolutions.

Table 12 Parameters and main characteristics of the SPH simulation.

The results of the inter-particle spacing sensitivity analysis show that for all resolutions the correlation is greater than 95%, and that improves when the resolution is finer. However, a resolution greater than 0.1 m may not significantly improve the convergence of the results and may increase the computation time. Therefore, a resolution of 0.1 m with a particle number of 104 825 and a computation time of 0.94 h was selected.

Figure 31 shows the pressures acting on the road embankment, and Fig. 32 shows the pressure evolution at points a–f, obtained in the first 90 s of loading, corresponding to the wave colliding with the slope of the embankment.

Figure 31 Initial sea level and control points considered for the earthquake–tsunami loading analysis.

Figure 32 Loading variation stress as a function of time at points a–f.

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.

Figure 33 shows the hydrostatic and hydrodynamic loading conditions applied at 50 s. The resulting vertical and horizontal displacements and the safety factor (i.e. capacity over demand) obtained from steady-state loading conditions are shown in Figs. 34 to 36.

Figure 33 Loading conditions at points a–f, at a time of 50 s.

Figure 34 Vertical displacements contours in metres (a) after seismic loading and (b) after tsunami loading.

Figure 35 Horizontal displacement contours in metres (a) after seismic loading and (b) after tsunami loading.

Figure 36 Safety factor versus liquefaction (a) after seismic loading and (b) after tsunami loading.

The analysis shows that the main effects induced by wave loading are the horizontal displacement of the beach slope, reaching 0.8 m near point f. No major impacts were seen in the road embankment, where both vertical and horizontal displacements and the liquefaction safety factor show minor variations. However, there was a loss of granular material due to excessive ground deformation in the beach area. This effect was simulated following a time step approach, where zones of the mesh that reach a deformation of 2 m were removed from the model. The hydrodynamic loads were applied to the remaining mesh, which allowed the incremental unloading of the unstable areas in the sand found on the beach slope. The coupled effect of liquefaction and a tsunami could potentially lead to a loss of functionality of the road associated with an approximate 0.7 m displacement of the slope.

4 Conclusions

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 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 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 give a 20 % reduction in the height of the structure.

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 seismic loading causes the greatest effects in pore pressure increase and in ground failure due to vertical and horizontal displacements in the embankment.

The sequential approach presented allows for soil displacements and strength to be accurately quantified, as well as pore pressure increase derived from an earthquake. The effects also couple with the tsunami arrival, which is not captured in decoupled models. The evaluation of these potential cumulative impacts provides additional information for the design and planning of more sustainable and resilient transportation infrastructure. The method presented is applicable to any coast, as long as there is sufficient information to characterize the site and structures, such as the seismic environment, geotechnics, bathymetry and structural systems. The degree of detail of the information required is of great importance to reduce uncertainty in the results.

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 ignores 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 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 (Chávez et al., 2021; Silva et al., 2021). In the case of Manzanillo, a pile-driven road is probably a better option.