A first estimation of tsunami hazard of the Pacific coast of Costa Rica from local and distant seismogenic sources

Costa Rica has been affected by several local and distant tsunamis in the past, but the historical information is scarce and incomplete. Its Pacific coast stretches for over a thousand kilometers, and tsunami hazard has never been evaluated for its full extent. Numerical modeling of tsunami propagation and inundation is a useful tool to assess tsunami hazard, particularly in cases with limited historical information available. Here, we perform a first estimation of tsunami hazard for the Pacific coast of Costa Rica from seismogenic sources, by numerical propagation of 57 local and distant tsunamis to a depth of 20 m. The results of our study identified tsunami sources that are particularly threatening for Costa Rica and determined locations with higher tsunami hazard. For the analysis, the Pacific coast of Costa Rica was divided into segments and subsegments based on differences in continental slope morphology. Subsegments with higher tsunami heights were Southwest Nicoya Peninsula and West Osa Peninsula, and in a lesser extent North Guanacaste, North Nicoya Peninsula, and Central Pacific. Regions with long and gentle slopes and narrow continental shelf were affected by higher tsunami waves, due to more efficient tsunami energy transmission to the shelf and reduced energy loss while traveling through a narrow shelf. On the opposite, steeper continental slopes reflected most of the tsunami energy, causing smaller tsunami heights nearshore, regardless of the shelf width. Nevertheless, other effects played a major role, like curved coastlines that focused tsunami energy, wave refraction, interference, and trapped edge waves. Distant tsunamis dominated the threat, with tsunamis coming from the Tonga-Kermadec and the Colombia-Ecuador Trenches causing the greatest heights due to directivity, and arrival times of about 15 h and 75 min, respectively. Local tsunamis had short arrival times but a localized impact, mainly at the shoreline in front of the generation region but were also affected by tsunami focusing, wave refraction, and edge waves. Outer rise and Osa sources caused the lowest impact within local sources. These results provide a guide for emergency planners to prioritize coastal locations and tsunami sources for tsunami preparedness actions and warning protocols.


Introduction
Costa Rica is a major international destination for sun-andsand, ecological, and medical tourism (Matarrita-Cascante and Brennan 2010;Warf 2010;Van Noorloos 2011;Nost 2013). It has 1228 km of coastline, of which 1016 km is on the Pacific coast and about 200 km on the Caribbean coast (Arozarena Llopis et al. 2015).
The Middle America Trench (MAT) has narrow coupling zones, fast subduction rates, and well-defined segmentation (Protti et al. 1994;Protti et al. 1995), making the simultaneous rupture of several segments unlikely. Goff et al. (2014) defined mega-tsunamis as those with "initial wave height/ amplitude at source exceeding 100 m / 50 m respectively". Such mega tsunamis were not regarded as likely to originate in Central America in the conclusions of the Experts Meeting Tsunami Hazard in Central America: Historical Events and Potential Sources (UNESCO/IOC 2018).
However, Costa Rica experiences frequent earthquakes, some of which provoke tsunamis, as a consequence of the interaction of three tectonic plates: the Cocos Plate, the Responsible Editor: Christoph Voelker Caribbean Plate, and the Panama Block. From the 39 tsunamis recorded in Costa Rica since 1746, 17 were caused by local earthquakes at either shores, one of which led to at least three deaths on the Caribbean coast (Chacón- Barrantes et al. 2021). There have been no deaths associated to tsunamis on the Pacific coast, but runups up to 7.3 m have been observed, with the highest occurring during the 1992 Nicaragua Mw = 7.7 and the 1950 Nicoya Mw = 7.8 tsunamis.
Some tsunamis can cause damage far away from their source, as happened in Japan and  This paper presents a review of tsunami hazard for the Costa Rican Pacific coast due to local and distant earthquake-generated tsunamis. Analysis of the hazard derived from landslide generated tsunamis is beyond the scope of this paper but will be considered in a future study. For example, von Huene et al. (2004) estimated maximum wave height of 27 m above the prehistoric Nicoya Slump (Fig. 2b), in the case it slid all at once. Although this maximum wave height was estimated by equations rather than numerical modeling, the estimated area of the Nicoya Slump is about 1350 km 2 (von Huene et al. 2004), which would imply a considerably high tsunami.
Along the Pacific coast of Costa Rica, there are 685 beaches and 273 towns or cities that are exposed to tsunamis, both from near and far field sources. The goal of this work is to get a first estimation of locations with highest earthquakegenerated tsunami hazard along the Pacific Coast of Costa Rica. This will be done by a deterministic approach of aggregated scenarios like the one employed by Álvarez-Gómez et al. (2013) or Tinti and Armigliato (2003). For this first approach, the distant tsunami sources employed are a combination of pre-calculated unit sources available in the numerical model. Seismic sources at the Middle America Trench (MAT) existing in literature are used to model local tsunamis. It will be analyzed the interaction of tsunamis with the continental slope and shelf offshore each location, in terms of the segmentation defined by Harders et al. (2011). The definition of coastal locations with highest hazard will allow local authorities to prioritize tsunami preparedness.

Tsunami interaction with continental platform and shelf
There are general relationships that can be established between continental platform, continental slope angle and width, and tsunami shoaling. Levin and Nosov (2009) defined the viscous linear dissipation length as "the distance, over which the amplitude of a long wave in a viscous liquid becomes e times smaller", which is a function of wave frequency and 3/2 power of depth. They also defined the non-linear dissipation length as "the distance along which the wave amplitude becomes two times smaller" (Levin and Nosov 2009). They established this length as a function of depth square and wave amplitude. Once the tsunami reaches the continental Fig. 1 Historical tsunamis registered on the Pacific Coast of Costa Rica (squares and triangles) (NOAA/NCEI 2021) and hypothetical distant tsunami scenarios used here. The COMMIT propagation grid at the Pacific Ocean is shown with a magenta rectangle slope, it behaves non-linearly. Thus, continental shelves are more efficient on tsunami dissipation than abyssal plains, not only because of shallower depths increasing the dissipation length, but also because of the higher power dependency on depth for non-linear processes. Consequently, tsunamis traveling over wider continental platforms will experience higher dissipation. Siva and Behera (2016) performed numerical onedimensional experiments of N-waves shoaling over different values of continental slope and depth of continental shelf. They found that for continental slopes of less than 11°, a decrease in slope drives to an increase on wave amplification. The authors mentioned that "as the continental slope becomes flatter, more energy is transmitted to the shelf and higher tsunami height is observed". For the two-dimensional problem the relationship between continental slope angle and tsunami amplification might not be straightforward.
In the following Section, we describe the geomorphology of the study area, its segments and subsegments, which will be used in the analysis of the tsunami response of each of them discussed in Section 5.4.

