Application of Scenario Earthquakes for Analysis of Seismically Triggered Landslide Hazard: A Case Study in Costa Rica

Aplicación de escenarios de terremotos para el análisis de la amenaza de deslizamientos disparados por sismos: un estudio de caso en Costa Rica

Dylan M. Seal1*, M. Anna Nowicki Jessee2, Michael W. Hamburger1, Paulo Ruiz3

1Indiana University, Department of Earth & Atmospheric Sciences, Bloomington, USA

2Indiana University-Purdue University Indianapolis, Department of Earth Sciences, Indianapolis, USA

3Universidad de Costa Rica (UCR), Escuela Centroamericana de Geología (ECG) y Red Sismológica Nacional (RSN: UCR-ICE), San Jose, Costa Rica

*Corresponding Author: seald@bc.edu

(Received: 15/03/2022; accepted: 03/05/2022)

Revista Geológica de América Central, 67, 2022, doi: 10.15517/rgac.v67i0.51700

ISSN: 0256-7024

ABSTRACT: In this study, we demonstrate the capabilities of hypothetical scenario earthquakes as a new tool for assessment of hazards associated with earthquake-triggered landslides. Costa Rica offers an ideal environment for demonstrating the utility of scenario earthquakes due to its diverse tectonic environments and associated widespread seismic hazard, rugged topography, and high landslide susceptibility. We investigate the relative influence of landslide proxies such as topographic slope, peak ground velocity (PGV), and compound topographic index (CTI), and earthquake source parameters such as magnitude and depth, on predicted landslide probability and fatality. We examine five distinct tectonic environments, including subduction events beneath the (1) Nicoya and (2) Osa peninsulas respectively, (3) intraplate earthquakes beneath the Central Volcanic Range (CVR) and (4) the Central Costa Rica Deformed Belt (CCRDB), and (5) back-arc thrust events on the eastern Caribbean coast. Our results demonstrate that the slope, PGV, and CTI thresholds necessary to produce landslide probabilities greater than 10% vary by tectonic environment. In all cases, we observe magnitude to be the primary control on the predicted maximum landslide probability and overall areal landslide coverage. We validate model predictions with observed landslide inventories from the 2009 Cinchona and 1991 Limon earthquakes, demonstrating a good fit, where over 70% of landslides occurring in zones of greater than 20% probability. We also use a global model of landslide impact to predict exposure and fatality ranges for each scenario earthquake of this study, revealing that moderate-sized earthquakes in the CCRDB and CVR and large subduction megathrust earthquakes each pose a significant hazard to Costa Rica’s population.

Keywords: earthquakes; coseismic-landslides; scenario; natural hazards; Costa Rica.

RESUMEN: Se muestran las capacidades de la modelación de escenarios de terremotos hipotéticos como una herramienta para la evaluación de la amenaza de deslizamientos cosismicos. El entorno geológico - tectónico de Costa Rica y sus características orográficas resultaron ideales para probar esta técnica. El modelo trabaja con: pendiente del terreno, velocidad máxima del suelo (VMS), índice topográfico compuesto (ITC), magnitud y profundidad del sismo disparador para pronosticar la probabilidad de deslizamientos y sus efectos. Se examinaron 5 ambientes tectónicos, eventos de subducción en las penínsulas de (1) Nicoya y (2) Osa respectivamente, (3) eventos intraplaca bajo la Cordillera Volcánica Central (CVC) y el (4) Cinturón Deformado del Centro de Costa Rica y (5) eventos en el tras-arco en el Caribe. Los resultados muestran variabilidad en los valores de umbral de pendientes VMS e ITC dependiente del ambiente modelado, para producir probabilidades de disparo de deslizamientos superiores al 10 %. La magnitud del sismo es siempre el principal controlador de la probabilidad máxima de ocurrencia de deslizamientos. Se validaron las predicciones del modelo comparando los resultados con inventarios de deslizamientos observados en los terremotos de Cinchona (2009) y Limón (1991), donde más del 70% de los deslizamientos ocurrieron en zonas con una probabilidad > 20%. Finalmente utilizamos un modelo global de impacto de deslizamientos para predecir los rangos de exposición y letalidad para cada escenario, revelando que los terremotos de tamaño moderado en CDCCR y CVC y los grandes terremotos de subducción representan un peligro potencial significativo para la población de Costa Rica.

Palabras clave: terremotos; deslizamientos-cosismicos; escenarios; peligros; Costa Rica.

Introduction

Seismically triggered landslides are a secondary effect of earthquakes that have potential to result in significant fatalities and damage to infrastructure in regions of high relief and seismotectonic susceptibility. When landslides are triggered by an earthquake, they can often overshadow damage caused by ground shaking (Bird and Bommer, 2004). Marano et al. (2009) estimated that 21.5% of all fatal earthquakes result in fatalities from non-shaking secondary causes, such as landslides. Seal et al. (2020) documented 220,284 fatalities directly attributed to EQIL (Earthquake-Induced Landslides) for the time period 1772 to 2020. Similarly, observations demonstrate that the number of earthquake fatalities is positively correlated with the presence of subsequently triggered landslides (Budimir et al., 2014). These observations emphasize the need for robust methods of estimating EQIL distribution and impact modeling for natural hazard mitigation and planning efforts.

There are currently a number of distinct methods for estimating the hazard and spatial likelihood of earthquake-induced landslides on a regional or global scale, which are primarily based on either probabilistic or deterministic seismic hazard assessment techniques. Probabilistic seismic hazard analysis (PSHA) aims to quantify shaking hazard from all possible earthquake events that may affect a region, along with inherent uncertainties related to magnitude, location, and shaking intensity, to compute the probability of exceeding some critical level of ground shaking at a given location over some time interval (Baker, 2008). Alternatively, deterministic seismic hazard analysis (DSHA) focuses on the potential impacts that may result from a single earthquake event. Often this may be a “design earthquake” that is used for evaluating safety of critical structures, such as dams or nuclear power plants (e.g., Rao and Choudhury, 2021).

In this study we concentrate on deterministic scenario earthquakes using the U.S. Geological Survey ShakeMap system (Worden and Wald, 2016), which is used to provide estimates of ground shaking for potential future seismic events at any specified magnitude and hypocenter location. Using these deterministic hypothetical ground shaking estimates, we utilize the global statistical landslide model of Nowicki Jessee et al. (2018), which, in addition to ground motion, incorporates several topographic and geologic variables to predict the landslide probability distribution for any earthquake. Importantly, these scenario earthquakes may be computed and analyzed using the global landslide model anywhere in the world, regardless of historical earthquake and observed landslide inventory availability.

Past studies have used a number of different methods to compute ground shaking and impact estimates for hypothetical scenario earthquakes. Robinson et al. (2018) and Robinson (2020) predicted ground shaking for scenario earthquake ensembles consisting of 90 different seismic events in Nepal and Bhutan, which were then combined with building type and population data to predict potential fatality and impacts to infrastructure. The USGS ShakeMap system has also been used to create scenario earthquakes at multiple locations of interest across the United States, including a M 7.0 earthquake on the Hayward Fault to examine potential impacts from ground shaking in the San Francisco Bay area of California (Detweiler and Wein, 2017). Other studies have applied different deterministic methods to predict ground shaking and related impacts for scenario earthquakes in major European cities, including Istanbul, Lisbon, Thessaloniki, as well as Baku, Azerbaijan (Spence et al., 2007; Babayev et al., 2010).

In addition to their value in analyzing shaking-related impacts, several studies have used scenario ground motion estimates with various statistical and physical models to predict earthquake-induced landslide occurrence and associated impacts. Robinson et al. (2016) used the fuzzy logic-based statistical coseismic landslide model of Kritikos et al. (2015) in conjunction with a deterministic scenario M 8.0 earthquake on the Alpine Fault of New Zealand to create a landslide hazard map for the region. Probabilistic ground motion estimates for earthquakes of variable recurrence intervals have also been used in conjunction with Newmark’s (1965) sliding block model to identify areas of instability and predict landslide-related impacts in multiple regions of Italy, Russia, Mexico, and Spain (Bozzano et al., 2013; Konovalov et al., 2019; Martino et al., 2018; Martino et al., 2020; Rodriguez-Peces et al., 2011; Salinas-Jasso et al., 2017). Likewise, deterministic methods for estimating ground shaking for individual scenario earthquakes have been alongside Newmark’s model to spatially predict landslide occurrence and impacts in China, India, and Spain (Kumar et al., 2021; Liu et al., 2018; Nayek and Gade, 2021; Rodriguez-Peces et al., 2011). Using our approach, the USGS ShakeMap system was used in conjunction with the global landslide model of Nowicki Jessee et al. (2018) to investigate the impact that both predicted ground shaking and landslide distribution may have on population and key infrastructure using three scenario earthquakes in Nepal (GHI, 2018).

