Tsunami run‐up along the Zihuatanejo region, Mexican Pacific Coast: A modelling study

A numerical hydrodynamic model was applied to Zihuatanejo Bay on the Mexican Coast of the Pacific Ocean to study the heights and arrival times of tsunamis produced by hypothetical 7.0 and 8.0 Mw magnitude earthquakes with epicenter in the Cocos Plate, at 16.9718°N–101.973°W, as well as the flooding of Zihuatanejo Town by the tsunamis. The travel times of the tsunamis to Zihuatanejo Coast are between 18 and 19 min, maximum sea surface elevations of 1 and 10 m at the coast are estimated following a 7.0 and 8.0 Mw earthquake respectively with a tsunami inland penetration up to 400 m in the second case affecting restaurants, hotels, banks, schools, businesses, museums, and even the municipal market of the town. Therefore, a densely populated area would be affected by a tsunami triggered by an 8.0 Mw earthquake.


| INTRODUCTION
Zihuatanejo is a tourist town located on the Mexican Pacific Coast (Figure 1), in front of the Orozco Fracture Zone that separates the Cocos and Rivera Plates (Nava, 1998). Therefore, this is an area with high seismicity (Eissler & McNally, 1984). The Catalogue of Tsunamis in the Western Coast of Mexico (Sánchez & Farreras, 1993) documents that 49 tsunamis arrived on the west coast of Mexico from 1732 to 1985, 16 of which were considered to have a distant origin, and 33 were considered to have a local origin. The wave heights of the recorded or observed distant-origin tsunamis have not exceeded 2.5 m, suggesting that the flood risk of distant tsunamis is low (García & Suárez, 1996).
Almost half of the local-origin tsunamis that occurred after 1952 caused considerable destruction. The tsunami that occurred in November 1925 reached a maximum height of 11 m in Zihuatanejo Town. Therefore, tsunamis of local origin represent a relevant hazard to the Mexican Pacific Coast due to their greater impact (Centro Nacional de Prevención de Desastres; CENAPRED, 2014).
The aim of this study is to estimate the wave heights, arrival times, form, and flood produced by tsunamis in the Zihuatanejo Area on the Mexican Coast of the Pacific Ocean caused by 7.0 and 8.0 Mw earthquakes with epicenters at 16.9718 N-101.973 W (Figure 1). Zihuatanejo was selected as study area because the literature shows that the region has been historically affected by tsunamis (Table 1) (SMN, 2017), respectively. Moreover, Zihuatanejo is a tourist region with high population growth, thus, a tsunami represents a potential risk of human and material losses.