Study area
The Middle America Trench (MAT) extends from Mexico to Costa Rica, where the Cocos Plate subducts beneath the Caribbean Plate. Six different tectonic-style segments are defined for the entire MAT according to Harders et al. (2011). Those segments were defined based on different Cocos plate morphologies and continental slope angles and widths (Table 1), and are not the same as the segments defined by Protti et al. (1994Protti et al. ( , 1995 based on tectonic coupling. The later was employed to define worst-case scenario seismic sources by a group of experts that met in 2016 under the leadership of IOC/UNESCO (UNESCO/IOC 2018). These seismic sources are listed in Section 3.3.2.
The first three tectonic-style segments are located offshore Guatemala, El Salvador, and Nicaragua. For Costa Rica, there are three tectonic-style segments of MAT based on the geomorphology of the continental slope (Table 1 and  This segmentation of the MAT, and thus of the coast, is based on the evidence of massive landslides along the continental slope. The presence of those landslide scars and their characteristics are directly related to subducting floor roughness (von Huene et al. 2004). They are defined as tectonic erosion processes and dominate the convergent margin of the Middle America Trench (Harders et al. 2011). Evidence for tectonic erosion in the Neocene has been found in the Central Costa Rica segment (Ranero and von Huene 2000;Vannucchi et al. 2003).  Here, we subdivide those three segments in terms of their geomorphology, as detailed in Table 2.
Nicoya Gulf is not in front of the continental slope neither oriented to the trench. Nicoya Gulf mouth is 16 km wide, and the gulf is 60 km long. About 14 km inside the Gulf is located the San Lucas Island at the western margin and Puntarenas city in a sandbar at the eastern margin of the Gulf, the distance between them is only 5 km. The maximum depth at Nicoya Gulf is about 30 m, and its average depth is about 10 m. In this work, the tsunami heights will be analyzed at 20-m depth, due to the absence of high-resolution coastal bathymetry along most of the coastline. Then, the hazard assessment for the Nicoya Gulf would require a different approach including the modeling with high-resolution bathymetry, which is not currently available. Therefore, Nicoya Gulf was not included in this study, but the authors strongly recommend the evaluation of its tsunami hazard particularly for its islands and cities, like Puntarenas.

North Nicoya segment
The North Nicoya segment has the softest subducting oceanic floor in the entire Middle America Trench (Harders et al. 2011). As a result, its continental slope is mainly straight and steep and defined by mid and lower slope slides and slumps, with an average 4-6°dip and a 53-km-wide (low variance) slope (Harders et al. 2011) (Table 2 North Guanacaste subsegment).
It can be divided into three subsegments in terms of its geomorphology. The first is North Guanacaste (1.1, yellow in Fig. 2b), characterized by a highly indented coast, with peninsulas and bays of varying sizes showing no main orientation. The continental shelf in this area is wide, decreasing from north to south, and ranges from 100 km in the northernmost end to 35 km in the south. In front of this subsegment, the Cocos Plate presents a smooth subducting surface, which is referred to as the Smooth Domain (Denyer and Kussmaul 2012).
The second and smallest subsegment is the North Nicoya Peninsula (1.2, light purple in Fig. 2b). It is also very indented, but only at the meso-scale, and is generally clearly north-west oriented. The continental shelf becomes narrower in this area, also decreasing from north to south, ranging from 35 km in the northern part of the segment to 20 km south of Tamarindo bay.
The third subsegment is the West Nicoya Peninsula (1.3, light green in Fig. 2b), its coastline is the straightest of all segments analyzed. The continental shelf in this area is extremely narrow, with widths around 8-10 km, and the continental slope is the steepest of the entire North Nicoya  Segment. Coastal uplift in this area has been documented by Fisher et al. (1998), Marshall et al. (2001), and Fisher et al. (2004). Strong earthquakes (Mw > 7.5) in the Nicoya Peninsula occur approximately every 50 years (Protti et al. 2001). During the interseismic period, the Caribbean and Cocos plates are coupled, and the failure occurs when the coupling is overcome (Protti et al. 2010). During the approximately 50 years of stress accumulation, the coast of the Nicoya peninsula undergoes subsidence, and when an earthquake occurs a general coastal uplift is observed (Protti et al. 2010). As a consequence, the Peninsula has several marine terraces and raised platforms (Marshall 2007) indented with large bays like Papagayo and Tamarindo (T in Fig. 10) and pocket beaches like Playa Garza (Gr in Fig. 10) and Playa Samara (Sa in Fig.  10). Uplifted beaches are very common in this segment, especially in the central area between playa Tamarindo and Playa Ostional (O-N in Fig. 10).

Seamounts-Central segment
The Seamounts-Central segment extends from Punta Guiones (PG) to Playa Tortuga (To, on the right margin of the mouth of the Térraba river, Fig. 2b), excluding the Gulf of Nicoya. Offshore, and especially from the southernmost part of the Nicoya Peninsula to the Quepos promontory (Qp in Fig. 2), the Cocos plate becomes rough through the rough-smooth limit in this segment (Lonsdale and Klitgord 1978;Meschede and Frisch 1998;Denyer and Kussmaul 2012). This segment has rotational landslides in its upper slope (their scars in the shelf break are shown with dashed lines in Fig. 2), whereas rock avalanches and slumps are present in the mid-and lower slopes (Harders et al. 2011). These landslides and slumps were produced by subducting seamounts and decreased the steepness of the continental slope compared to the general trend in the Guanacaste area (von Huene et al. 2004). Its average dip slope and width are 3.5°and 56 km, but the width varies from 44 to 73 km (Harders et al. 2011).
The Seamounts-Central segment can be divided into three zones: the Southwest Nicoya Peninsula (2.1, orange in Fig.  2b), the Southeast Nicoya Peninsula (2.2 brown in Fig. 2b), and the Central Pacific (2.3, fuchsia in Fig. 2b). Since the southern part of the Nicoya Peninsula is located over the subducting seamounts of the "rough domain" (Marshall 2007), it was not included with the rest of the Peninsula in the first segment. The different Quaternary uplift rates between the North and South Nicoya peninsula may be linked to differences in subducting plate roughness, dip, and thickness (Marshall 2007).
The Southwest Nicoya Peninsula is morphologically heterogeneous because of the presence of two geological patterns: one formed by sedimentary rocks, and the other characterized by oceanic basalt outcrops (Denyer and Alvarado-Induni 2007). There is a large landslide scar (von Huene et al. 2004) in the shelf break in front of this subsegment (dashed line in Fig. 2b) which lowers its steepness greatly, especially when compared with the North Nicoya segment (von Huene et al. 2004;Harders et al. 2011). Coastal morphology in this area reflects the semicircle of the landslide detachment area (Fig. 2b). The northern part of this subsegment is WNW-ESE oriented, with scattered capes and bays, and surrounded by steep terrains, even in some major bays (e.g., Samara Bay, Sa in Fig. 10). In general terms, the southern part is open, straight, and NNW-SSE oriented. It is characterized by long sandy beaches and wide coastal plains (B: Bejuco, SM: San Miguel de Jabilla, Cy: Coyote, Ca: Caletas, Pe: Pencal, M: Manzanillo, all in Fig. 10). Continental shelf width in this subsegment ranges from 8 to 10 km in the northernmost area to as much as 17 km in the south.
The Southeast Nicoya Peninsula subsegment (2.2, brown in Fig. 2b) is in the province of Puntarenas, which extends southward along the rest of the Costa Rican Pacific Coast. It is SW-NE oriented, as opposed to the general NW-SE trend of the Costa Rican Pacific Coast, alternating sedimentary and oceanic volcanic rock outcrops (Denyer and Alvarado-Induni 2007) which form very steep coastal slopes. Given the orientation of this segment, mostly perpendicular to the trench, continental shelf width is variable, increasing from south to north (from the outermost part of the peninsula, to the innermost part), ranging from 8 to 10 km in the south, to almost 40 km for the northernmost (innermost) part of this subsegment. Continental slope in this subsegment is also characterized by the presence of a massive landslide, which reduces considerably the average angle of dip (dashed line in Fig. 2b).
The Central Pacific subsegment (2.3 fuchsia in Fig. 2b) shows different combinations of continental slope and shelf morphologies. The continental slope in front of this subsegment is affected by two massive landslides offshore Herradura and Quepos (dashed lines in Fig. 2b) that decrease the shelf width in those places. One of these landslides also affected the previous subsegment (Southeast Nicoya Peninsula). The eastern part of this subsegment (close to Playa Tortuga, To in Fig.  2b) is characterized by a wide shelf, no landslides, and a very short and abrupt continental slope.
The Seamounts-Central segment is also geologically and geomorphologically different from the Nicoya and Cocos Ridge-Osa Peninsula regions. While uplifting in the Nicoya and Osa Peninsula is infrequent but involves large displacements, uplifting between those peninsulas involves small but frequent earthquakes (Von Huene et al. 1995). The coast is characterized by differential Quaternary uplift rates due to faults perpendicular to the coastline that divide the coast into several blocks with fluvial and marine terraces (Fisher et al. 1994;Fisher et al. 1998;Marshall 2007). As a result, there are relatively extensive beach systems broken by only two remarkable rocky promontories (Herradura and Quepos), where several small pocket beaches can be found. Wide coastal plains flank the shores between Esterillos (WE in Fig. 10) and Quepos (Fig. 2b), while steep slopes come closer to the coast from playa Savegre (Sv in Fig. 2) to the south end of the segment. Except for the two promontories, no notable rocky outcrops are present in the coastal segment to the south of Playa Hermosa and the north of Playa Dominical (Do in Fig.  2). This subsegment extends NW-SE, with very few exceptions such as the Quepos and Herradura (Fig. 2) promontories.

Coco's Ridge-Osa Peninsula
This segment is associated with different continental slope characteristics which have led to the division in four subsegments. The first one, the West Osa Peninsula (3.1, dark red in Fig. 2b), is associated with a relatively wide continental shelf, and the absence of major landslides on the continental slope. It represents a transition area between the Central Costa Rica landslides, and the extremely narrow continental slope in front of the Osa Peninsula, which is also characterized by landslides (Fig. 2b). The width of the continental shelf decreases from NW to SE, reaching a width of 4 km near to Punta Chancha (Ch in Fig. 10), where the next subsegment begins.
The South Osa Peninsula (3.2, dark blue in Fig. 2b) is characterized by the narrowest continental shelf of the whole Costa Rican Pacific coast, where it hardly reaches 2 km wide in some points (near Matapalo Cape, * in Fig. 10). This narrowness is due in part to the existence of slumps and rock avalanches in the lower-middle slope, mainly due to the subduction of crests owing to the Coco's Ridge (Harders et al. 2011), which have made the shelf break to retreat backwards. The continental slope of this segment is one of the steepest of the entire study area.
Dulce Gulf mouth is 14 km wide, and the gulf is 54 km long divided in 23 and 31-km sections oriented south-north and southeast-northwest, respectively; it is characterized by high depths, reaching 100 m in some points. For the Dulce Gulf subsegment (3.3, dark purple in Fig. 2b), the continental slope is also affected by massive slumps and avalanches, as the previous subsegment. The slope is narrow and steep with an average width of 8 km and angles of 6-20° (Harders et al. 2011). The continental shelf is practically inexistent, given that the shelf break is almost placed at the gulf mouth.
The last subsegment, Punta Burica (3.4, cyan in Fig. 2b), is associated with a very narrow continental shelf (almost inexistent) together with steeper continental slope, where no major landslides have been described. Its coastline is mostly straight with exception of a cape and a bay south of Puesto La Playa (PLP in Fig. 10) This subsegment ends in correspondence to the Panama Fracture Zone, PFZ in Fig. 2 (Denyer et al. 2009).

Numerical model
ComMIT (PMEL/NOAA 2016) model was used, which is a graphical user interface (GUI) for MOST numerical model (Titov and González 1997). It performs offline nesting between three user-defined inundation grids (A, B, and C), and its incorporated 4 arc-minute propagation grid covering the entire Pacific Ocean. ComMIT employs linear shallow water equations in the propagation grid, and non-linear shallow water equations in the inundation grids, which are nested online between them. The seismic source can be defined by a linear combination of pre-calculated unit seismic sources or by the simulation of user-defined seismic sources.

Model setup
The model setup performed by the authors consisted of three online nested grids of 60, 12, and 4 arcsec resolutions (Grid A, Grid B, and Grid C respectively in Fig. 2a). Grid C covers Quepos bay, which tide gauge has recorded 15 tsunamis (TG in Fig. 2a). The records of the following tsunamis were used for model verification in Chacón-Barrantes and Gutiérrez-Echeverría (2017) The bathymetry of the grids was obtained combining GEBCO-08 data up to a depth of 200 m and digitized nautical charts for depths less than 200 m. The lack of high-resolution and updated coastal bathymetry is a major handicap for Costa Rica; nonetheless, nautical charts are reliable enough to permit characterization of hazard nearshore. Local effects, like the response of small bays to tsunamis, cannot be analyzed with this resolution; however, the model allows the inundation of dry land in all grids.
The topography was constructed based on GEBCO-08 data. While coastal LIDAR topographic data up to 1-m resolution are available, the low resolution of the bathymetric data available did not justify the use of this data. The LIDAR topographic data will be used in future work to set up highresolution inundation models for selected locations where higher-resolution bathymetric data will be measured directly.
The simulations were carried on for 12 h after tsunamis' first arrivals to account for shelf resonance and trapped waves. The maximum tsunami height was chosen between all the seismic scenarios modeled at each node of grid B. These joint maximum tsunami heights were used to analyze tsunami interaction with the continental slope and platform; only grid nodes with depths larger than 20 m were considered due to the lack of reliable high-resolution coastal bathymetry. Then, for each beach, a maximum tsunami height value was assigned from the closest nodes.

Tsunami sources
In this study, we used 57 worst-case scenarios, aggregating them to obtain maximum nearshore tsunami heights like performed by Álvarez-Gómez et al. (2013) or Tinti and Armigliato (2003). Worst-case scenarios consist of simplified tsunami sources, characterized by one or more rupture planes with homogeneous slip, instantaneous and simultaneous ruptures. Probabilistic or hybrid probabilistic-deterministic approaches such as Brizuela et al. (2014) and Zamora and Babeyko (2020) were not considered here but might be employed in future work.

Distant tsunamis
The Pacific Ocean has experienced more tsunamis than any other ocean basin throughout the period for which data is available. Around this basin, five regions have generated earthquakes with Mw≥ 9.0, four of which were recorded by gauges in Costa Rica: 1952Kamchatka, 1960Chile, 1964 Alaska, and 2011 Japan (gray triangles in Fig. 1) (NOAA/NCEI 2021). There are not records of effects of the Cascadia tsunami (1700) in Costa Rica, probably due to energy directivity.
Therefore, 36 hypothetical earthquakes were used as tsunami sources from all around the Pacific Basin with a standard magnitude of Mw = 9.3 (red and blue fault planes in Fig. 1), identified as the worst-case scenario, in a similar way as done by Borrero et al. (2015). These hypothetical earthquakes were originated at the following trenches: Mx: Middle America (offshore Mexico), JF: Juan de Fuca (Cascadia), GA: Gulf of Alaska, Al: Aleutian, Ku: Kuril, Jp: Japan, Nk: Nankai, Bo: Bonin, Ma: Marianas, NH: New Hebrids, TK: Tonga-Kermadec, Ch: Chile, PC: Peru-Chile, and SA: South American. The initial tsunami conditions for all hypothetical scenarios were obtained by linear combinations of precomputed unit sources, with homogeneous slip distribution.

Local tsunamis
For local tsunami sources were considered only sources existing in literature, as the definition of sources was beyond the scope of this paper. In June 2016, UNESCO/IOC organized the Experts Meeting "Tsunami Hazards in Central America: Historical Events and Potential Sources" (UNESCO/IOC 2018), with the objective "Identification of credible sources of tsunamis that could significantly impact the Pacific and Caribbean coasts of Central America and that can be used for tsunami modelling, evacuation mapping, planning and exercises". As an outcome of this meeting, four tsunami sources were identified along the Middle America Trench (Fig. 3b) which were considered worst-case scenarios for tsunami sources along the Pacific coast of Central America by over 20 national and international experts attending the meeting. The seismic parameters of the sources are listed in Table 3. More detailed information on the definition of the sources can be found in the Experts Meeting's Report (UNESCO/IOC 2018). Although the use of such simplified sources might seem unrealistic, they were selected to simplify the modeling process and the hazard assessment.
Also, the 16 largest tsunami sources considered in Zamora and Babeyko (2015) with Mw ≥ 7.3 and the worst-case scenario given by Álvarez-Gómez et al. (2012) with Mw = 7.9 were simulated here (Table 3 and Fig. 3a) for a total of 21 local sources. Some of the authors of the mentioned studies participated in the Experts Meeting or were consulted for the definition of those sources (UNESCO/IOC 2018).

Results and discussion
The values of the maximum tsunami heights given here are used only for comparative purposes, as they were calculated at a depth of 20 m. Coastal tsunami heights were not estimated using Green's Law, since it is not valid for all the locations along the studied coastline. Each coastal location would respond differently to the tsunami wave, depending on local bathymetry and geomorphology, which effects were not considered here.

Distant tsunami scenarios and directivity
Tsunamis from two regions posed the highest hazard for Costa Rica due to the direction of their maximum energy: the Tonga-Kermadec Trench and South America Trench (Figs. S8, S9, and S11). The maximum nearshore tsunami heights for Costa Rica were of 12.5 and 12.3 m and were caused respectively by one scenario offshore Colombia (SA2) and one scenario at Kermadec Trench (TK4) about 1000 km north-east from New Zealand. These scenarios are shown in Fig. 4a,b together with deep ocean amplitude plots (insets) showing the directions of maximum energy. Other sources at Tonga-Kermadec trench, TK1, TK2, TK3, TK5, and TK6, caused maximum offshore tsunami heights of 4.0, 3.1, 4.6, 1.3, and 2.9 m respectively (Figs. S8 and S9). The other source at the South American Trench (SA1) caused a maximum offshore tsunami height of 3.3 m (Fig. S11).
These results agreed with historical records of three tsunamis originated at the Colombia-Ecuador Trench and three tsunamis originated at the Tonga-Kermadec Trench that were recorded in Costa Rica; particularly, the 1906 Colombia tsunami was widely witnessed in the country. On the opposite, large tsunamis originated in Chile (1960, 2010, 2014, and 2015), Japan (2011), Kamchatka (1952 and2007), and the Aleutian Arc (1964 and 1957) did not flood Costa Rican shores, except for Cocos Island. Nonetheless, some of these tsunamis caused important coastal currents that resulted in damaged ships and coastal erosion.
The Colombia scenario implied the shortest arrival time for distant tsunamis, of about 1 h and 15 min, posing the highest hazard of all the scenarios considered here. The tsunami came from the south and was refracted to the Central Pacific and Nicoya Peninsula (Fig. 5). Historically, three tsunamis from that region have arrived at Costa Rica: 1906Colombia Mw = 8.8, 1979, and 2016 Ecuador Mw = 7.8.  For one Kermadec scenario, the strongest energy beam was focused directly on Costa Rica, arriving from the southwest and first hitting the Northwest Nicoya Peninsula subsegment. The wave was then refracted to the Southwest Nicoya   (Fig. 6) and the rest of the coastline. The long travel time of more than 15 h would allow timely evacuation but still would pose an important hazard. There are three records of tsunamis originated at the Tonga-Kermadec subduction zone that arrived in Costa Rica: 7.9 cm registered in Quepos after a Mw = 8.0 earthquake in Tonga in 2006 and sea level perturbations at Wafer Bay, Coco's Island, after the 2009 Samoa earthquake, Mw = 8.1. While this manuscript was under review, a tsunami was originated at Kermadec Islands on 4 March 2021 after a Mw = 8.1 earthquake; the tsunami was registered at the tide gauges in Chattam Bay (Coco's Island) and Quepos with maximum wave heights of 5.8 and 17.6 cm, respectively (Chacón- Barrantes et al. 2021). It was also witnessed as sea level perturbations and strong currents at Wafer Bay, Coco's Island, lasting for about 3 h with wave periods of about 7 min; a video was posted in YouTube at https://youtu.be/cLqGJIE1glo (in Spanish). The modeling of this tsunami also reproduced correctly the tsunami records at Quepos and Coco's Island tide gauges, and a separated manuscript is being prepared with these results and records.
Historically, there are no records of mega-earthquakes (Mw ≥ 9.0) generated at those subduction zones. IOC/ UNESCO organized experts' meetings in 2018 and 2019, to define tsunami worst-case scenarios originated at the Tonga-Kermadec and Colombia-Ecuador subduction zones, respectively. The reports of both meetings were published while this manuscript was on its final phase, and therefore it was not possible to include those sources in this work. The experts considered that the Colombia-Ecuador subduction zone can generate earthquakes up to Mw = 8.7 with high probability, and Mw = 8.9 with an intermediate probability (IOC/ UNESCO 2021). The experts considered that Tonga-Kermadec subduction zone could generate earthquakes up to Mw = 9.3 with high probability, and up to Mw = 9.7 with intermediate probability (IOC/UNESCO 2020). The sources identified at both experts' meetings will be modeled in future work, as some of the sources identified for Tonga-Kermadec seem to direct their energy towards Costa Rica: high probability source C (Mw = 8.9), intermediate probability sources ABC (Mw = 9.3), ABCD (Mw = 9.5), ABCDE (Mw = 9.6), BCD (Mw = 9.4), CD (Mw = 9.2), CDE (Mw = 9.3), and HABCDE (Mw = 9.7) (IOC/UNESCO 2020). Emergency planners should consider tsunamis originated at both regions as special cases in evacuation plans and tsunami warning procedures.
In the supplementary material are shown the maximum tsunami heights nearshore (left panels) and in deep ocean (right panels) for all the distant sources simulated. Also, there is a table with maximum heights for each source. Sources from México (Fig. S1) caused maximum offshore tsunami heights up to 3.3 m, Juan de Fuca and Gulf of Alaska sources (Fig. S2) up to 0.46 m, Aleutian Arc (Fig. S3), Kuril (Fig. S4), Nankai, and Bonin (  (Fig. S10) and Perú-Chile (Fig. S11) sources up to 3.94 m. All those sources do not pose a large hazard for Costa Rica mainland due to their directivity. However, Coco's Island is more exposed to tsunamis from those sources, as was seen during the 2010 Chile and 2011 Japan tsunamis (Chacón-Barrantes and Gutiérrez-Echeverría 2017). Fortunately, vulnerability is much lower there due to its low population and high relief; thus, it was not included in the study area.

Local tsunami scenarios
Maximum tsunami heights in meters for six of the 21 local scenarios are shown in Fig. 7, together with their rupture areas, when possible. Modeling results from all the 21 local sources can be seen in Figs. S12 to S15 and Table S2.
Local tsunamis produced much lower nearshore heights than distant tsunamis, because of the moderate earthquake potential of the Middle America Trench. For all local scenarios simulated, the impact was localized and not a single scenario caused tsunami heights over 1 m for the entire coastline. Only four local scenarios had maximum offshore heights over 3 m: GUANICA and NICOBANO from the 2016 Central America Experts Meeting (UNESCO/IOC 2018), Fig. 7a,b respectively, and PM2 and C2 from Zamora and Babeyko (2015), Fig. 7c,e respectively.
The maximum height was 3.48 m corresponding to NICOBANO scenario, at the Southwest Nicoya Peninsula subsegment that is located within the deformation region, also due to energy focusing over the concave shape of the continental slope. As the edge of the deeper rupture plane is aligned with Punta Mala (PM in Fig. 7), the tsunami energy focused there, causing also higher tsunami heights.
Despite its oblique incidence, the tsunami caused by the GUANICA scenario got refracted by the concave shape of the shelf break (Fig. 8a), focusing the energy on the Southwest Nicoya Peninsula subsegment as well. This scenario also caused wave refraction and interference at Papagayo Gulf (Fig. 8b) and Central Pacific subsegment (Fig. 8c). However, tsunami heights were much smaller south of Punta Chancha.
C2 scenario caused the higher tsunami heights in the Central Pacific subsegment as it was the closest to the generation region. Moreover, the location of this fault plane mostly over the shelf encourages trapped waves, which will be discussed in Section 5.4.2.
Seismic faults with uniform slip underestimate tsunami heights particularly for local tsunamis (Geist 2002;Ruiz et al. 2015). However, it was beyond the scope of this paper the definition of local seismic sources. As mentioned before, the purpose was to get a first approach of the tsunami threat from already identified sources in the literature. As local tsunamis did not dominate the hazard for Costa Rica, this approach is sufficient at this stage. In the future, the authors might consider the use of several methods of slip distribution,

Regions with higher joint maximum tsunami heights for distant and local tsunamis
To determine the joint maximum tsunami height, it was chosen the maximum value among all the scenarios simulated here, for each cell in grid B with depths larger than 20 m. The results for distant and local tsunamis are shown in Fig.  9a,b, respectively. The joint maximum tsunami height for all scenarios is shown in Fig. 10.
Two subsegments presented the highest joint maximum tsunami heights within the Pacific coast of Costa Rica: Southwest Nicoya Peninsula and West Osa Peninsula (Fig.  10). Those subsegments also presented higher tsunami heights both for distant and local tsunamis (Fig. 9).
For distant tsunamis, most of the coastline had nearshore maximum tsunami heights of less than 10 m, except for Southwest Nicoya Peninsula and West Osa Peninsula that had maximum tsunami heights up to 12.5 m (Fig. 9a). Other locations with tsunami heights over 7.5 m were Salinas Bay, Papagayo Gulf, and close to the mouth of Nicoya Gulf.
Nearshore maximum tsunami heights for local tsunamis were lower than 2.5 m for most of the coastline. Nearshore tsunami heights up to 3.5 m were obtained at Salinas Bay and Papagayo Gulf, Southwest Nicoya Peninsula, close to Nicoya Gulf mouth and West Osa Peninsula (Fig. 9b). The moderate maximum tsunami heights for local scenarios are good news considering the arrival times of between 10 and 30 min.
The hazard specific to each segment considering all the 57 modeled scenarios will be analyzed in the following section.

Hazard analysis by region
The beaches were classified in five groups based on the joint maximum nearshore tsunami height offshore of each beach. The ranges of these groups were defined arbitrarily to facilitate the prioritization of the community preparedness work at the communities nearby or at those beaches.
Beaches on the lower groups are not exempt from tsunami hazard. According with IOC/UNESCO (2014), tsunami heights over 0.30 m represent a marine threat, tsunami heights over 1 m coastal flooding threat, and tsunami heights over 3 m are a major tsunami threat.
This study classifies beaches and not coastal towns and cities as it does not consider vulnerability. Also, beaches provide a more static framework than towns and cities, that might expand, be relocated or established in the future.
The groups were defined as: & Group 5, purple, maximum nearshore tsunami height larger than 10 m, & Group 4, red, maximum nearshore tsunami height between 7.5 and 10 m, & Group 3, orange, maximum nearshore tsunami height between 5 and 7.5 m, & Group 2, yellow, maximum nearshore tsunami height between 2.5 and 5 m, & Group 1, green, maximum nearshore tsunami height less than 2.5 m.
The beaches are plotted in Fig. 10 with the color of their respective group and maximum nearshore tsunami height in meters. Their vulnerability still needs to be assessed to define their tsunami risk.
Due to the low resolution of the model setup, local effects were not fully considered. Therefore, higher resolution inundation modeling is strongly recommended at locations showing larger near-shore tsunami heights to assess their geophysical responses and tsunami hazard levels.

North Nicoya segment
The region with smaller maximum tsunami heights was the northern part of North Guanacaste subsegment (1.1), mostly because Santa Elena Peninsula (SEP in Fig. 10) shadowed the tsunami energy. This, together with the wider continental shelf, caused smaller heights at Santa Elena Gulf (SEG in Fig. 10), where all beaches were classified in group 2. However, at the entrance of Salinas Bay (BS in Fig. 10), nearshore tsunami heights were larger than 5 m due to reflection at the Nicaraguan coast; then, all beaches within this bay were classified in group 3.
On the other hand, wave interference caused larger tsunami heights at Papagayo Gulf (PG in Fig. 10 and Fig. 5d), at the southern part of North Guanacaste subsegment (1.1). Due to refraction on the continental slope, tsunami waves reached the end of Santa Elena Peninsula and Cabo Zapotal at the same time (SEP and CZ in Fig. 5d). Then, both sides of the wave traveling along the coast from the north and the south arrived at the same time to Naranjo (N in Fig. 10) beach and interfered positively with each other, causing higher tsunami heights. Naranjo and other few beaches belonging to Santa Elena National Park were classified in group 4; fortunately, most of these beaches are not open to the public. The rest of the beaches in the southern side of Santa Elena Peninsula were classified in group 3.
The smaller shelf width and irregular coastline of the southern part of North Guanacaste subsegment and of North Nicoya Peninsula subsegment contributed to larger tsunami heights at both places, due to interactions of the tsunami with the coastal bathymetry. At the entrance of Culebra Bay, at the southern part of North Guanacaste subsegment (Fig. 10), nearshore tsunami heights were more than 8 m; then, all beaches within this bay were classified in group 4.
All the beaches at the subsegment 1.2, the North Nicoya Peninsula, were categorized as being in hazard groups 4 and 3, except for Conejo (Co in Fig. 10) on group 2, which currently can only be accessed by boat. Beaches in hazard group 4 are within Potrero Bay and Grande and Tamarindo. Potrero is a 2km-wide pocket beach located in a bay where several other pocket beaches are present and has registered several tsunamis, which might suggest additional local effects not considered here. Landward, the terrain is characterized by low steepness, especially in its northern part. Beaches in the neighbor Brasilito, Grande, and Tamarindo bays (T in Fig. 10) were in group 3.
The beaches in the West Nicoya Peninsula subsegment (1.3) were classified mostly in group 3. This subsegment had moderated tsunami heights; on one hand, these heights are higher than the northern part of North Guanacaste, as it faces the trench and then is more exposed to tsunamis than the latter. On the other hand, this subsegment had smaller heights than the southern part of North Guanacaste and North Nicoya Peninsula. This subsegment has the straightest coastline and the steepest continental slope of the segment; both factors contributed to the smaller tsunami heights compared to the forementioned subsegments, agreeing with theoretical predictions detailed in Section 2.
Beaches with higher heights at the West Nicoya Peninsula subsegment were Ostional and Nosara (O-N in Fig. 10), in group 4. These beaches are open and sandy shores near to main river mouths and have smoother underwater landscapes, which might amplify the tsunami even more once it approaches to the shore.

Seamounts-Central segment
The Southwest Nicoya Peninsula (2.1) together with the West Osa Peninsula (3.1) subsegments had the highest concentration of beaches assigned to hazard groups 5 and 4. For the Southwest Nicoya Peninsula, this is due to the combination of two factors: the wider and less steep continental slope that increase the tsunami energy transmitted to the shelf (Siva and Behera 2016), together with the concave shape of the continental slope that focuses tsunami energy to beaches in the central part of the curved slope (Saloor et al. 2020). Beaches in group 5 in this subsegment are shown in Fig. 10: Camaronal (Cm), Punta Islita (PI), Corozalito (Cr), Bejuco (B), San Miguel (SM), Coyote (Cy), Caletas (Ca), and Pencal (Pe). Beaches in group 3 are located between Punta Guiones (PG in Fig. 2b) and Carrillo (C), and between Santa Teresa (M) and Cabo Blanco (CB in Fig. 2b).
The Central Pacific subsegment (2.3) had beaches in groups 4 and 3 and only a couple of very small non-visited beaches in group 2. Within beaches in group 4, the most visited are Jacó, Hermosa, West Esterillos (WE in Fig. 10) and the beaches of the Ballena National Park (BNP).
Tsunami heights caused by the Colombia SA2 scenario were extracted along the 25-m isobath offshore subsegments 2.2, 2.3, and 3.1, and their time evolution is shown in Fig. 11b. Most of the tsunamis simulated here had similar behavior.
Tsunamis reached both Cabo Blanco (CB) and Punta Llorona (PL) almost at the same time due to refraction (dashed horizontal line in Fig. 11b). This section of the continental platform between Nicoya and Osa Peninsulas is very particular. The very narrow continental shelf offshore Cabo Blanco and Punta Chancha (Ch in Fig. 11a), together with Tortuga Island (IT) and the scars from rotational slides offshore Nicoya Gulf and Quepos (Qp), caused four sets of trapped edge waves.
The first set almost covered the Southeast Nicoya Peninsula subsegment (2.2). It extended from Cabo Blanco to Tortuga Island (CB and IT in Fig. 11a) having wave heights over 2 m for up to 10 h. Most tsunamis arrived obliquely to this subsegment ( Fig. 5c and Fig. 6d), which increased the apparent width of the continental shelf decreasing the wave height, despite the low steepness wide continental slope. Most of the beaches in this subsegment (2.2) were assigned to group 3. The only two beaches in group 4 are Tambor and Pochote in Ballena Bay (Ba), in Fig. 10, due to local effects that were not fully modeled here.
The second set of trapped edge waves extended from Tortuga Island to Punta Leona (IT and Le in Fig. 11a) and had wave heights over 2 m for up to 13.5 h. The largest wave heights and over the longest timespan offshore this section of the coastline were located at Nicoya Gulf mouth (between Negritos Islands and Punta Torre, IN and PT, delimited by gray dashed lines at Fig. 11). Further studies are needed to determine if these waves are transmitted inside the Gulf; as mentioned before, Nicoya Gulf is very shallow, which requires higher resolution bathymetry, not available now.
The third set of trapped edge waves extended from Punta Leona (Le) to Quepos (Qp), and the fourth set covered Coronado Bay, from Quepos to Punta Llorona (PL, Fig.  11a). Both had wave heights over 2 m for up to 10 h. The trapped waves at Coronado Bay had the largest periods of up to 1:56 h.
These edge waves caused larger joint maximum tsunami heights at Punta Mala and Marino Ballena National Park (PM and BNP in Fig. 10 and Fig. 11) and caused Hermosa, West, Center and East Esterillos, Palo Seco, Cocal, Guapil (WE, CE, EE, PS, CQ, G in Fig. 10) and beaches within the National Park to be classified in group 4.
The geomorphology of Coronado Bay is similar to that of Tehuantepec shelf, México, where the 2017 tsunami caused trapped waves that lasted for 3 days (Adriano et al. 2018;Chacón-Barrantes 2018;Melgar and Ruiz-Angulo 2018). A comprehensive study of these edge waves in Coronado Bay shelf, like the one performed by Melgar and Ruiz-Angulo (2018), is beyond the scope of this paper, but it should be done to have a more accurate threat estimate. These trapped waves might be dangerous not only for the coastal towns but also for the many fishermen and touristic boats that offer whale watching tours, snorkeling, and visits to Caño Island.

Coco's Ridge-Osa Peninsula segment
The West Osa Peninsula subsegment (3.1) presented the highest tsunami heights of the segment and the highest of the entire coastline together with the Southwest Nicoya Peninsula subsegment (2.1). Subsegment 3.1 had all beaches in groups 5 to 3; beaches in group 5 are Ganado (GA), Ganadito (Ga), Drake, Colorada (Cl), and Llorona (Fig. 10), and the only beach in group 2 is Salsipuedes (SS), at Corcovado National Park. The trapped waves mentioned above at Coronado Bay caused these larger tsunami heights. Also, these waves got refracted around Caño Island (IC in Fig.  6c) interfering positively behind it, causing larger tsunami heights for its eastern side and for mainland beaches eastern from the island (arrow in Fig. 6d). This effect was stronger for tsunamis arriving from the east, like those originated at the Tonga-Kermadec Trench (Fig. 4b) or the C2 scenario (Fig. 7e) from Zamora and Babeyko (2015). However, it was also present for tsunamis arriving from the south-southeast and northnorthwest, like the Colombia-Ecuador scenarios ( Fig. 4a and Fig. 5d) or scenarios close to Nicoya Peninsula (Fig. 7a, b, or d and Fig. 8d), respectively. Subsegments 3.2, 3.3, and 3.4 had beaches mostly in group 2 due to the very steep continental slope that reflected most of the tsunami energy, despite their extremely narrow continental shelf. Directivity of the most threatening distant tsunamis (Colombia and Tonga Kermadec trenches) may be also involved in these low to moderate heights. For both scenarios, tsunamis arrived almost perpendicularly to the coastline of subsegments 3.2 and 3.4 (Colombia scenarios from the east and Tonga scenarios from the west, Fig. 4a and Fig. 5b, respectively), hindering wave shoaling, like happened for subsegment 2.1. The South Osa Peninsula subsegment (3.2) had beaches only in group 2, as it mostly has a straight coastline, but the other two subsegments had a few beaches in group 3.
In Dulce Gulf subsegment (3.3), beaches in group 3 are located at both sides of the gulf's mouth: Carbonera and Sombrero (* in Fig. 10) at the west side and beaches between Cocal Amarillo and Río Claro de Pavones (C-P in Fig. 10) at the east side. The rest of the beaches in the mouth were in group 2, as well as beaches at the end of the gulf, with a region of beaches in group 1 between them. The Gulf normal modes are worth to study with higher resolution bathymetry in future work to discard amplification effects like happened at Palu Bay, Indonesia, in 2018. These normal modes seemed responsible for the well-defined segmentation of beach threat group within the gulf.
In Punta Burica subsegment (3.4), there are only two beaches in group 3: Clarita (Cl in Fig. 10) and Puesto La Playa (PP in Fig. 10). These are explained by refraction on the cape south of the former. This is a very localized effect, offshore these two very long beaches.

Conclusions
This investigation analyzed the hazard from earthquake generated tsunamis for the Pacific coast of Costa Rica, using numerical modeling of the propagation of 57 distant and local tsunami scenarios around the Pacific basin. The maximum tsunami height at 20-m depth was selected among all the scenarios to define the hazard level of each coastal location. Coastal vulnerability was not analyzed since it was beyond the scope of this paper.
Due to directivity, the largest tsunami hazard for the Costa Rican Pacific coast was posed by earthquake sources along the South American and Tonga-Kermadec subduction zones. The former poses a special hazard due to a reduced arrival time of less than 90 min. The results of the modeling performed here of tsunamis coming from Middle America (offshore Mexico), Juan de Fuca (Cascadia), Gulf of Alaska, Aleutian, Kuril, Japan, Nankai, Bonin, Marianas, New Hebrids, Chile, and Peru-Chile trenches showed that Costa Rica does not lie along the direction of maximum energy for those sources. These results agreed with historical records of distant tsunamis in Costa Rica.
Simulation results showed that continental shelf width and continental slope steepness and width exerted a major influence on nearshore tsunami heights. In general, tsunamis traveling over steeper slopes were mostly reflected and transmitted less energy to the shelf. If the shelf is extremely narrow, as happened for the South Osa and Punta Burica subsegments, the continental slope acts almost like a rigid wall that reflects most of the tsunami energy.
In regions where submarine landslides lowered the continental slope steepness and narrowed the shelf, most of the tsunami energy is transmitted to the shelf. In a narrow shelf, the tsunami dissipates less energy before reaching the coast, and then tsunami heights in the nearshore will be large. This happened for the Southwest Nicoya Peninsula subsegment, and at Jaco and Palo Seco, in the Central Pacific subsegment.
However, tsunami interactions with continental slope geomorphology and coastal geomorphology played a major role on tsunami heights as well, like happened for the two subsegments having the largest tsunami heights of the study area. For the Southwest Nicoya Peninsula subsegment, the concave shape of its shelf focused the energy of most of the tsunamis simulated here. For the West Osa Peninsula subsegment, the combination of wave refraction, interference, and trapped edge waves leads to the higher tsunami heights obtained there, particularly east of Caño Island. Trapped edge waves were caused by the narrow shelves at Nicoya and Osa Peninsulas and the geomorphology of the continental slope, platform, and coast between them. Therefore, these effects also influenced tsunami heights at the Central Pacific subsegment.
In general, irregular coastlines led to higher tsunami heights, as happened for the North Guanacaste and North Nicoya Peninsula subsegments, despite coastal effects were not fully considered due to model resolution. Less irregular coastlines led to smaller heights as happened for Western Nicoya Peninsula subsegment and Punta Burica subsegment. Coastal straightness usually reflects active uplifting, which is often related to narrow continental shelf, straight shelf break (in absence of slope landslides), and steep continental slopes, which resulted in lower tsunami heights.
Patchy distribution of maximum tsunami heights along Dulce Gulf subsegment might be related to normal modes. Nicoya Gulf was not considered in this study due to its shallow depths, but tsunami heights obtained at the eastern margin of its mouth might imply moderate or high tsunami heights inside the gulf, together with trapped waves. The authors recommend further studies of both gulfs with high-resolution bathymetry to study their tsunami response, including normal modes.
Off the western margin of Nicoya Gulf, the Southeast Nicoya Peninsula subsegment had moderated tsunami heights due to its main orientation, perpendicular to the trench.
The results obtained here, as well as the beach classification defined, should be considered by stakeholders to increase tsunami preparedness in coastal communities. It is recommended to assess the vulnerability of those locations and to define a tsunami risk index. Furthermore, it is strongly recommended to survey coastal bathymetric data to simulate tsunami inundation of the locations with higher tsunami risk to define the extent of the potential tsunami inundation area. The definition of that area allows to create tsunami evacuation maps, plans, and local protocols.
The bathymetry employed for the modeling was obtained from nautical charts and global bathymetric data (GEBCO-Group 2019). Even though that data allows the characterization of the tsunami hazard nearshore, the results may vary slightly if better and newer bathymetric data is em pl oyed, which is advi sable. Additionally, only seismic sources were considered in this study. Other tsunami sources, such as landslides, might originate larger tsunamis than the simulated here. We recommend the use of non-seismic tsunami sources in future studies of tsunami threat for Costa Rica.
Self-evacuation should be emphasized in the country, together with appropriate definition of differentiated evacuation zones for local and distant tsunamis. Knowing the specific hazard will save lives, as for local tsunamis people will know they have little time to evacuate but the distances are moderate. Also, landslide generated tsunamis should be included in the future for a complete tsunami hazard assessment, as they might represent a larger hazard than the earthquake generated tsunamis.