In this study, we create suites of scenario earthquakes to investigate the potential for EQIL in Costa Rica. This Central American country is ideally suited to this study based on its proximity to a major subduction plate boundary, its historical record of frequent and severe earthquake and secondary impacts, widespread landslide susceptibility due to rugged topography and highly variable climate, and local observational data that may independently validate our landslide model predictions (Fig. 1). We systematically explore a range of plausible earthquake source parameters (magnitude and depth) for each scenario location, based on the historical coseismic landslide data of Ruiz et al. (2019), which enables us to examine the range of probable impacts that each region may experience after inevitable future earthquakes occur. Using this approach, areas of high landslide susceptibility can be readily identified for any region of seismic potential. We argue that this methodology can serve as a powerful tool for aiding in the identification of susceptible areas requiring additional hazard and mitigation planning efforts, and for comparison of landslide likelihood across a variety of tectonic, geographic, and geologic environments, critically in regions lacking high-quality data.

Seismotectonic and geographic setting

With its diverse tectonic environments and associated widespread seismic hazard, rugged topography, and extensive areas of high landslide susceptibility, Costa Rica offers an ideal environment for demonstrating the capability of scenario earthquakes to coseismic landslide risk assessment. We examine five environments to represent the range of hazard presented by the tectonic, climatic, and lithologic regions present in Costa Rica (Fig. 1). These include (1) subduction megathrust events along the Cocos-Caribbean plate boundary beneath the northern Nicoya peninsula and (2) the southern Osa peninsula, (3) intraplate events beneath the Central Volcanic Range (CVR) and (4) the Central Costa Rica Deformed Belt (CCRDB), and (5) back-arc thrust events on the eastern Caribbean coast of Costa Rica.

Cocos-Caribbean Subduction Zone

Along the western coast of Costa Rica lies the Cocos-Caribbean convergent plate boundary where the Cocos plate subducts beneath the Caribbean plate at an average velocity of 83 mm/yr (DeMets et al., 2010). The morphology of the subducting Cocos slab varies greatly along strike, from areas of rugged seamount topography to others with relatively smooth seafloor (von Huene et al., 2000). Each region of high seafloor relief coincides with clusters of observed seismicity, one on the southeast coast of the northern Nicoya peninsula and one beneath the southern Osa peninsula (von Huene et al., 2000). We analyze subduction-related seismicity and associated landslide hazard for the northern (Nicoya) and southern (Osa) flanks of Costa Rica’s subduction system, as each have historically demonstrated capability of producing large magnitude 7.0+ subduction megathrust earthquakes that trigger landslides (Ruiz et al., 2019; Tajima and Kikuchi, 1995).

The Wadati-Benioff zone is observed to be highly variable along strike, extending to a maximum depth of 180 km in the north, and to no greater than 50 km in the south, hypothesized to be indicative of along-strike variation in the age of the subducting Cocos plate (Protti et al., 1994; Suarez et al., 1995). The northern segment of the subduction zone last ruptured in 2012 during the Ww 7.6 Nicoya earthquake at a depth of 15.4 km, approximately 10 km offshore of the peninsula (Ruiz et al., 2019; Yue et al., 2013). Seismological data from this event reveal that many areas adjacent to the megathrust rupture, including an offshore locked patch of the plate interface, experienced no significant coseismic slip or aftershocks and thus remained locked. This suggests the potential for similar-sized earthquakes on these fault surfaces in the future (Yue et al., 2013). The southern segment of the subduction zone ruptured directly beneath the western coast of the Osa peninsula and triggered landslides during the 1983 Golfito Ww 7.1 earthquake (Ruiz et al., 2019). This occurred in close proximity to the epicenters of the 1904 Ms 7.6 and 1941 MS 7.5 subduction earthquakes (Tajima and Kikuchi, 1995), where the aseismic Cocos Ridge is thought to subduct beneath the Caribbean plate (Adamek et al., 1987).

Central Volcanic Range

The Central Volcanic Range (CVR) contains a series of active Quaternary volcanoes within the central region of Costa Rica and is associated with a zone of crustal seismicity spanning the length of the country. The area is a primary focus of seismic hazard research due to close proximity with the capital city San Jose and related populated centers such as Cartago. The region is very susceptible to slope failure following seismic events due to its high relief, high mean annual precipitation (3 - 6 m/yr), extensive hydrothermal alteration and weathering, and high rates of erosion (Ruiz et al., 2019). Among the twenty earthquake-induced landslide events presented in the Ruiz et al. (2019) catalog, ten occurred near volcanic centers, the deadliest of which was the 2009 Cinchona Ww 6.1 earthquake near Poás Volcano, approximately 50 km north of San Jose. That event resulted in 30 fatalities directly attributed to landsliding. Importantly, the CVR is characterized by multiple networks of active shallow crustal faults, concentrating historical seismicity near Poás Volcano and the Central Valley near Cartago (Fernández and Rojas, 2000).

Central Costa Rica Deformed Belt

The Central Costa Rica Deformed Belt (CCRDB), located along the western margin of the Panama microplate, is a zone of active intraplate deformation that spans the width of central Costa Rica and coincides spatially with subduction of the Cocos Ridge (Fisher et al., 1994; Marshall et al., 2000). Marshall et al. (2000) suggested that the crustal shortening resulting from convergence between the Caribbean plate, and rough and thick oceanic crust of the Cocos Ridge, is primarily accommodated through motion along the many strike-slip faults of the CCRDB. Carvajal-Soto et al. (2020) proposed that faults in the CCRDB have the potential to rupture as Ww 7.0+ earthquakes. The Agua Caliente Fault is a left-lateral strike-slip fault within the CCRDB that is especially problematic to the population of central Costa Rica given its proximity with major population centers such as Cartago and San Jose. The 1910 Cartago MS 6.4 earthquake, which resulted in 600 fatalities in Cartago (population of 12,000 in 1910) and destroyed many buildings, took place along the Agua Caliente Fault (Alonso-Henar et al., 2013).

Arc thrust events on the eastern Caribbean coast of Costa Rica.

Back-Arc Convergence

The convergent tectonics of the Cocos-Caribbean boundary are linked with a back-arc thrust regime along the eastern Caribbean coast of Costa Rica and Panama, hypothesized to contain a dense network of imbricate thrust and strike-slip faults (Suarez et al., 1995). In addition to active seismicity, the region hosts the Talamanca Cordillera, the highest mountain range of Central America, and extensive unconsolidated sediment along the coast. This combination makes the region susceptible to widespread earthquake-triggered landslides and liquefaction (Suarez et al., 1995).

Historically, the region has experienced large earthquakes with devastating secondary impacts, notably the 1991 Limon Ww 7.6 earthquake, which triggered landslides over an area of approximately 2000 km2 in close proximity with Puerto Limon and other population centers (Mora and Mora, 1994). Geophysical observations suggest that subduction of the Cocos Ridge increases plate coupling, resulting in part of the Cocos plates relative motion being transferred and accommodated by back-arc crustal shortening (Suarez et al., 1995). The North Panama Deformed Belt (NPDB; Fig. 1) also inhabits the region directly offshore of southeastern Costa Rica and may be attributed to historical ground shaking in the back-arc. The largest earthquake recorded within the NPDB occurred in 1882 and is interpreted to have a magnitude of 8.3 (Mendoza and Nishenko, 1989). While the exact nature of the NPDB convergence is unresolved, evidence suggests that active subduction of the Caribbean plate beneath the Panama microplate may be taking place (Camacho et al., 2010).

Methods

ShakeMap Scenario Earthquakes

Reliable estimates of strong ground motion are essential for estimating the likely spatial distribution of earthquake-triggered landslides. We use the ShakeMap software (Worden and Wald, 2016) to compute spatial estimates of ground shaking for both observed and anticipated future scenario earthquakes. Scenario earthquakes are a powerful approach to deterministic seismic modeling as they allow for in-depth hazard assessment in areas with limited records of historical seismicity or lacking in high-quality seismic data. ShakeMap allows for the creation of scenario earthquakes at any arbitrary location, depth, magnitude, and fault geometry. In this study, we create scenarios across the range of tectonic environments present in Costa Rica to assess the predicted distribution and impact of both ground shaking and earthquake-induced landslides.

In order to ground our predictions in observational data, we use the database presented by Ruiz et al. (2019), which provides a detailed record of all landslide-inducing earthquakes in Costa Rica dating back to 1772. Specifically, we focus on two events with relatively well constrained inventories of earthquake-triggered landslides: the 1991 Limon earthquake (Hernandez et al., 1992; Mora and Mora, 1994) and the 2009 Cinchona event (Ruiz et al., 2019). We statistically compare our model predictions with observational data from each inventory, which is presented in the discussion. Lastly, we use the historical seismicity and fault systems identified in Ruiz et al. (2019) to place scenario earthquake events at geophysically plausible locations, magnitudes, and depths to represent the hazards associated with each major tectonic region of Costa Rica.