| STUDY AREA
Zihuatanejo is a semi-elliptical bay, with a mouth width of 1 km, width of 1.5 km, and length of 2 km ( Figure 1). The Bay of Zihuatanejo contains four beaches: "Playa Principal" located in the downtown area of Zihuatanejo, which is 400 m long. The main port is located on this beach; it is also the most populated area and receives most of the tourism. "Playa La Madera" is only 300 m long and receives small to moderate wind waves. "Playa La Ropa" is the longest in the bay, extending for 1 km, and is characterized by mediumheight wind waves. Finally, "Playa Las Gatas" is the smallest beach in the bay, with a length of 250 m, and is located away from the city; it receives almost no swell due to the construction of a breakwater by pre-Hispanic natives (Ixtapa-Zihuatanejo, 2015). There is a shallow lagoon west of "Playa Principal" called Salinas lagoon. Furthermore, there is a water channel that crosses a large part of the city between "Playa Principal" and "Playa La Madera." The beaches in Zihuatanejo Bay on average have a slope of 5.2%, and a height of 1.2 m. The beaches are composed of 42% fine-grained sand, 50% medium-grained sand, and 8% coarse-grained sand. This indicates that are low energy beach systems. The peak period in the region is between 15 and 18 s, with a peak direction of SW, and SSW. Only the low energy wave reaches inside the Zihuatanejo Bay. In Zihuatanejo Bay, the sea/land breezes and the mesoscale winds produce wind waves with a peak period of 7 s, and a wave height between 0.5 and 1.0 m. The mean tidal range is 1.1 m (Secretaría de Turismo; SECTUR, 2014).
The bathymetry in the region ( Figure 2) reaches a maximum depth of 4,000 m in the Cocos Plate area and a depth of 25 m at the entrance to Zihuatanejo Bay.
Only one study has shown the computed areas of inundation due to a possible tsunami in Zihuatanejo, which was by Farreras, Ortiz, and Ramírez (1999). These results appear in other publications like (CENAPRED, 2001;Farreras & Ortiz, 2007). They classified the different flood F I G U R E 1 Study area. (a) Tectonic setting along the coasts of southern Mexico (Modified from Ramírez & Urrutia, 1999), location of Zihuatanejo Bay, epicenter of the two hypothetical earthquakes considered in the present study, and epicenter of the 1985, 1995, 2003, and 2017 earthquakes. (b) Main beaches of Zihuatanejo Bay T A B L E 1 Information about tsunamis with a local origin observed or registered (*) in Zihuatanejo or Western Coast of Mexican (CENAPRED, 2014;CMT, 2020;Lander, Whiteside, & Lockridge, 2003;Sánchez & Farreras, 1993;SMN, 2014bSMN, , 2017UNAM Seismology Group, 2015;U.S. Geological Survey [USGS], 2003Zobin, 1997) Date (

| MODEL FOR SIMULATING TSUNAMI GENERATION AND FLOODING
There are several models to study tsunamis; some of the most used models are TUNAMI code (Goto & Ogawa, 1997;Imamura, Yalciner, & Ozyurt, 2006), MOST (Method Of Splitting Tsunami) which are used by National Oceanic and Atmospheric Administration (NOAA) (Titov & Gonzales, 1997), COMCOT (gridcoupled tsunami model) (Liu, Yeh, & Synolakis, 2008;Wang, 2009), and Delft3D-FLOW (Robke, Schuttrumpf, & Vott, 2018). In addition to TUNAMI, MOST, COMCOT, and Delf3D-Flow, also Mike21 Flow Model FM developed by DHI has been used for tsunami simulation (Gayer, Leschka, Nöhren, Larsen, & Günther, 2010;Masina, Archetti, Besio, & Lamberti, 2017). In previous tsunami-modelling literature, several authors agree that the tsunami propagation phase should be described as though it were the movement of long waves and applied the shallow water theory in their study (Fernádez, Ortiz, & Mora, 2004;Goto & Ogawa, 1997;Imamura et al., 2006;Kowalik, 2012;Liu et al., 2008;Shuto, 1991;Titov & Gonzales, 1997;Wang, 2009). The shallow water theory correctly represents the propagation of tsunamis at depths greater than 50 m, although the theory is not accurate near the coast, where wave breaking occurs (Bryant, 2008). In this approximation the height of the wave compared to its length is considered to be very small, the vertical acceleration of the water particles is considered to be negligible compared with the gravitational force, and the curvature of the trajectory of water particles is considered to be very small. Consequently, the vertical movement of water particles has no effect on the pressure distribution, allowing the use of the hydrostatic approach, and the horizontal velocity of the water particles is considered to be vertically uniform (Goto & Ogawa, 1997). When analyzing the equations that all the models mentioned above use to describe the propagation of tsunamis, they are very similar and the only differences are that some use as variable the velocity (Fernádez et al., 2004;Titov & Gonzales, 1997) while others use the variable transport (Goto & Ogawa, 1997;Imamura et al., 2006;Liu et al., 2008;Wang, 2009). In addition, some models were developed using spherical coordinates, which is necessary to study regional and distant tsunamis because, at those scales, the propagation of the tsunami is affected by the Earth's curvature, and other models use Cartesian coordinates to simulate local tsunamis. The greatest differences between the models are the discretization of the equations for the numerical solution and their grid domains.
In this study, we used the MOHID (Hydrodynamic Model) with a MOHID Studio Professional License, Academic Usage (Mohid Mohid Studio, 2019) to simulate the propagation of a tsunami and the subsequent flooding. This numerical model uses finite volume approximation to discretize the primitive equations of momentum and continuity in a structured grid. In this approach, the discretized form of the equations is applied to the cell used as the control volume. Therefore, the actual form of solving the equations is independent of the cell geometry and generic vertical coordinates can be used (Martins, Leitao, Silva, & Neves, 2001). The equations are horizontally discretized over a staggered Arakawa-C grid. Temporal discretization is conducted by a semi-implicit alternating F I G U R E 2 Bathymetry of the study zone. (a) Regional area.

| Digital elevation model and bathymetry
To generate the terrain model of Zihuatanejo and its environs used in flood modeling, high-resolution digital terrain model derived from LIDAR data (Morelli, Pazzi, Frodella, & Fanti, 2018) with horizontal and vertical resolutions of 5 and 0.01 m, respectively were used, which were property of the INEGI (National Institute of Statistics and Geography of México). The database files from the INEGI digital elevation models used were E14C22A4_MS, E14C22B3_MS, E14C22D2_MS, and E14C22E1_MS (INEGI, 2015). Only points on the ground having a height smaller than 20 m over the mean sea level were considered in this study, as the highest wave height observed during historical tsunamis in Zihuatanejo was 11 m. To obtain the final bathymetry to be used by the model, combined information from Geomar bathymetry (Geomar, 2018) data and bathymetric charts of the regions near the Zihuatanejo Coast obtained from the Mexican Navy (Nautical Number chart 512.2 and 512.4) were used. Information from the ETOPO (Earth Topography) database was used for the deep-water regions (NOAA, 2015). The bathymetric and topographic data were homogenized and referred to the same datum, the mean sea level. The coastline features were also obtained from the bathymetric charts of the Mexican Navy.
Four grids with increasing resolution were used, at the open boundary between nested domains, a Flather (1976) radiation method for the barotropic flow was applied. The first grid, with a latitude and longitude resolution of 1/133 (825 m), is used for modeling the tsunami generation on the ocean floor; the second and third grids are transitional grids with latitude and longitude resolutions of 1/400 (275 m) and 1/1200 (91 m), respectively, which downscale the tsunami conditions from the origin grid to a local grid near the coast. The fourth grid has a latitude and longitude resolution of 1/3600 (30 m) and is used for flooding simulation in Zihuatanejo Town (Table 2; Figure 3).

| Initial conditions of the model
In this work, the initial sea surface displacements for hypothetical 7.0 and 8.0 Mw earthquakes along the Cocos Plate with epicenters at 16.9718 N-101.973 W, 86 km from Zihuatanejo (near the epicenter of the April 18, 2014 earthquake 7.2 Mw, the last hardest submarine earthquake near Zihuatanejo [Servicio Sismológico Nacional; SSN, 2020]), were calculated using the Okada (1985) formulae with the okada85 function (Beaudecel, 2019) for a fault plane with a depth, strike, dip, and rake of 20 km, 300 , 14 , and 100 , respectively ( Figures 5 and 6), values similar to those reported by the Harvard Centroid-Moment-Tensor Project catalog (CMT, 2020) for the September 19, 1985 earthquake 8.1 Mw. The fault length, width and average slip were calculated via the earthquake moment magnitude (Levin & Nosov, 2016), assuming incompressibility of water, and considering the initial water surface elevation from the mean sea level to be equivalent to the seabed displacement field induced by the slip on the fault plane. The movement of the seafloor was considered to be instantaneous, and the water surface followed the changes to the seafloor (Wang, 2009).
The simulated earthquakes have magnitudes of (a) 7.0 Mw, as this is similar to the magnitude of the earthquake that occurred on November 16, 1925, which caused a tsunami that reached the Zihuatanejo area with registered wave heights between 7 and 11 m (Sánchez & Farreras, 1993)  4 | RESULTS

| Tsunami propagation in the open sea
Time series of numerically simulated tsunami wave height were extracted at three virtual stations (Table 3) with different depths (Figure 4). As the tsunami wave propagates along the open ocean ( Figures 5 and 6), it faces the changes in the bathymetry of the region, which induce changes in the wave velocity and height. When the tsunami reaches shallow water, such as point C, the height increased from 0.2 to 0.4 m for the 7.0 Mw earthquake, and from 1.4 to 2.9 m for 8.0 Mw earthquake. Although these differences appear considerable, the behavior of the waves is similar for both cases; they just have different scales. The tsunami waves undergo deceleration and shrinkage, causing an increase in their height and overlap with the returning waves, as shown in Figures 5 and 6.   Zihuatanejo downtown, 3 located in the "Playa Principal" and 1 in "Playa la Madera" according to Figures 7 and 8; the touristic and densely populated area. Tsunami arrival time and wave height at Zihuatanejo beaches were extracted at these locations.

| Tsunami height and arrival time
The results of the arrival times at the virtual stations indicate that the average times taken for the first tsunamis wave to reach Zihuatanejo Bay are 19 and 18 min for the 7.0 and 8.0 Mw earthquakes, respectively. On the time series for the sea surface height at the coast for the tsunami generated by the 7.0 Mw earthquake (Figure 7), three significant sea surface peaks exceeding 1 m, at 20, 33, and 38 min were observed; and three significant troughs during the water drawdown, with a height of −1 m, at 24, 29, and 42 min. The third sea surface peak is the highest, reaching a height of 1.2 m, followed by the lowest water drawdown of −1.4 m. This indicates a possible case of resonance in the bay. Station 4, located on the La Madera beach at 1.7 m above mean sea level, is not affected by the tsunami waves induced by the hypothetical 7.0 Mw earthquake (Figure 7).
In the second case, the sea surface elevation exceeds 10 m, two peaks are between 5 and 6 m, and other 5 peaks reach 2 m (Figure 8). Constant negative values are shown as the sea surface height was low enough to reach the beach bed. The positive-phase peaks appear at 20, 33, and 38 min. The first peak is almost twice as high as the following two consecutive peaks due to the form and height of the initial conditions associated with the 8.0 Mw earthquake, reaching a wave height greater than 10 m; the third peak is slightly higher than the second. The sea surface height of Station 4 is with respect the ground elevation.
F I G U R E 6 Initial sea surface condition and propagation phase of the tsunami for hypothetical 8.0 Mw earthquake. Time elapsed from the earthquake In both cases, after 45 min from the beginning of the tsunami, the trend on the time series of the sea surface height shows that the wave height starts to decreases in amplitude, but its period change from 7 min to a period of 15 min.
A sequence of images that represents the flooding of the populated zones was created for the 7.0 and 8.0 Mw earthquakes using the results of the model, which show the most susceptible areas to flooding, the safest zones, and the routes followed by the flood. For the 7.0 Mw earthquake (Figure 9), the flood generated by the three peaks during the positive phase does not extend beyond the beach. Therefore, wave height of 1 m is insufficient to flood the city. However, for the 8.0 Mw earthquake (Figure 10), the first wave arrives 20 min and 30 s after the beginning of the tsunami with a height of 10 m, flooding the areas near the beach at Station 2 and extending 70 m inland, with an elevation of 10.8 m with respect to the mean sea level, and 3.8 m with respect to the ground elevation. Twenty-two minutes after the beginning of the tsunami, occurs the maximum extent of inland flooding of 400 m, near the municipal market, with an elevation of 6 m with respect to the mean sea level, and 0.9 m with respect to the ground elevation.
The use of MOHID allowed running nested models with a high resolution, which allow to have a greater detail to study the flood zones, and can also calculate the wave propagation in the first grid to model the impacts of the first wave that impacts the coast, as well as those that follow.
In the tsunami simulation for the hypothetical 7.0 and 8.0 Mw earthquakes with epicenters at 16.9718 N-101.973 W, on the Cocos Plate, the propagation and flood zone results were obtained. The initial water surface F I G U R E 9 Flood sequence produced by the tsunami generated by the 7.0 Mw earthquake in Zihuatanejo City. Time elapsed from the earthquake F I G U R E 1 0 Flood sequence produced by the tsunami generated by the 8.0 Mw earthquake in Zihuatanejo City. Time elapsed from the earthquake elevation was computed using the same fault parameters of strike, dip, rake, reported for the September 19, 1985 earthquake 8.1 Mw by the Harvard Centroid-Moment-Tensor Project (CMT, 2020), one of the most serious natural disaster to date in Mexico's history. The hypothetical magnitudes of the earthquake were chosen from 7.0 and 8.0 Mw, because at these magnitudes the highest wave height produced by a tsunami in Zihuatanejo has been reported. The earthquake of 7.0 Mw on November 16, 1925, caused that the sea level suddenly rose 6-7 m, invaded the streets and swept away the houses. This flood lasted between 10 and 15 min (Sánchez & Farreras, 1993). And the earthquake of 8.0 Mw on September 19, 1985, produce that the water first receded to the end of the 200 m pier, a depth of 2.5-3.0 m, then returned flooding the land 2.5-3.0 m above normal levels (Lander et al., 2003).
The arrival times of the tsunamis generated by the hypothetical 7.0 and 8.0 Mw earthquakes at Zihuatanejo are quite similar, ranging between 18 and 19 min. The main difference between the tsunamis is that the initial surface elevation associated with 8.0 Mw earthquake affect a larger area than the 7.0 Mw earthquake (Figures 5 and 6). According to the time series obtained at the stations located in the open sea (Figure 4), the waves of the tsunamis are long with a period of 10 min, and they present similar behaviors for both cases, but at different scales. The waves are also amplified as they reach the coast.
When comparing the results obtained with those of Farreras et al. (1999), it should be noted that they use the fault parameters of the 8.0 Mw earthquake of September 19, 1985, with the real location of the epicenter, while in this work, the epicenter changes place. They have a 50 km × 175 km fault size and a 90 cm slip, while we obtain a 43 km × 85 km slip and a 4.2 m slip. They find that the maximum computed water level is 2.4 m that reaches the 2 m land contour above the mean sea level, which coincides with the flood zones obtained for the case of our hypothetical 7.0 Mw earthquake, mainly the area near the beach.
The flood region obtained for the 8.0 Mw earthquake matches the area described by the numerical simulation of the tsunami of extreme events on November 16, 1925and June 22, 1932by Farreras et al. (1999 with a maximum computed water level of 8 m that reaches the 4 m land contour above the mean sea level.
Finally, the use of finer bathymetric data, as well as information from Digital Terrain Model helps to have a better and more realistic calculation of possible flood zones. In this case it is found that the maximum extent of the flood produced by the tsunami is 400 m inland, with an elevation of 6 m with respect to the mean sea level, and 0.9 m with respect to the ground elevation, affecting restaurants, hotels, banks, schools, businesses, museums, and even the municipal market of the city, which are all within a highly populated area.

| CONCLUSIONS
The study area has historically been reported as a tectonically active zone and, due to its location, it is affected by tsunamis; the wave heights recorded in 1925, 1985, 2003, 2014, 2015, and 2017 were between 7 and 11 m, 3, 0.6, 0.3, 0.6, and 0.4 m, respectively.
According to the time series of the arrival of the tsunami generated by a hypothetical 7.0Mw earthquake, the highest peak is not the first, as expected; rather, the fourth peak in both the positive and negative phases has a higher amplitude. Moreover, after the first wave train, which has the most important peaks, the wave height decreases while the wave period increases, indicating a possible case of resonance in the bay. Therefore, it is important to calculate the natural modes of oscillation in the bay.
In this study, possible flooding in Zihuatanejo Bay is only assessed for tsunami with local origin generated in the area of the Cocos Plate by a hypothetical 8.0 Mw earthquake with an epicenter at 16.9718 N-101.973 W. The tsunami wave takes 18 min to arrive. The maximum height at the coast of 10 m occurs after 20 min and 30 s, and the maximum flooding occurs after 22 min, extending 400 m inland and affecting restaurants, hotels, banks, schools, businesses, museums, and even the municipal market of the city, which are within areas with a high population.
The Zihuatanejo region is impacted on some occasions by hurricanes, distant waves produced by storms in the southern hemisphere and high astronomical tides. A scenario of maximum risk is to consider the effect of the swell of hurricanes, plus distant waves and spring tides, this type of conditions should be considered in future numerical experiments to improve the delimitation of safety zones for the population.

ACKNOWLEDGEMENTS
This work was supported by the Posgrado en Ciencias del Mar y Limnología of the Universidad Nacional Autónoma de México grant. This work was also supported by the CONACYT, México; with a master's degree scholarship under the PCT-UNAM (Posgrado en Ciencias de la Tierra of the Universidad Nacional Autónoma de México) granted to Víctor Kevin Contreras-Tereza and a master's degree scholarship under the PCML-UNAM (Posgrado en Ciencias del Mar y Limnología of the Universidad Nacional Autónoma de México) awarded to Rosalinda Monreal-Jiménez. Special acknowledgments are extended to INEGI, who provided the LiDAR information, to GEO-MAR, who provided the bathymetry information, and to Bently-MOHID Studio team, particularly to Luis Fernandes, who provided insight, knowledge, and expertise that enhanced this research. We are grateful to Jorge Castro for improving the study figures.