Earthquake-Induced Landslide Modeling

In addition to assessing the distribution and impact of ground shaking, ShakeMap scenario earthquakes can be used in conjunction with existing earthquake-induced landslide models to investigate controls, distribution, and impacts of ground failure from future earthquakes. We use our ShakeMap scenario earthquake ground motion estimates, in conjunction with the global empirical landslide model presented in Nowicki Jessee et al. (2018), to estimate the spatial probability distribution of landslides after an earthquake. The correlation and spatial distribution of landslide probability with each ground failure proxy can be examined to determine the relative influence of each factor on earthquake-induced landsliding. We also compute landslide areal percentages for each grid cell using the relationship outlined in Nowicki Jessee et al. (2018), interpreted to represent the portion of each cell predicted to contain landslides. We sum all areal percentages for each scenario to obtain the landslide hazard index (LHI; Eq. (1)), where LA is the predicted percentage of areal coverage for each grid cell and n is the number of pixels with predicted ground shaking (from each ShakeMap). LHI is a metric that we use to quantify and compare the overall landslide area triggered by each respective earthquake.

(1)

Selection of Scenario Earthquakes

The selection of a plausible earthquake hypocenter and magnitude is essential to realistically assessing the seismic hazard associated with a given fault. Therefore, we primarily utilize historical earthquake-induced landslide data from the database presented in Ruiz et al. (2019) as a basis for the location and depth of scenario earthquakes created in each of the five tectonic environments analyzed in this study. Historical earthquake moment magnitudes (Mw), and hypocenters not included in Ruiz et al. (2019) are derived from the USGS PAGER-CAT (Allen et al., 2009) and ANSS ComCat (USGS, 2017) systems. In order to investigate the range of hazard possible from a large interplate subduction megathrust earthquake, we examine events located along the Cocos-Caribbean plate interface. The location of the plate interface is obtained from the

Slab 2.0 model (Hayes et al., 2018), to a maximum depth of 50 km, which we interpret to be the maximum depth for events on the subduction plate interface. This geometry is used as the basis for identifying earthquake hypocenters along two transects of the plate boundary, beneath both the northern Nicoya and southern Osa peninsulas. So as to encompass the widespread landslide hazard near population centers of central Costa Rica, we create scenario earthquakes within the CVR, on the Angel Fault of Poás Volcano, northwest of San Jose at the epicenter of the 2009 Cinchona Ww 6.1 (Ruiz et al., 2019), and within the CCRDB, on the Agua Caliente Fault, directly south of Cartago at the epicenter of the 1910 Cartago MS 6.4 earthquake (Alonso-Henar et al., 2013). Back-arc scenario earthquakes are placed at the epicenter of the 1991 Limón Mw 7.6 earthquake, on the eastern Caribbean coast of Costa Rica directly south of Puerto Limón.

Fatality Modeling

Because of the relatively small area and concentration of population in the central valley of Costa Rica (Alonso-Henar et al., 2013), where high seismic hazard and numerous landslide-susceptible landscapes are present, earthquake-induced landslides in Costa Rica often occur in close proximity with population centers. Ruiz et al. (2019) recorded 98 fatalities that can be directly attributed to seismically triggered landslides, which underscores the need for improved methods for assessing the spatial distribution of landslide exposure. We use Landscan population data (Bright et al., 2014) at 30 arc-second (~1 km) resolution to compute population exposure for each scenario earthquake. We calculate the probabilistic landslide exposure value (PLEV) detailed in Nowicki Jessee et al. (2020) using these population density estimates, alongside scenario landslide probability, using Eq. (2):

(2)

where Po is the population count per grid cell, PL is the landslide probability per grid cell, and n is the number of pixels with predicted ground shaking for a given scenario earthquake. Nowicki Jessee et al. (2020) contains two regressions that relate PLEV to predicted fatality, one that is trained using observational data from 91 historical earthquakes, and another that incorporates the United Nations Human Development Index (UNHDI) as a proxy for vulnerability. We also compare PLEV and landslide fatality observations for these historical earthquakes (Nowicki Jessee et al., 2020) with our scenario earthquake predictions to determine the context of potential earthquake-induced landslide hazard in Costa Rica relative to global impact observations.

Results

Controls on Earthquake-Induced Landslides

Subduction Megathrust Earthquakes

Large subduction megathrust earthquakes along the Cocos-Caribbean interface have potential to produce devastating impacts and widespread landsliding throughout Central America, prominently Costa Rica’s Pacific coast. Historical subduction seismicity has been concentrated beneath the northern (Nicoya) and southern (Osa) peninsulas (Ruiz et al., 2019; Tajima and Kikuchi, 1995; von Huene et al., 2000; Yue et al., 2013). We therefore choose to analyze scenario earthquakes along the Cocos-Caribbean plate interface of each region to elucidate the relative influence of earthquake source and geomorphologic factors on ground failure in Costa Rica. Both regions are especially vulnerable, due to their proximity to the plate boundary and rugged relief. In order to constrain physical controls on landslides triggered by subduction-related seismicity, we compare the relative influence of ground failure proxies, including topographic slope, peak ground velocity (PGV), and compound topographic index (CTI), on the predicted landslide probability distribution for five suites of scenario earthquakes (e.g., Fig. 2). The level and distribution of ground shaking caused by rupture on the subduction interface is highly dependent upon earthquake location, magnitude, and depth. Earthquake depth and location are closely linked by the geometry of the subduction interface and dictate whether the subduction interface ruptures offshore (beneath the Pacific Ocean) or beneath mainland Costa Rica. This in turn controls shaking intensity and severity of impacts to population. Example landslide probability maps for M 7.5 subduction earthquake scenarios at 40 km depth beneath the Nicoya and Osa peninsulas are shown in Figure 3.

We first compare scenario earthquakes beneath the Nicoya and Osa peninsulas at M 7.5 and 40 km depth. Analysis of the effect of these proxies on predicted landslide probabilities reveals numerous similarities between scenario megathrust events in the north and south of Costa Rica (Fig. 4). Foremost, we observe that areas of high landslide probability are strongly associated with both high slope and PGV, and likewise predicted landslide probabilities are low when either one of PGV and slope do not exceed some threshold value. Landslide probabilities above 10% typically require a threshold value of ~5° slope and ~6 cm/s PGV in Nicoya, and ~5° slope and ~4 cm/s PGV in Osa (Fig. 4). Interestingly, the Osa peninsula suite of scenarios produce higher predicted maximum PGV and landslide probability values than the corresponding Nicoya scenarios, despite a similar range of slope values. The Cocos-Caribbean interface beneath the Osa peninsula is situated further inland at each depth than under Nicoya, which may account for the higher observed ground motions and landslide probabilities of the Osa scenario suite. CTI (relative wetness) plays a more complicated role in controlling predicted landslide probabilities, with higher probabilities generally occurring at low CTI values, below ~5, when the PGV the slope thresholds are met (Fig. S3 and S4 in the Supplementary Materials). When CTI is compared to slope, CTI appears to exhibit only a secondary control on predicted landslide probability. These results suggest that in each megathrust setting of Costa Rica, slope and PGV are the dominant factors that control predicted landslide probability, with CTI playing a secondary yet significant role.

The magnitude and depth at which a seismic event occurs on the Cocos-Caribbean interface have a direct control on the location of the earthquake epicenter, and thus predicted coseismic landslide probability. To evaluate this relationship, we create two suites of scenario earthquakes along transects of the northern Nicoya and southern Osa portions of Costa Rica’s subduction system, at plausible depths and magnitudes determined by historical seismicity (see methods) to encapsulate the respective landslide hazard of each region. In order to better understand the control that magnitude has on predicted landslide probability in the Nicoya and Osa regions, we also compare scenario earthquakes at 20 km depth along the Cocos-Caribbean interface, with M 7.0, M 7.5, and M 8.0. In both subduction settings, when hypocenters are held constant, the areas of high predicted landslide probability increase with magnitude (Fig. 5). Similarly, we observe the LHI and maximum landslide probability to increase with increasing magnitude, each being largest for the M 8.0 subduction scenarios (Supplementary Materials). The LHI of each event increases nearly linearly with magnitude, whereas the maximum landslide probabilities increase logarithmically. We note that LHI and maximum landslide probability are not identical between the Nicoya and Osa subduction segments for scenario events of the same magnitude and depth, likely representing the differences in subduction plate interface location and topography associated with each region.

For all subduction scenarios, we find that earthquakes beneath the Osa peninsula produce larger LHI values, and thus are predicted to trigger more landslides, than equivalent earthquakes beneath the Nicoya peninsula. In all likelihood, the primary driver of these differences is the location of the Cocos-Caribbean plate interface. The subduction interface is located further inland at each respective depth beneath the Osa peninsula relative to Nicoya, resulting in a larger land area proximal to the earthquake epicenter where the highest intensity of ground shaking is located. As each peninsula is characterized by slope and hydrologic conditions necessary for landslide occurrence, the larger relative area of land subjected to high levels of ground shaking would explain the larger LHI values observed for the Osa scenario suite.

In addition to magnitude, depth also plays a significant, albeit more complex role, in controlling the predicted landslide probability distribution. Because we compute subduction interplate scenario earthquakes along the Cocos-Caribbean plate interface, changes in depth correlate with changes in epicenter location, moving progressively towards mainland Costa Rica as depth increases. To test the role of earthquake depth on predicted landslide probabilities, we hold magnitude constant at M 7.5 and vary depth (and thus position along the plate interface at 20, 40, and 50 km.

We observe that scenarios at 20 km depth produce a significantly larger cumulative area of high landslide probabilities, above 20%, than the 40 and 50 km scenarios, which produce marginally larger areas with landslide probability less than 20% (Fig. 5). Likewise, the 20 km events produce the largest maximum landslide probability in both the Nicoya and Osa regions, though the second largest maximum landslide probability results from the 40 km scenario in Osa, and the 50 km event in Nicoya. We conclude that increasing earthquake depth on the subduction interface results in a larger cumulative area of low landslide probability below 20%, and a smaller area affected by high landslide probabilities.

For all scenarios, we observe LHI to increase with depth, and note that LHI values are nearly identical between regions for 40 and 50 km events, only differing significantly for the 20 km set of scenarios, where LHI is higher in Osa (Fig. 6). This likely reflects the position relative to the coast. These results demonstrate that while an increase in depth decreases the maximum landslide probability and areas of highest landslide probability, it conversely increases LHI, and thus the overall spatial likelihood of landslide occurrence is positively correlated with earthquake depth on each subduction interface.

Shallow Crustal Earthquakes in the Central Volcanic Range and Central Costa Rica Deformed Belt

Moderate to large earthquakes in the CVR and CCRDB have produced devastating shaking and secondary impacts throughout recorded history, dating back to 1772 (Ruiz et al., 2019). We analyze two suites of scenario earthquakes located on the Angel (CVR) and Agua Caliente (CCRDB) faults, near the epicenters of the 2009 Cinchona Ww 6.1 and 1910 Cartago MS 6.4 earthquakes, respectively (see methods). At both locations, we compare the influence of each ground failure proxy on landslide probability for M 6.5 earthquakes at 5 km depth (Fig. 4; Fig. S1 and S2 in the Supplementary Materials). We then investigate the relative influence of earthquake source parameters, including earthquake magnitude and depth, on maximum landslide probability and overall areal landslide likelihood (Fig. 5 and 6).

To the first order, landslide occurrence is governed by the interaction of ground shaking with local topographic and climatic variables. Therefore, we analyze the relative influence that ground motion, slope, and relative wetness has on landslide probability for events located along the Angel and Agua Caliente faults of the CVR and CCRDB, respectively (Fig. 4). We find that both high PGV and slope are required for high predicted landslide probability, with the largest landslide probabilities being produced exclusively when both PGV and slope are highest. Secondarily, CTI correlates negatively with landslide probability, albeit with much more variability than the correlation of PGV or slope with probability (Fig. S1 and S2 in the Supplementary Materials). We observe high landslide probabilities to be generally concentrated in areas of low CTI, with heavy dependence on the PGV and slope values present. Very high landslide probability in the CVR (P > 50%) are associated with areas of relatively high (> ~10°) slope and PGV (> ~25 cm/s). Moderate to high landslide probabilities (P > 10%) require slopes greater than ~3° and PGV above ~8 cm/s. Moderate to high landslide probabilities, above 10% probability, require slopes greater than 3° and PGV above 8 cm/s. For both the Cinchona and Cartago locations, very high landslide probabilities are most concentrated at CTI values below 6; however, CTI appears to exhibit less of a distinct control on probability than either slope or PGV.

Here we investigate the influence of earthquake source parameters on ground failure, varying earthquake magnitude and depth for suites of scenario earthquakes at two locations in the CVR and CCRDB. Each of these two shallow intraplate scenario earthquake suites includes M 5.0, M 5.5, M 6.0, and M 6.5 at 5 km depth, as well as M 6.5 events at 10 km and 15 km depth. For both suites of earthquakes, we find, not surprisingly, that magnitude is the primary control on the distribution of predicted landslide probability. Increasing earthquake magnitude always results in larger areas affected by each respective landslide probability value, as well as higher maximum landslide probability and LHI values (Fig. 5). Of these two intraplate suites, the Cinchona events always produce a significantly higher maximum landslide probability than the corresponding Agua Caliente scenario, while LHI is nearly identical and only marginally higher for Cinchona (Fig. 6).

Despite the relatively small (~40 km) distance between the two faults, they exist in significantly different climatic regions. The topography of the CVR divides Costa Rica into two distinct orographic regions; the Pacific slope to the southeast of the CVR, which receives 3 – 4 m/yr precipitation with a dry season from December to April, and the Caribbean slope to the northwest of the CVR, which receives 4 – 6 m/yr precipitation, with no dry season (Ruiz et al., 2019). Previous studies have demonstrated that monthly precipitation correlates positively with landslide occurrence in each region (Quesada-Román et al., 2021). Observations from the 2009 Cinchona earthquake demonstrated that 95% of the mapped landslides were located on the Caribbean slope (Ruiz et al., 2019). As the Angel Fault is located on the Caribbean slope and the Agua Caliente Fault is located on the Pacific slope, distinct climatic and hydrologic differences may exist between these two scenario locations, which may be partially accounted for in the Nowicki Jessee et al. (2018) landslide model by CTI. Predicted differences in maximum landslide probability and LHI likely reflect the differential distribution of precipitation, land-cover types, steep slopes, and lithologies between the two locations.

Depth also plays an important secondary role in determining the intensity and extent of ground shaking, and thus controls landslide occurrence. We observe that increasing depth causes the LHI to increase, where the deepest set of scenarios, at 15 km depth, are predicted to trigger the most extensive areas of landsliding (Fig. 6). On the other hand, maximum landslide probability decreases as the scenario earthquake hypocenter deepens for both locations, as does the cumulative area of landslide probabilities above 20% (Fig. 5; Supplementary Materials). The relationship of LHI and depth identified in the CCRDB and CVR match those of each subduction scenario suite.

Validation of Results with Observed Landslide Inventories

As the Nowicki Jessee et al. (2018) landslide model is based on global data, it is crucial that the model predictions are grounded in local observational data from Costa Rica. In this study, we use mapped landslide inventories from the 2009 Cinchona Ww 6.1 and 1991 Limon Ww 7.6 earthquake to validate model predictions in the regional geophysical context of the CVR and back-arc thrust systems of Costa Rica. We compare the observed mapped landslide polygons of each inventory with predicted landslide probability distributions (Fig 7). For both inventories, we find that landslide observations fit the model predictions well, with landslides generally occurring in areas of high landslide probability, above ~30%. Detailed statistical discussion of model validation for these two events can be found in the discussion.

Cinchona Inventory

The 2009 Cinchona Ww 6.1 earthquake resulted from a 72 km² rupture of the strike-slip Ángel fault on the eastern flank of Poás Volcano, an active volcano located within the CVR, approximately 40 km northwest of San José (Fig. 1; Ruiz et al., 2019). Three months after the earthquake, 4,846 coseismic landslides were mapped over an area of 519 km2 (Fig. 7) using a combination of airborne light detection and ranging (LiDAR), orthophoto data, and field observations (Ruiz et al., 2019). In addition to Poás Volcano, landslides were also observed on two dormant volcanoes, Platanar and Barva, located directly west and east of Poás, respectively (Ruiz et al. 2019). We compare predicted landslide probabilities with these mapped landslide polygons (Fig. 7) to evaluate the performance of the global landslide model. As predicted by the landslide model, nearly all landslides are located in areas with very high peak ground velocity (PGV) and steep slope (Fig. S5 in the Supplementary Materials). We observe the predicted landslide probability distribution to match observed landslide occurrence well in the areas adjacent to the epicenter. Nearly all mapped landslides are found in areas of landslide probability greater than 20%, excluding a small area directly west of the epicenter where disrupted surface landcover inhibited aerial landslide observations. We note a few small areas of predicted high probability (e.g., red patches west of -84.7° W) where landslides were not observed. However, few other areas of high landslide probability were landslide-free, and the majority of observed landslides were in relatively high-probability zones.

Limon Inventory

Following the 1991 Limon Mw 7.6 earthquake, landslides were mapped (Fig. 7) using a combination of oblique aerial photography and observation (Hernandez et al., 1992; Mora and Mora, 1994), which were then georeferenced and digitized as part of this study. This digitized landslide map contains 1994 landslides over a total area of 190.31 km², with an estimated volume of 15.93 km³ of material being displaced. While no fatalities attributed to landsliding were recorded, there was significant damage to infrastructure during the earthquake, and a maximum uplift of 1.5 m was observed near Puerto Limon (Suarez et al., 1995). The earthquake epicenter was approximately 35 km south of Puerto Limon; however, the majority of landslides were observed in the rugged topography to the west and southwest of the city. We compare these mapped landslides to our predicted model probabilities (Fig. 7), which reveals that most landslides occur at landslide probabilities above 50%, and no landslides are found in areas of less than 20% probability. These results indicate that model predictions match landslide observations reasonably well for this earthquake, albeit with a significant degree of overprediction that will be discussed in the discussion. There are large areas of moderate to high probability where landslides were not observed (e.g., areas west and south of the epicenter), suggesting that the model may be over-predicting landslide occurrence in some areas.

Fatality Modelling

The results of our EQIL scenarios described above can also contribute to assessment of likely impacts. PLEV and fatality ranges for each scenario earthquake in Costa Rica are included in the Supplementary Materials (see methods). We examine these scenario impact estimates to investigate the role that magnitude and depth play in the potential societal impacts of an earthquake in each tectonic environment (Fig. 8; Fig. S6 in the Supplementary Materials). It is critical to note that the fatality regression used to obtain these estimates is based on global observations and contain high levels of uncertainty. Earthquake-induced landslide fatalities vary widely from event to event; therefore, our fatality calculations should be considered only as order-of-magnitude estimates.

Regardless of tectonic environment, increasing magnitude corresponds to a logarithmic increase in PLEV, and thus predicted fatality (Fig. 8). Interestingly, the impact of depth is not constant across or within tectonic environments, even for suites of scenarios with constant epicentral locations (Fig. S6 in the Supplementary Materials). For both subduction settings, as earthquakes deepen, predicted exposure and fatality grow, which is to be expected as increasing depth on the Cocos-Caribbean interface both increases the LHI and moves the epicenter closer to the more heavily populated central region of Costa Rica. By contrast, we observe that increasing depth in the CVR and CCRDB increases PLEV for events at the Cinchona earthquake epicenter yet decreases predicted exposure for scenarios on the Agua Caliente Fault. Multiple key differences between the Angel and Agua Caliente faults may contribute to this anomalous depth relationship, including: (1) larger areas of variable relief adjacent to the Angel Fault and thus differential landslide triggering between the two sites, (2) increased distance of the Angel Fault to major population centers such as San Jose, and (3) lithologic variations.

Critically, all scenario earthquakes investigated in this study demonstrate potential for significant landslide exposure to the population of Costa Rica. We observe that even for the smallest LHI scenarios considered in this study, M 5.0 and M 5.5 earthquakes in the CVR, PLEV is above 105, equating to a likelihood of at least several predicted fatalities. While the set of M 5.0 events in the CVR produce the smallest PLEV values, each of the M 5.5 events result in higher PLEV than M 7.0 subduction scenarios beneath the Nicoya and Osa peninsulas, and the M 5.5 event on the Agua Caliente fault produces only marginally lower PLEV than an M 7.5 subduction earthquake. This result is critical to hazard planning and disaster mitigation efforts as these moderate-sized crustal earthquakes in central Costa Rica more frequently trigger landslides than subduction-related earthquakes (Ruiz et al., 2019). Likewise, M 6.0 earthquakes in the CVR and CCRDB at 5, 10, and 15 km produce relatively small areas of landslide hazard (~103 LHI), yet result marginally higher PLEV than M 7.5 subduction earthquakes, which produce significantly larger landslide areas (~103.6 LHI). The same is true for M 6.5 earthquakes in the CVR, which result in marginally higher PLEV than M 8.0 subduction megathrust earthquakes, despite significantly less areal landslide coverage. The Limon M 7.6 Depth 10 km scenario results in nearly the same LHI value as M 7.0 subduction scenarios, but produces PLEV similar to M 6.0 CVR or M 7.5 subduction scenarios. We interpret these results to reflect the higher proportion of the Costa Rican population that resides in close proximity to the CCRDB and CVR, and the relative low population density around the Nicoya and Osa peninsulas, demonstrating that both moderate-sized shallow crustal earthquakes and large subduction megathrust earthquakes in Costa Rica have potential to result in devastating impacts.

Discussion

Earthquake-induced landslides have a long-recorded history of producing devastating impacts in Costa Rica (Ruiz et al., 2019) due to the combination of widespread seismic hazard, and close proximity of population centers with active faults and rugged topography. The five scenario earthquake suites presented here encompass the range of tectonic, geologic, and climatic conditions of Costa Rica. Our scenario landslide predictions confirm that multiple factors govern landslide likelihood, foremost peak ground velocity, topographic slope, earthquake magnitude, and secondarily earthquake depth and compound topographic index. Lithology and land cover are incorporated into the landslide model and are likely important constraints on ground failure; however the global databases used for these two proxies may not differentiate landslide susceptibility in sufficient resolution to distinguish areas of lesser and greater landslide probability. Importantly, all scenario earthquakes demonstrate the potential to generate significant levels of exposure to landslide hazard, and consequently the potential for landslide-related fatalities.

Comparison with Global Models

Our primary result, that landslide likelihood is highest with a combination of strong ground shaking, steep topographic slope, and low compound topographic index (CTI), is confirmed by multiple global and regional statistical analyses (Keefer, 2002; Harp et al., 1981; Martino et al., 2018; Bozzano et al., 2013; Salinas-Jasso et al., 2017; Liu et al., 2018). Liu et al. (2018) used scenario earthquake ground shaking estimates with Newmark’s model to demonstrate that zones of high landslide risk were concentrated in loess or mudstone lithologies in north-central China. Salinas-Jasso et al. (2017) used scenario earthquakes in northeastern Mexico, showing that only slopes steeper than 40° are likely to fail in dry conditions, though this slope threshold angle would lower as ground saturation increased.

Harp et al. (1981) analyzed the landslides triggered by the 1976 Guatemala M 7.5 earthquake to propose that landslide occurrence was governed by seismic intensity, lithology, slope, topographic amplification of ground motion, fracture systems, and the presence of pre-earthquake landslides. Costa Rica is tectonically and climatically similar to Guatemala; therefore, we expect that these factors are likely to exert a control on earthquake-induced landslide occurrence there as well. Following the 2009 Cinchona earthquake, Ruiz et al. (2019) observed that 95% of the coseismic landslides were in the Caribbean orographic region of Costa Rica, despite the epicenter being located ~4 km north of the Caribbean-Pacific orographic slope border, suggesting that increased ground saturation associated with higher rainfall may have influenced landslide distribution and occurrence. These global and regional observations are consistent with the results of this study, while highlighting the need to incorporate of detailed climatic and lithologic variables to encompass other important controls on coseismic landsliding that are likely present in Costa Rica.

We demonstrate that overall landslide extent increases alongside earthquake magnitude and depth for each of the five tectonic environments of Costa Rica analyzed (Fig. 7). Past scenario EQIL studies in India, Mexico, and China corroborate our findings related to magnitude, where the highest magnitude scenario earthquake of each study results in the highest number and area of predicted landslides (Kumar et al., 2021; Salinas-Jasso et al., 2017; Liu et al., 2018). These findings are also corroborated by global landslide observations presented in Keefer (2002), which demonstrate that both the area affected by landslides and maximum distance of landslides from the epicenter increase with earthquake magnitude. In addition, we identify the complex relationship between earthquake depth and landslide distribution, where depth can serve either to increase or decrease landslide likelihood in different tectonic environments, chiefly in subduction settings, where depth controls earthquake location.

To our knowledge, no studies have previously used scenario earthquakes to investigate the complex relations between earthquake depth and coseismic landslide occurrence. While a number of studies have used scenario earthquakes to examine human and infrastructure impacts of major earthquakes (Wyss, 2006; Robinson et al., 2018; 2020) this is the first study that has applied scenario earthquakes specifically to quantify predicted EQIL exposure and impacts.

Global Context of Costa Rican Impact Estimates

Nowicki Jessee et al. (2020) presents both a global average fatality model, and a model that incorporates the United Nations Human Development Index (UNHDI) as a proxy for vulnerability, which we use to predict exposure and fatality ranges for each scenario earthquake. As summarized in Figure 8 and Supplementary Materials Tables S2 and S3, our median predicted fatality values are relatively low, ranging from 3 to 15 (UNHDI model) or 7 to 65 (global average model) fatalities. It is important to note there are very large uncertainties in these models. For example, the upper bound (95% confidence interval) estimates range from 2,320 to 15,600 fatalities (UNHDI) or 4,565 to 49,546 (global) fatalities, demonstrating the high levels of uncertainty inherent in our fatality estimation method – and the potential for devastating impacts if a major landslide were to strike a populated area.

In order to place Costa Rican fatality estimates into a global context, we compare our scenario probabilistic landslide exposure value (PLEV) and predicted fatality results with PLEV and observed landslide fatalities for the 91 historical earthquakes included in the Nowicki Jessee et al. (2020) catalog. The catalog includes events dating back to 1923 for which fatalities directly attributed to EQIL were recorded (Fig. 8). We observe that all of the PLEV values for the 23 scenario earthquakes presented here are above 105, placing our predicted exposure and fatality estimates in the zone of “significant impact” identified by Nowicki Jessee et al. (2020). On the other hand, our predicted fatality values for Costa Rica plot within the lower half of globally observed landslide fatalities. The largest Costa Rican PLEV values calculated in this study are most similar in their estimated exposure to global events including the 2010 Chile M 8.8 and 1923 Japan M 8.1 earthquakes, which resulted in 4 and 688 landslide fatalities, respectively (Seal et al., 2020). Importantly, we note that predicted fatalities for Costa Rica have large uncertainty bounds due to the highly variable nature of EQIL impacts, evident in the large spread of global fatality data in Figure 8.

Limitations of the Scenario Earthquake Method

Inherent uncertainties arise as a result of the scenario earthquake methodology utilized in this study, introduced into results primarily through three different sources, including the USGS ShakeMap system, the global statistical landslide model, and earthquake-induced landslide fatality model. The ShakeMap system uses an earthquake source radiation model, combined with a number of ground-motion prediction equations (GMPEs) to estimate the spatial distribution of ground shaking following an earthquake for a range of different tectonic and lithologic environments (Worden et al., 2020). The GMPEs applied here are based on global seismic wave attenuation relationships from similar tectonic environments, which do not take into account the particular three-dimensional wave propagation and site conditions present in Costa Rica. Alongside epistemic uncertainties contained within these GMPEs from estimation and simplification of various seismic attenuation relationships, uncertainty may also arise from the arbitrary selection of scenario earthquake hypocenter and fault rupture geometry. Distance from fault rupture is a critical parameter in all GMPEs used by ShakeMap; therefore, any deviation from the actual dimension, orientation, and location of a seismogenic fault introduces uncertainty into the distance term and resulting ground motion estimate. Resulting uncertainty in the distance parameter increases with magnitude, becoming the dominant uncertainty in ground motion estimates for large earthquakes above M 5.5 (Wald et al., 2008). One way in which we introduce this type of uncertainty is through our approximation of each fault as a rectangle, scaled in size by magnitude, where the earthquake hypocenter is located arbitrarily in the center of that rectangle. In reality, earthquakes seldom rupture in this idealistic manner, and unaccounted for effects from the non-rectangular rupture shape, as well as rupture directivity, may drastically change the actual distribution of ground shaking. For all events in this study, we try to reduce this uncertainty by obtaining earthquake fault geometries and locations from past studies and historical earthquakes. Nonetheless, any real event would undoubtedly differ in significant ways from our idealized source model. During an actual earthquake, near-surface conditions that are unaccounted for can have a significant impact on ground shaking, creating irreducible aleatory uncertainties in ShakeMap results (Worden et al., 2020).

Validation: Cinchona and Limon Earthquakes

In order to evaluate the performance of our model predictions for the two events with well documented landslide inventories (the Cinchona and Limon earthquakes), we test our predictions using three approaches: (1) a visual comparison of mapped landslide occurrence with model predictions (Fig. 7); (2) a statistical comparison of landslide occurrence with predicted landslide probability (Supplementary Materials, Table S4 and Fig. S7, S8, S9, S10); and (3) by examining prediction accuracy using a range of landslide probability thresholds (Supplementary Materials Fig. S8).

As discussed in the results and apparent in Figure 7, there is generally a strong spatial correlation between areas of high predicted landslide probability and mapped landslides. The statistical comparison shown in Supplementary Materials Table S1 and Figure S7 helps quantify this correlation. Our analysis demonstrates (1) that ~70% and 100% of landslides are located in areas with greater than 20% probability for Cinchona and Limon, respectively, and (2) the relative proportion of landslide-positive pixels systematically increases with probability. This suggests that our model appropriately accounts for landslide vulnerability for these two earthquakes, even if the absolute probabilities do not match the numbers of observed landslides. In order to formally test prediction accuracy, we initially apply a threshold of 50%, where any location with predicted probability 50% and above is considered likely to contain a landslide, and a location with probability below 50% is considered unlikely to contain a landslide. Each landslide inventory contains binary data (landslide positive or negative pixels); using this threshold, we calculate the number of landslide positive pixels correctly predicted by the model as: true positives (TP), incorrectly predicted as false positives (FP), correctly predicted landslide negative pixels (true negatives (TN)), and incorrectly predicted false negatives (FN). We find that there are a significant number of FN pixels when using the 50% threshold for both events. This suggests that this threshold may be arbitrarily high – in other words, that this threshold results in “missing” a significant number of true landslide events. We thus examine the statistical match between predictions and observations using a “balanced accuracy” criterion to identify the optimum threshold for each event, as defined in equation 3, where TPR is the true positive rate and TNR is the true negative rate.

(3)

We find that a lower probability threshold of 20% provides an optimum accuracy for the Cinchona event, while a significantly higher threshold of 55% is optimal for the Limon event (Supplementary Materials Fig. S8 and S10). At all thresholds, the balanced accuracy is significantly higher for the Cinchona event than for Limon. The significant over-prediction for the Limon earthquake may be attributed to a number of factors unaccounted for by our landslide model, including (1) complex lithologic, land-cover, or climatic variabilities and effects at a scale finer than the 0.25 – 1 km2 resolution used for model predictor variables, (2) differential weathering of geologic units resulting from the tropical climate of the Limon region, (3) numerous faults and folds that inhabit the compressional back-arc tectonic environment of Costa Rica, which can drastically alter the orientation of bedding (even within the same geologic unit), causing a unit to become more susceptible to failure if the bedding and topographic slope orientations coincide, and (4) the global landslide model being trained predominantly using earthquakes in high-relief mountainous terrains that may not apply reliably to the coastal, moderate-relief environment of Costa Rica’s Caribbean coast.

The global landslide model used in this study (Nowicki Jessee et al., 2018) was trained using 23 landslide inventories from historical earthquakes spanning a wide range of tectonic and climatic settings. Therefore, the relationships established between landslide occurrence and each landslide predictor variable are based unavoidably on averaged global observations that may not be entirely applicable to Costa Rica. Although our model predictions are demonstrated to be a reasonably good fit to the 1991 Limon and 2009 Cinchona coseismic landslide observations (Fig. 7), over-prediction is apparent for each event, especially the Limon earthquake. However, Nowicki Jessee et al. (2018) demonstrates that this magnitude of over-prediction is not an uncommon occurrence when the global model is applied to individual earthquakes. As shown in the supplementary materials (Fig. S8 and S10), the degree of interpreted under and over-prediction can be greatly altered by the probability threshold imposed for analysis. Note that it is not necessary to impose a probability threshold to use the model outputs, and the optimal probability threshold will depend on the relative level of conservativeness that may be appropriate for a given hazard planning or mitigation application. Uncertainty in model predictions may also be decreased through the use of more regional and locally applicable landslide proxies, higher resolution inventories, and incorporating various factors that decrease slope stability but are unaccounted-for in the current landslide model, such as anthropogenic alterations to landscapes (e.g., the construction of highway roadcuts).

It is crucial that these results are effectively disseminated from scientists to policymakers to ensure that necessary hazard mitigation measures are employed in susceptible regions. This study and data therein are made available to the Comisión Nacional de Prevención de Riesgos y Atención de Emergencias Costa Rica. Using the results of our study, the Comité Asesor Técnico de Geotécnia of this organization can more effectively bridge the gap between scientists, land planners, and policy makers to reduce risks related to seismicity and EQIL.

Limitations of the Landslide Fatality Model

The global landslide fatality model of Nowicki Jessee et al. (2020), which was used to calculate a range of potential fatalities for each scenario earthquake, was trained using observed landslide fatalities for 91 historical earthquake events, and integrates a proxy for vulnerability. While it is valuable as an order-of-magnitude estimate of the potential for casualties from earthquake-triggered landslides, we recognize that our fatality estimates include high levels of uncertainty, both due to the model’s calibration for global use, and intrinsic high variability in the relationship between exposure and number of resulting fatalities. High variability in casualty observations stems from the many unaccounted-for factors that control fatality occurrence, such as landslide type and time of day (Nowicki Jessee et al., 2020).

Uncertainty also results from the use of the United Nations Human Development Index (UNHDI), which generalizes the vulnerability of a population to a single value of economic development. Realistically, Costa Rica likely contains a diverse array of demographic groups with highly variable vulnerability levels. Due to these large sources of uncertainty and irreducible variability, it is impossible to precisely predict the exact number of fatalities resulting from any given hypothetical earthquake. Our fatality predictions are therefore solely order-of-magnitude estimates, and great caution should be used in the applications and context in which our results are presented. Despite these limitations and the inevitable imperfect match between our model predictions and landslide observations, the scenario approach outlined in this study represents a powerful tool to reduce natural hazard risk in road and city planning.

Utility and Future Applications of the Scenario Earthquake Method

Notwithstanding these limitations and uncertainties, we believe this study demonstrates the great potential of the scenario earthquake methodology for deterministic hazard assessment of earthquake secondary effects. Through our computation of five scenario earthquake suites across the range of tectonic environments present in Costa Rica, we suggest that large subduction megathrusts beneath the Nicoya and Osa peninsulas, and moderate-sized shallow intraplate earthquakes along faults within the central portion of the country, present the greatest potential for landslide occurrence and impacts to Costa Rica’s population. Our methodology solely utilizes globally available data, and therefore can easily be applied to any region of interest throughout the world. Perhaps the major utility of using scenario earthquakes for risk assessment is that it enables disaster mitigation and planning stakeholders to easily identify regions with the highest potential for EQIL exposure in a given area, which may then be combined with vulnerability estimates to better inform policy makers and aid organizations to prioritize use of risk reduction resources. Scenario hazard results may be especially valuable for regions with limited records of historical seismicity or lack of high-quality observational data. Though not investigated in this study, scenario earthquakes can also be readily used to assess impacts to critical infrastructure, including power stations, hospitals, and roadways, and can thus be used to identify specific areas that are most vulnerable to damage and would require resources following an earthquake. A number of the uncertainties detailed above may be reduced in future studies by using increasingly higher-resolution ground failure proxy data, or by incorporating locally or regionally calibrated earthquake-induced landslide and fatality models into our methodology.

Conclusions

We investigate the relative influence of landslide proxies such as topographic slope, peak ground velocity (PGV), and compound topographic index (CTI), together with earthquake source parameters such as magnitude and depth, on predicted landslide probability and fatality for five distinct tectonic environments in Costa Rica. These include subduction megathrust events beneath the northern Nicoya and southern Osa peninsulas, intraplate events in the Central Volcanic Range and Central Costa Rica Deformed Belt, and back-arc thrust events on the eastern coast of Costa Rica. We observe distinct threshold values of ~3° slope and ~ 8 cm/s PGV, and CTI values < ~6 in each subduction setting, while a combination of higher (~5°) slope and lower (~6 cm/s) PGV, and CTI (< ~5) in the CVR are required to produce landslide probabilities above 10%. In all tectonic environments, magnitude is the primary control on landslide occurrence, and is positively correlated with maximum landslide probability and overall areal landslide coverage. Depth plays a secondary but significant role, where deeper earthquakes decrease the maximum landslide probability while increasing overall areal landslide likelihood. We validate our global model results with local landslide inventories for two historical earthquakes, demonstrating a reasonably good fit between model predictions and observations, where over 70% of observed landslides are located in zones of 20% or greater landslide probability. Further validation suggests that the relative proportion of landslide-positive pixels systematically increases with predicted landslide probability, but that the model overpredicts landslide probability as a whole for both events. We also predict exposure and fatality ranges for each earthquake scenario, revealing that M 6.5+ events in the CVR and CCRDB pose the greatest hazard to population due to close proximity with many densely populated cities, while great (M ≥ 8.0) subduction megathrust earthquakes additionally pose high levels of hazard throughout the country. To further constrain our global landslide model predictions in the regional context of Costa Rica, future work will incorporate high-resolution geologic and topographic data, along with additional observed landslide inventories.

Acknowledgements

We thank Andrew Marden, Natalia Rodriguez, Joan Valverde, and Lea Vecchiarelli for their contribution digitizing mapped landslides from the 1991 Limon earthquake and assistance using mapping software. We are grateful to Elizabeth Sherrill for help with fault rupture geometry calculations and helpful suggestions regarding ways to analyze results. We also thank Kate Allstadt, Eric Thompson, Kristin Marano, and David Wald for their assistance with software, and continual support and insightful feedback at all stages of this project. Constructive reviews by Dr. Hernán E. Martínez were helpful in improving the manuscript content. Analysis was completed using Matlab (ver. R 2020b; MATLAB 2020). Maps were created using QGIS (ver. 3.14; QGIS 2020). Funding for this work was provided by the United States Geological Survey National Earthquake Hazards Reduction Program Grant No. G19AP00027 and the National Science Foundation Grant No. NSF OISE 1658648.

References

Adamek, S., Tajima, F., & Wiens, D. A. (1987). Seismic rupture associated with subduction of the Cocos Ridge. Tectonics, 6(6), 757-774.

Allen, T. I., Marano, K. D., Earle, P. S., & Wald, D. J. (2009). PAGER-CAT: A composite earthquake catalog for calibrating global fatality models. Seismological Research Letters, 80(1), 57-62.

Alonso‐Henar, J., Montero, W., Martínez‐Díaz, J. J., Álvarez‐Gómez, J. A., Insua‐Arévalo, J. M., & Rojas, W. (2013). The Aguacaliente Fault, source of the Cartago 1910 destructive earthquake (Costa Rica). Terra Nova, 25(5), 368-373.

Babayev, G., Ismail-Zadeh, A., & Mouël, J. L. L. (2010). Scenario-based earthquake hazard and risk assessment for Baku (Azerbaijan). Natural Hazards and Earth System Sciences, 10(12), 2697-2712.

Baker, J. W. (2008). An introduction to probabilistic seismic hazard analysis (PSHA). White paper version, 1, 72.

Bird, J. F., & Bommer, J. J. (2004). Earthquake losses due to ground failure. Engineering Geology, 75(2), 147-179.

Bozzano, F., Esposito, C., Martini, G., Martino, S., Prestininzi, A., Rinaldis, D., ... & Mugnozza, G. S. (2013). Earthquake-reactivated landslide scenarios in Southern Italy based on spectral-matching input analysis. Bulletin of Earthquake Engineering, 11(6), 1927-1948.

Budimir, M. E. A., Atkinson, P. M., & Lewis, H. G. (2014). Earthquake-and-landslide events are associated with more fatalities than earthquakes alone. Natural hazards, 72(2), 895-914.

Camacho, E., Hutton, W., & Pacheco, J. F. (2010). A new look at evidence for a Wadati–Benioff zone and active convergence at the north Panama deformed belt. Bulletin of the Seismological Society of America, 100(1), 343-348.

Carvajal-Soto, L. A., Ito, T., Protti, M., & Kimura, H. (2020). Earthquake potential in Costa Rica using three scenarios for the Central Costa Rica deformed belt as western boundary of the Panama microplate. Journal of South American Earth Sciences, 97, 102375.

DeMets, C., Gordon, R. G., & Argus, D. F. (2010). Geologically current plate motions. Geophysical Journal International, 181(1), 1-80.

Detweiler, S. T., & Wein, A. M.. (2017). The HayWired earthquake scenario. U.S. Geological Survey Scientific Investigations Report 2017–5013. https://doi.org/10.3133/sir20175013

Fernández, M., & Rojas, W. (2000). Faulting, shallow seismicity and seismic hazard analysis for the Costa Rican Central Valley. Soil dynamics and earthquake engineering, 20(1-4), 59-73.

Fisher, D. M., Gardner, T. W., & Marshall, J. S. (1994). Kinematics associated with late Cenozoic deformation in central Costa Rica: Western boundary of the Panama microplate. Geology, 22(3), 263-266.

GeoHazards International (GHI). (2018). Earthquake Scenarios, Nepal. https://www.geohaz.org/earthquake-scenarios-nepal

Harp, E. L., Wilson, R. C., & Wieczorek, G. F. (1981). Landslides from the February 4, 1976, Guatemala earthquake (No. 551.3 HAR). US Government Printing Office.

Hayes, G. P., Moore, G. L., Portner, D. E., Hearne, M., Flamme, H., Furtney, M., & Smoczyk, G. M. (2018). Slab2, a comprehensive subduction zone geometry model. Science, 362(6410), 58-61.

Hernández, G., Vahrson, W., & Ruiz, A. (1992). Deslizamientos producto del terremoto (4-22-91)/Landslides produced by the 22 Apr 1991 earthquake. Escuela de Ciencias Geográficas, Universidad Nacional de Costa Rica.

Jibson, R. W., Harp, E. L., & Michael, J. A. (2000). A method for producing digital probabilistic seismic landslide hazard maps. Engineering geology, 58(3-4), 271-289.

Keefer, D. K. (1984). Landslides caused by earthquakes. Geological Society of America Bulletin, 95(4), 406-421.

Konovalov, A., Gensiorovskiy, Y., Lobkina, V., Muzychenko, A., Stepnova, Y., Muzychenko, L., Stepnov, A., & Mikhalyov, M. (2019). Earthquake-induced landslide risk assessment: An example from Sakhalin Island, Russia. Geosciences, 9(7), 305.

Kritikos, T., Robinson, T. R., & Davies, T. R. (2015). Regional coseismic landslide hazard assessment without historical landslide inventories: A new approach. Journal of Geophysical Research: Earth Surface, 120(4), 711-729.

Kumar, S., Gupta, V., Kumar, P., & Sundriyal, Y. P. (2021). Coseismic landslide hazard assessment for the future scenario earthquakes in the Kumaun Himalaya, India. Bulletin of Engineering Geology and the Environment, 1-17.

Liu, J., Shi, J., Wang, T., & Wu, S. (2018). Seismic landslide hazard assessment in the Tianshui area, China, based on scenario earthquakes. Bulletin of Engineering Geology and the Environment, 77(3), 1263-1272.

Marano, K. D., Wald, D. J., & Allen, T. I. (2010). Global earthquake casualties due to secondary effects: a quantitative analysis for improving rapid loss analyses. Natural hazards, 52(2), 319-328.

Marshall, J. S., Fisher, D. M., & Gardner, T. W. (2000). Central Costa Rica deformed belt: Kinematics of diffuse faulting across the western Panama block. Tectonics, 19(3), 468-492.

Martino, S., Battaglia, S., D’alessandro, F., Della Seta, M., Esposito, C., Martini, G., Pallone, F., & Troiani, F. (2020). Earthquake-induced landslide scenarios for seismic microzonation: Application to the Accumoli area (Rieti, Italy). Bulletin of Earthquake Engineering, 18(12), 5655-5673.

Martino, S., Battaglia, S., Delgado, J., Esposito, C., Martini, G., & Missori, C. (2018). Probabilistic approach to provide scenarios of earthquake-induced slope failures (PARSIFAL) applied to the alcoy basin (South Spain). Geosciences, 8(2), 57.

MATLAB. (2020). 9.9.0.1570001 (R2020b). The MathWorks Inc.

Mendoza, C., & Nishenko, S. (1989). The north Panama earthquake of 7 September 1882: Evidence for active underthrusting. Bulletin of the Seismological Society of America, 79(4), 1264-1269.

Mora, S., & Mora, R. (1994). Los deslizamientos causados por el terremoto de Limón: Factores de control y comparación con otros eventos en Costa Rica. Revista Geológica de América Central, Volumen Especial Terremoto de Limón, 139-152.

Nayek, P. S., & Gade, M. (2021). Seismic landslide hazard assessment of central seismic gap region of Himalaya for a M w 8.5 scenario event. Acta Geophysica, 1-13.

Newmark, N. M. (1965). Effects of earthquakes on dams and embankments. Geotechnique, 15(2), 139-160.

Nowicki Jessee, M. A., Hamburger, M. W., Allstadt, K., Wald, D. J., Robeson, S. M., Tanyas, H., Hearne, M., & Thompson, E. M. (2018). A global empirical model for near‐real‐time assessment of seismically induced landslides. Journal of Geophysical Research: Earth Surface, 123(8), 1835-1859.

Nowicki Jessee, M. A., Hamburger, M. W., Ferrara, M. R., McLean, A., & FitzGerald, C. (2020). A global dataset and model of earthquake-induced landslide fatalities. Landslides, 17(6).

Protti, M., Gu, F., & McNally, K. (1994). The geometry of the Wadati-Benioff zone under southern Central America and its tectonic significance: Results from a high-resolution local seismographic network. Physics of the Earth and Planetary Interiors, 84(1-4), 271-287.

QGIS.org. (2020). QGIS Geographic Information System. QGIS Association. http://www.qgis.org

Quesada-Román, A. (2021). Landslide risk index map at the municipal scale for Costa Rica. International Journal of Disaster Risk Reduction, 56, 102144.

Rao, V. D., & Choudhury, D. (2021). Deterministic Seismic Hazard Analysis for the Northwestern Part of Haryana State, India, Considering Various Seismicity Levels. Pure and Applied Geophysics, 178(2), 449-464.

Robinson, T. R. (2020). Scenario ensemble modelling of possible future earthquake impacts in Bhutan. Natural hazards, 103(3), 3457-3478.

Robinson, T. R., Davies, T. R. H., Wilson, T. M., & Orchiston, C. (2016). Coseismic landsliding estimates for an Alpine Fault earthquake and the consequences for erosion of the Southern Alps, New Zealand. Geomorphology, 263, 71-86.

Robinson, T. R., Rosser, N. J., Densmore, A. L., Oven, K. J., Shrestha, S. N., & Guragain, R. (2018). Use of scenario ensembles for deriving seismic risk. Proceedings of the National Academy of Sciences, 115(41), E9532-E9541.

Rodriguez-Peces, M. J., Garcia-Mayordomo, J., Azanon, J. M., & Jabaloy, A. (2011). Regional hazard assessment of earthquake-triggered slope instabilities considering site effects and seismic scenarios in Lorca Basin (Spain). Environmental & Engineering Geoscience, 17(2), 183-196.

Ruiz, P., Carr, M. J., Alvarado, G. E., Soto, G. J., Mana, S., Feigenson, M. D., & Sáenz, L. F. (2019). Coseismic landslide susceptibility analysis using LiDAR data PGA attenuation and GIS: the case of Poás volcano, Costa Rica, Central America. In F. Tassi, O. Vaselli & R. Mora Amador (eds), Poás Volcano. Active Volcanoes of the World (pp. 79-118). Springer. https://doi.org/10.1007/978-3-319-02156-0_4

Salinas-Jasso, J. A., Ramos-Zuñiga, L. G., & Montalvo-Arrieta, J. C. (2019). Regional landslide hazard assessment from seismically induced displacements in Monterrey Metropolitan area, Northeastern Mexico. Bulletin of Engineering Geology and the Environment, 78(2), 1127-1141.

Seal D. M., Jessee A. N., Hamburger M. W., & Allstadt K. E. (2020). Comprehensive Global Database of Earthquake-Induced Landslide Events and Their Impacts. U.S. Geological Survey data release, https://doi.org/10.5066/P9NWIRZZ

Spence, R., So, E., Ameri, G., Akinci, A., Cocco, M., Cultrera, G., Franceschina, G., Pacor, F., Pessina, V., Lombardi, A. M., Zonno, G., Carvalho, A., Campos, A., Coelho, E., Pitilakis, K., Anastasiadis, A., Kakderi, K., Alexoudi, M., Ansal, A., Erdic, M., Tonuk, G., & Demircioglu, M. (2007). Earthquake disaster scenario prediction and loss modelling for urban areas. IUSS Press.

Suárez, G., Pardo, M., Domínguez, J., Ponce, L., Montero, W., Boschini, I., & Rojas, W. (1995). The Limón, Costa Rica earthquake of April 22, 1991: Back arc thrusting and collisional tectonics in a subduction environment. Tectonics, 14(2), 518-530.

Tajima, F., & Kikuchi, M. (1995). Tectonic implications of the seismic ruptures associated with the 1983 and 1991 Costa Rica earthquakes. Geologic and tectonic development of the Caribbean Plate boundary in southern Central America, 295, 327.

USGS. (2017). ANSS comprehensive earthquake catalog (ComCat). https://earthquake.usgs.gov/data/comcat/

von Huene, R., Ranero, C. R., Weinrebe, W., & Hinz, K. (2000). Quaternary convergent margin tectonics of Costa Rica, segmentation of the Cocos Plate, and Central American volcanism. Tectonics, 19(2), 314-334.

Worden, C. B., & Wald, D. J. (2016). ShakeMap Manual.Technical Manual, users guide. USGS.

Wyss, M. (2006, May). The Kashmir M7. 6 shock of 8 October 2005 calibrates estimates of losses in future Himalayan earthquakes. In B. Van de Walle & M. Turoff (eds), Proceedings of the 3 rd International ISCRAM Conference (pp. 1-5). ISCRAM.

Yue, H., Lay, T., Schwartz, S. Y., Rivera, L., Protti, M., Dixon, T. H., Owen, S., & Newman, A. V. (2013). The 5 September 2012 Nicoya, Costa Rica Mw 7.6 earthquake rupture process from joint inversion of high‐rate GPS, strong‐motion, and teleseismic P wave data and its relationship to adjacent plate boundary interface properties. Journal of Geophysical Research: Solid Earth, 118(10), 5453-5466.