1 Introduction

One of the main natural hazards that threaten human life and property is the landslides triggered by earthquakes (Kirschbaum et al. 2015; Lin et al. 2017). Central America is exposed to intense rainfall (3000–6000 mm/year), weathering, hydrothermal alteration, erosion, and fluctuating temperatures (Ruiz et al. 2019a). These conditions make this region vulnerable to multiple geohazards resulting from earthquakes, but mainly landslides (Bommer and Rodríguez 2002; Quesada-Román et al. 2019). In Costa Rica, there has been recurrent intermediate magnitude earthquakes near to volcanoes (5.0 < M < 7.7). Some of these earthquakes have triggered destructive coseismic landslides during the last decades with an average rate of one damaging earthquake every three years (Climent et al. 2008; Linkimer et al. 2018). Compared to other regions of the world, Central America is known to have a disproportional and high number of landslides triggered by earthquakes (Keefer 1984; Rodríguez et al. 1999; Bommer and Rodríguez 2002; Ruiz et al. 2019a).

The Central American landslide susceptibility is only comparable with some sectors of Southeast Asia and Oceania (Shi and Karsperson 2015). Some examples of comparable earthquakes are the 2011 Sikkim Mw 6.8 earthquake, the 2013 Lushan earthquake in Sichuan Ms 7.0 which caused 3883 landslides (Chen et al. 2020; Zhang et al. 2021), the 2015 Gorkha, Nepal Mw 7.8 earthquake which triggered more than 2064 landslides (Tian et al. 2019; Chen et al. 2020; Maharjan et al. 2021), and the relevant event of the 2008 Wenchuan, China Mw 7.9 earthquake which triggered more than 56,000 landslides (Dai et al. 2011; Li et al. 2020). Volcanoes in Costa Rica are in an active tectonic setting with a high seismicity rate (Alvarado 2021). In particular, the latest highly destructive event had a magnitude of 6.2 Mw Cinchona earthquake in 2009, which triggered approximately 4600 landslides (Barrantes et al. 2013).

There are several methods to model landslide susceptibility. The approaches include the weight-of-evidence analyses, index-based, data-overlay, neural network analysis, linear multivariate regression, boosted regression tree analysis highlight (e.g., Reichenbach et al. 2018; Arabameri et al. 2019; Chen et al. 2020; Li et al. 2020; Frigerio et al. 2021), discontinuous deformation analyses (DDA) (e.g., Chen et al. 2021; Zhang et al. 2021), remote sensing and high DEM resolution (e.g., Maharjan et al. 2021; Smail et al. 2021), and deep learning (e.g., Li et al. 2021). Two widely used landslide susceptibility methods in Latin America are the Mora-Vahrson-Mora (MVM) and the morphometric methods (Hengl and Reuter 2008; Barrantes et al. 2011; Quesada-Román and Barrantes 2016, 2017). The MVM method uses slope, lithological conditions, and humidity as control factors to susceptibility (Mora et al. 1992). In contrast, the morphometric method was summarized by Quesada-Román and Barrantes (2017) and uses the analysis of the slopes based on morphometric parameters such as depth and density dissection, relief energy, and total erosion, assigning weights and values for the susceptibility to occurrence of landslides. These are index-based or data overlap methods that allow the identification of areas with potential to manifest landslides and their triggers normally are heavy rains, earthquakes, or a combination of both. These methods have been used individually to evaluate landslide susceptibility near Poás volcano and, in this paper, we merge some of their variables to propose a more robust and accessible hybrid method.

During the last years, different researchers have studied the Cinchona earthquake and its subsequent landslides. First, a study using an airborne imagery mapped the landslides within the walls of the main active crater of Poás volcano and the geologic description of the area affected by landslides in route 126 at the eastern side of the Poás volcano (GVN 2009). In addition, a sedimentological analysis described the debris flows related to this seismic event, including a preliminary map of the landslides in the Poás massif (Alvarado 2010). Recently, other studies presented susceptibility models and landslides inventories (e.g., Barrantes et al. 2013; Quesada and Barrantes 2016; Ruiz et al. 2019a), and a national landslide risk index by municipality (Quesada-Román 2021). However, there is a lack of susceptibility evaluations as an integrated combination of complementary methods, which could help with interpretations in cases of limited information, especially for areas with scarce geological mapping.

We propose the evaluation of landslide susceptibility in the NW sector of the Poás volcano through the integration of the MVM method used by Ruiz et al. (2019a) and the morphometric parameters analysis. This research emphasizes coseismic landslides based on the evaluation of the seismic potential to propose the earthquake trigger for the modeling, looking for hypothetical large earthquakes. We seek to identify the worst-case scenarios that would generate landslides through a seismic analysis. Besides, we complement the seismic analysis with geomorphological mapping to define the final coseismic landslides zoning. It is intended to develop a landslide zoning that relates their greater or lesser incidence with the identified geomorphological units. Our approach can provide a reference for future disaster prevention, mitigation, and reconstruction zoning and a basis for landslide risk management in the NW sector of Poás volcano and/or other regions with analogous volcano-tectonic characteristics in tropical environments.

2 Study area

Costa Rica’s Central Volcanic Arc was formed by the subduction of the Cocos plate under the Caribbean plate at relative plate movements of ~ 8 cm per year (DeMets et al. 2010). This subduction has resulted in active volcanism and the faulting of Quaternary stratovolcanoes (Marshall et al. 2000). Volcanic activity in current volcanic massifs began in the Late Cenozoic, and it currently consists of fumarole emissions and hot lakes within the craters (Fernandez 2013). This mountain range is a chain of andesitic stratovolcanoes oriented to the northwest, consisting of five main massifs: Platanar, Poás, Barva, Irazú, Turrialba, and several pyroclastic cones (Alvarado 2011).

The study area (Fig. 1) is located at the NW of the Poás volcano, within an active fault belt that has been called the Central Costa Rica Deformed Belt (CCRDB). Regional tectonic stress and/or volcanic processes control the neotectonic structures surrounding the volcanic massif. Around the Poás volcano, the Viejo-Aguas Zarcas, Sabanilla, El Angel, and Carbonera faults are strike-slip faults-oriented NNW-SSE, while the San Miguel fault is a thrust fault with a general N70ºW orientation. Several authors mentioned that at least six historic shallow earthquakes of intermediate magnitude (5.0 ≤ M ≤ 6.5) have occurred near the Poás volcano in the last 242 years, which have generated destructive and deathly landslides in communities such as Sarchí, Grecia, Naranjo, Fraijanes, Bajos del Toro, and Vara Blanca (Montero 2001, 2014; Montero et al. 2010; Ruiz et al. 2014).

Fig. 1
figure 1

Map of the study area

Full size image

This area mainly comprises the upper Toro River catchment, which drains toward the Caribbean basin. The limit was defined according to the available geological information (Ruiz et al. 2019b) and represents the west sector of the area covered by LiDAR data and high-resolution orthophotographs after the Cinchona earthquake in 2009. A small proportion that is not covered by the LiDAR images has been complemented with the BID-cadaster program information with a resolution of 10 m in the Z axis (Epypsa 2011). The study area covers part of the mesoseismic area where the landslides triggered by the Cinchona earthquake occurred, the NW flank of the Poás volcano and part of the NE flank of the Platanar volcano. Most of the land use is covered by protected forest use; however, there are some important population settlements and agricultural sectors. In addition, there are important hydroelectric dams for electricity production in the country, and it is also an important tourist destination in Costa Rica, with one of the tops visited national parks: the Poás volcano National Park (Fig. 1; Pérez-Umaña et al. 2019).

The most important geological formations of the Poás volcano are the units Paleo Barva, Volcán Platanar, Andesites La Paz, Achiote, Lavas Río Cuarto, and the Volcán Poás Summit, composed mainly of andesitic-basaltic lava, lahars, ignimbritic, and pyroclastic deposits (Fig. 2). These geological formations, their lithology, age, and evolution are described by Alvarado et al. (2011) and Ruiz et al. (2019b) and have been used in this paper as reference to explore the influence of the lithology and the geological units age, with the susceptibility to present landslides. This effect has been previously described by Ruiz et al. (2019a), which associate higher susceptibility to older geological units.

Fig. 2
figure 2

Geological formations of the study area based on Ruiz et al. (2019a)

Full size image

3 Data and methods

3.1 Geomorphological mapping

We developed a geomorphological map at 1:25,000 scale to help understand the genesis, morphology, and chronology of the landscape of the NW of the Poás volcano and relate it to the analyses of the main seismic sources and finally to the zones prone to earthquakes. First, three digital elevation models (DEM) and a terrain analysis were created to obtain more detail of the geomorphological units. One of these DEMs have a resolution of 10 m, obtained from cadastral information (Epypsa 2011), the other one has a 20 m resolution, obtained from aerial photographs of the 2005 CARTA project, and the last one is from LiDAR images, obtained from the Costa Rican Institute of Electricity (ICE), with three points per m2 and the resulting DEM presents a resolution of 50 cm on the X and Y axes, and 15 cm on the Z axis. We analyzed the terrain models and chose the LiDAR imagery model to define the geomorphological units. These LiDAR images were acquired between March and April 2009.

The terrain models allowed the interpretation of detailed local geomorphological units, defined mainly by means of slopes, geometry, and river drainage network analyses, and using, as reference, the regional geomorphology of the Geomorphological Atlas of Costa Rica at a scale of 1:100,000 (Bergoeing and Brenes 2017). Finally, the landform, its geometry, morphometry, and its relation to the geological and tectonic context were interpreted to map the genesis, dynamics, morphology, and evolution of the different landforms and its processes using various digital graphic techniques to develop the final cartographic product (Gustavsson et al. 2006; Bishop et al. 2012; Otto et al. 2018). The resulting landforms were classified into fluvial, gravitational, structural, and volcanic genesis.

3.2 Seismic analysis

The maximum earthquake magnitude (Mmax) is defined as the largest possible earthquake for a seismogenic zone or region (e.g., Kijko 2004). The seismic sources were identified firstly based on previous neotectonics studies (Denyer et al. 2003; Montero et al. 2010) from which we selected the Viejo-Aguas Zarcas, Carbonera, Sabanilla, Ángel, and San Miguel as the faults that represent the greatest hazard to the study area. Then, the historical and instrumental seismicity of these faults were inspected according to the location of epicenters, magnitudes, and depth of the earthquakes. In addition, we inspected the traces of the tectonic structures and faults from aerial photographs and LiDAR data and two epicenters were chosen which asserts possible worst scenarios in the fault traces near to the interest area.

For the Mmax approximation, we used the seismic zoning proposed by Alvarado et al. (2017), specifically the crustal seismic zone of the Central Volcanic Range and the Costa Rica’s Central Valley (C6). Based on the seismological catalog of the National Seismological Network of Costa Rica (RSN) (Arroyo-Solózano and Linkimer 2021; Red Sismológica Nacional de Costa Rica 2021), between 1975 and 2020, and the portion of the historical catalog from Rojas et al. (1993), between 1723 and 1975, we used the Gutenberg–Richter (1944) recurrence law (logN = a−bM), in which N is the number of earthquakes greater or equal to magnitude M, and the a and b values are constants. The seismicity was modeled using a homogeneous catalog in Mw, the minimum magnitude was set to 3.0, the completeness analysis was taken according to the results from Arroyo et al. (2017) and using the Gardner and Knopoff (1974) declustering method. We based this analysis on 7825 seismic events from 1 to 40 km depth, determining the seismic parameters a and b-values by the Weichert (1980) maximum likelihood approximation. Moreover, we based the Mmax on the Makropoulos and Burton (1983) method, which is a pseudo-graphic formulation that approximates the Mmax by transforming the accumulated amount of the seismic moment into magnitude (e.g., Weatherill 2014). This approach was executed using the Hazard Modeler Toolkit (Weatherill 2014) of the OpenQuake software (GEM 2020).

We also evaluated the mean recurrence interval (MRI) of magnitudes above 5.0 Mw for the seismic zone C6. This is usually the time to elapse between two subsequent earthquakes of a given magnitude rupturing on a given seismogenic source (e.g., Wang 2007; Salazar 2018). This can be described by the Gutenberg–Richter relationship, if M is solved in terms of the frequency of the events (1/t) based on the estimated a and b-values. Besides, we estimated the accumulated probability of occurrence exceedance (POE) of the Mmax earthquake observed and inferred in the analysis for this seismic source, based on Poisson’s theory, defined as P(t) = 1−e−ƛt (Cárdenas et al. 2010; Rivas et al. 2014). This is approximated from the annual exceedance probability, which is very similar to ƛ (annual exceedance rate of earthquakes of a given magnitude). When graphing the Gutenberg-Richter relationship in terms of the earthquake frequency, ƛ is equal to N, obtaining N = 10a−b * M, getting the accumulated probability of occurrence for the magnitude of interest and evaluating these results for different periods.

3.3 Landslide susceptibility modeling

The model proposed by Ruiz et al. (2019a) is based on the MVM model (Mora and Vahrson 1991, 1993; Mora et al. 1992). This approach considers four control factors elements: (1) Lithological susceptibility (Sl), (2) Slope susceptibility (Sp), (3) Humidity content susceptibility (Sh), and (4) the landslide trigger event (DT), which uses the attenuation of the Peak Ground Acceleration (PGA) produced by an earthquake. This model has the advantage of specializing in earthquakes as landslide triggering events, and it can be applied to different shallow earthquakes with different characteristics in location, depth, and magnitude. Likewise, Quesada-Román and Barrantes (2017), based on morphometric parameters, used the dissection density (d), the dissection depth (p), the relief energy (I), and the total erosion (ET) to assess landslide susceptibility (SPL).

We propose to use the Ruiz’s model as a base but incorporating the morphometric parameters (M): dissection density (d) and dissection depth (p), replacing the lithological susceptibility parameter (Sl). This substitution has been defined to avoid giving a double weight to lithology since both the dissection and density and dissection depth have been observed in previous studies as closely related to the lithology of the geological units and their ages. In addition, this is suggested since the access to the drainage network and topographic models are commonly easier to obtain than geological detailed information in some areas. Equation 1 shows the resulting integration for the landslide susceptibility (SPL) proposed:

$${\text{SPL}} = \left( {{\text{Sp}}*{\text{Sh}}*M } \right)*{\text{DT}}$$
(1)

where \(M = \left( {d + p} \right)/2\).

3.3.1 Control factors

Landslide susceptibility comprises multiple factors. The quantity and accuracy of the control factors determine the quality of landslide susceptibility maps (Tien et al. 2016). We chose the slope, humidity content, and morphometric parameters as the most significant. To assess the slope susceptibility (Sp) (Fig. 3a), the generation of a digital terrain elevation model (DEM), explained in Sect. 3.1, and its subsequent calculation of slopes were required. The slope classes are based on expert criteria from the MVM method, which define seven slope classes (0; 0.1°–4°; 4°–8°; 8°–16°; 16°–35°; 35°–55° and > 55°), assigning a weight between 0 and 6, with 0 being the least susceptible and 6 the most susceptible.

Fig. 3
figure 3

Control factors weighted for landslides susceptibility. a Slope angle (Sp), b humidity content susceptibility (Sh), c Morphometric susceptibility (M): dissection density (d) + depth density (p)

Full size image

For the humidity content susceptibility (Sh) (Fig. 3b), the original estimation is obtained from a water balance based on monthly precipitation records. We propose to vary the calculation of this factor using the humidity provinces (Herrera 1985), specifically their water indices, assigning weights from 1 to 5 as follow: 1 = 0%; 2 = 0–20%; 3 = 20–100%, 4 = 100–300% and 5 = 300–600%). This percentage relates the mean annual precipitation (mm) (P) and the annual potential evapotranspiration (ETP), as Im = (P/ETP − 1) *100. These provinces are based on weather types for Costa Rica, determined by annual water balances, representing a mean of the humidity conditions. These kinds of indices are useful in areas where baseline information and meteorological stations are scarce and to consider altitude’s influence in the humidity regimen.

To include the morphometric parameters (M) (Fig. 3c), the rivers were first identified from the DEM obtained by LiDAR (Sect. 3.1). Subsequently, the area was divided into quadrants with individual surfaces of 0.1 km2, the centroid of each of these was obtained, and its respective ID assigned. The dissection density (d) was calculated obtaining the concentration of rivers per unit area, while the dissection depth (p) is explained as the height difference between the talweg and the maximum height of the analyzed quadrant. Consequently, the values for each quadrant were interpolated using the Kriging method and classified into five classes by the statistical method of quantiles. To combine the two morphometric variables (M), the range of values obtained for both depth and dissection density was normalized between 0 and 1 (Oyana 2021) and combined by means of the multiplication by a map algebra process. Finally, five classes were established using quantiles classification. This categorization represents the selected morphometric parameter (M), being consistent with the other susceptibility parameters established by Mora and Vahrson (1991, 1993), being 1 the least susceptible and 5 the most susceptible.

3.3.2 Trigger landslide event

The seismic trigger (DT) is determined from the evaluation of the seismic potential and hypothetical earthquakes. These earthquake scenarios were modeled as a simple fault by means of the Input Preparation Toolkit of the OpenQuake platform (GEM 2020) according to the epicenters of the hypothetical scenarios obtained from the analyses of the seismic potential explained in Sect. 3.2. The deterministic scenarios modeling using OpenQuake allowed to estimate the PGA values of the earthquake rupture definition: coordinates, magnitude, depth, and rake of the hypothetical earthquake and dip, upper, and lower depth of the fault source. Then, the Vs30 values of the soils are added to incorporate the site effects in the calculation. We used the Global Vs30 Mosaic of the United States Geological Survey (USGS), based on topographic slope, with custom embedded maps (Heath et al. 2020) and classified it in the soil classes according to the NEHRP (2020) classification. Finally, the PGA values were computed in OpenQuake engine, using a rupture mesh spacing of 0.2 km, a grid region spacing of 1 km and the Ground Model Prediction Equation (GMPE) of Zhao et al. (2006). This GMPE was selected according to Benito et al. (2012), who proposed this equation because it fits better with the attenuation observations for Central America’s crustal earthquakes.

Once the PGA values were obtained, they are converted to a weight between 1 and 7 according to the MVM method. For this process, we took the PGA values obtained from OpenQuake engine for each point on the grid of the study area and transformed it to the weight ranges for the trigger factor, described before. This was carried out by the empirical formula with a better logarithmic fit, proposed by Ruiz et al. (2019a): DT = 1,2597* ln (PGA)−1,2517.

3.3.3 Susceptibility maps and landslide zonation

According to the hypothetical earthquakes, based on a worse-case scenario and selected according to the procedures explained in Sect. 3.1, the different values of each control and trigger factors (Sects. 3.3.1 and 3.3.2) were integrated by means of a map algebra, obtained by the multiplication of each of them. Then, they were classified into five classes, by means of equal intervals, and the susceptibility to coseismic landslides was obtained. The classes used were very high, high, moderate, low, and very low.

To verify the effectiveness of the modeling, the Cinchona 2009 earthquake was additionally modeled following the same methodology proposed for the hypothetical earthquakes, as a real case of comparison and a model validation of the susceptibility classes. It was also modeled with OpenQuake engine as simple fault, based on El Angel Fault trace and according to its earthquake location and the rupture parameters based on Barquero (2009): 6.2 Mw, 4.6 km depth, rake −36.86°, dip 47.3° and upper and lower seismic depth of 0 and 15 km, respectively.

We combined previous landslide catalogs (Barrantes et al. 2013; Ruiz et al. 2019a), which correspond to the Cinchona earthquake. Each of these catalogs contain more than 4000 landslides, but we just evaluated the landslides in the eastern sector of the study area (available landslides information), using the riverbed of the Rio Toro as reference. This landslide catalog was revised according to the aerial photographs from the Carta 2005 project and the slopes were selected where landslides occurred, eliminating the deposits of them. This approach focuses the study on evaluating zones where landslides could occur, not slip volumes, or slip deposit zones. Moreover, some landslides were repeated in both catalogs, so a revision and selection of them was necessary to avoid duplications. The resulting catalog consists of 606 landslides and a total slip area of 2.37 km2. The verification was made by means of relationships between landslides areas and susceptibility areas obtained by the model using the Cinchona earthquake as trigger.

Finally, to represent the landslide susceptibility zoning, the susceptibility maps were averaged by means of a map algebra processing. This process was done to obtain a global view of the susceptibility for the entire area, that contemplates the possibility of a big earthquake, both, north and south of the study zone. The result was classified into three classes: high, moderate, and low using Natural breaks classification (Jiang 2013), to simplify the final global analysis and seeking to group the previous classes of very high and high, and the classes of very low and low susceptibility.

4 Results

4.1 Geomorphological units

4.1.1 Endogenic landforms

Tectonic-structural and depositional volcanic forms represent endogenic landforms controlled mainly by the volcano-tectonic fracture of the Poás volcano (Fig. 4). The most relevant structural morphology identified corresponds to the fault scarp (E1) of the San Miguel Fault, located in the northern sector. The most prominent expression is a scarp 100–300 m high, which extends 15 km and affects the drainage of the Sarapiquí River. Upstream of the fault scarp, the deep canyons are antecedents of this tectonic structure, which has its maximum height in the central part, decreasing toward the east and west. There are also a series of paleo-landslides along the scarp.

Fig. 4
figure 4

Geomorphological map and main faults of the NW sector of the Poás volcano. Costa Rica

Full size image

As volcanic morphologies, we determined the following: Caldera (V1), Volcanic crater (V2), Maar (V3), Secondary volcanic cones (V4), Flow fields and pyroclastic deposits (V5), and low (V6), moderate (V7) and high (V8) volcanic slopes. The volcanic crater (V2) is the active crater of the Poás Volcano, the Maars (V3) are Laguna Hule and Río Cuarto. In the latter, relevant concentrations of CO2 accumulated deep in the lagoon have been determined (Alvarado et al. 2011). The secondary volcanic cones (V4) are Cerro Congo and Von Frantzius, which responded to old sources of emissions that are not currently active. The pyroclastic flow fields and deposits (V5) are associated with areas covered by ash deposits from the Poás Volcano. Finally, the volcanic slopes reflect the mountainous relief produced by the volcanic deposits. The low slope (V6) responds mainly to more distant deposits and products of lahars or pyroclastic flows, while those moderate (V7) and high (V8) slope conform the volcanic buildings of Poás and Platanar, made up mainly of lava and pyroclastic deposits.

4.1.2 Exogenic landforms

Exogenic landforms are constituted by erosive and depositional landforms of fluvial and gravitational origin (Fig. 4). The fluvial genesis forms (F) respond directly to the action of the rivers, especially the Río Toro and some of its main tributaries. Two fluvial morphologies are identified, the floodplain (F1) of the Río Toro in the central sector and to the SW and the river canyons (F2) which cross the mountainous sector. The floodplains correspond to the settlement areas of some of the main towns of Bajos del Toro. The origin of these floodplains is related to river sectors that have not yet down-cut enough through the geological substrate to overcome the depositional processes that occur during flooding. The river canyons have depths between 50 and 100 m deep and are located mainly to the NW, the main one is the Río Toro canyon.

We also determined three gravitational morphologies: the erosion slopes (G1), the colluvial deposit (G2), and landslides scarps remnant (G3). These morphologies mostly respond to erosion processes on steep slopes and the accumulation of these deposits in flat or horizontal sectors. The erosion slopes, or foothills (G1), are mainly related to areas such as hills at the foot of the volcanic slopes, which have been shaped by erosive processes. Colluvial deposits (G2) are accumulation zones close to rivers and erosion slopes. Finally, ancient landslide zones (G3) respond to large landslides that have generated a specific geomorphology; these landslides are related to the Congo and Von Frantzius secondary volcanic cones.

4.2 Seismicity and seismic potential

Six important historical earthquakes have affected the study area. Three of these of 6.0 Mw occurred in 1772-02-15, 1888-12-30, and 1911-08-28, one of 5.5 Mw in 1912-06-06, one more of 6.1 Mw in 1955-09-01 in the Viejo-Aguas Zarcas fault and the Cinchona earthquake on 2009-01-09 (6.2 Mw) (Fig. 5a). Based on these historic earthquakes and according to the location of the epicenters, zones near faults that have not presented big ruptures yet and the trace faults analyses, we follow the deterministic criterion in terms of locating the epicenter within or as close as possible to the study area (worst scenarios) to select the trigger earthquakes (Fig. 5a). Scenario 1 consists of locating the epicenter at the Viejo-Aguas Zarcas fault, within the study area, related to Viejo-Aguas Zarcas fault. The second scenario proposes to locate it where an earthquake could occur in the San Miguel Fault, assuming a dip of ~ 30° for this fault, in the NW sector of Cerro Congo. The highest Mmax determined is 6.8 Mw and the shallow depth 5 km, based on the depth distribution of the earthquakes into the C6 seismic zone, with great influence of the Central Costa Rica Deform Belt (Fig. 5b).

Fig. 5
figure 5

a Tectonic framework of the study area and location of the trigger events used (hypotheticals and Cinchona earthquake 2009). 1: Hypothetical earthquake in Viejo-Agua Zarcas Fault. 2: Hypothetical earthquake in San Miguel Fault. The name of the main faults are shown in Fig. 2. b Central Volcanic Cordillera seismic zone and its seismicity between 1723 and 2018. The region contained within the dotted line represents the Central Costa Rica Deformed Belt (CCRDB). The dashed line represents the simplified northeast boundary of the Central American Forearc Block along the Volcanic Arc Faults (VAF). NPDB is the North Panama Deformed Belt

Full size image

The recurrence law parameters (a and b-values) estimated for C6 were 3.59 and 0.82, respectively, and the inferred and observed Mmax (Mw) were 6.4 and 6.8, respectively. The magnitude frequency distribution observed in these parameters suggests that the MRI in seismic zone C6 for the observed Mmax (6.4 Mw) is ∼ 70 years, and for the inferred Mmax (6.8 Mw), it is ∼ 110 years (Fig. 6a). The POE was also estimated for the observed and inferred Mmax, which resulted in a probability of occurrence greater than 80% for the observed Mmax (6.4 Mw) in ∼ 75 years and for the inferred Mmax (6.8 Mw) in ∼ 155 years (Fig. 6b).

Fig. 6
figure 6

a Mean recurrence interval (MRI) for earthquakes of certain Mw for the C6 seismic zone. The Mmax observed and inferred are marked. b Probability of occurrence (POE) for Mmax observed and inferred in seismic zone C6

Full size image

4.3 Susceptibility maps

Control factors (Fig. 3) showed for Sp the highest weights (3, 4, 5 and 6) correspond to scarps and river canyons, as well as mountain slopes and erosive cuts of the Poás and Platanar volcanoes, and low weights (1, 2) are mainly associated with plains and foothills. For Sh, the results showed three humidity provinces: humid, very humid, and excessively humid, with weights of 3, 4 and 5, respectively. This shows differences between mountainous areas of higher altitude in the Poás and Platanar volcanoes, with respect to flat areas and in the Río Toro valley. In the case of M, areas of lower weight (1) are associated with low slopes (< 15°), those of moderate weight (2 and 3) to the foothills, and those of greater weight (4 and 5) are associated to volcanic slopes, as well as the deep canyons of the Río Toro.

The Vs30 values for the study area shows a direct influence of the topography, with a range between 900 and 245 m/s, and classified in soils B, CB, C and D of the NEHRP (2020) classification (Fig. 7a). The derived site effect is observed in the PGA maps (Fig. 7b, c), it highlights the effect of high acceleration around the fault trace and the lower Vs30 values correlated with less attenuation of the seismic movement. Scenario 1 earthquake rupture was computed assuming a strike slip fault movement with a rake angle of 0° and a dip angle of 90° (Fig. 7b). Scenario 2 earthquake rupture was computed assuming a trust fault movement, with a rake angle of 90°, and a dip angle of 30° (Fig. 7c). For both hypothetical earthquake scenarios, the magnitude was defined as 6.8 Mw, a depth of 5 km, and an upper and lower seismogenic zone of 0 and 15 km, respectively. The results of the acceleration distribution show that PGA values greater than 500 cm/s2 could occur in a radius of ~ 10 km, close to the selected epicenters, with maximum intensities expected for both earthquakes of up to IX on the Modified Mercalli Intensity scale (MMI) (Fig. 7b, c).

Fig. 7
figure 7

a Soils classification and Vs30 values for the study area surroundings, based on topographic slope, with custom embedded maps (NEHRP 2020; Heath et al. 2020). b Hypothetical earthquake in Viejo-Agua Zarcas strike-slip Fault PGA values (scenario 1). c Hypothetical earthquake in San Miguel thrust Fault PGA values (scenario 2)

Full size image

The DT values and the susceptibility maps obtained for each of the hypothetical scenarios are shown in Fig. 8. For scenario 1, the distribution according to the modeling of susceptibility to landslide was very low = 36.1%, low = 27.8%, moderate = 22.7%, high = 12.1% and very high = 1.3%. The areas classified as high or very high susceptibility are located near the epicenter, specifically in the steep slope zone without vegetation located west of the main Poás crater, as well as on the southeast flank of the Platanar Volcano and the Toro River canyons (Fig. 8a). For scenario 2, the distribution of the susceptibility categories is as follows: for very low susceptibility 31.8%, low = 28.8% = , moderate = 23.4%, high = 13.3% and very high = 2.7% (Fig. 8b).

Fig. 8
figure 8

DT values for hypothetical earthquakes and coseismic landslides susceptibility maps obtained based on its combination with the control factors (Fig. 3). a Trigger event in Viejo-Agua Zarcas strike-slip Fault. b Trigger event in San Miguel thrust Fault

Full size image

4.4 Susceptibility model validation

The Cinchona earthquake landslide catalog (described in Sect. 3.3.3) was used to verify the model (Fig. 9a). We evaluated the total slipped area percentage, the relationship between the slipped area by susceptibility category and the category area (Fig. 9). Table 1 shows the distribution of landslides, part of the statistical and percentage data with respect to the susceptibility categories, obtained according to the final landslides catalog. The zones with the greatest slipped area were high, moderate, and very high susceptibility categories. The relationship between the slipped area and the susceptibility category area allowed us to determine the percentage of the area of each category that slides (Table 1 and Fig. 9b, c). Figure 9b shows very small landslide areas in low and very low susceptibility, mainly corresponding to sectors with low slope and far from the epicenter.

Fig. 9
figure 9

a DT values and coseismic landslides susceptibility map for Cinchona earthquake (2009) in the model validation area. b Area (km2) from coseismic landslides per susceptibility category. c Area from landslides per susceptibility category divided by the area covered by each susceptibility category

Full size image
Table 1 Landslides distribution based on susceptibility category in the eastern sector of the study area
Full size table

There is also a proportional increase in the slipped area as the susceptibility category rises, with a greater amount of slipped area in the high and very high susceptibility categories. Figure 9c and Table 1 show that large areas of very low susceptibility have, in proportion to their size, a very low percentage of landslide area (less than 2%), whereas the areas in the very high susceptibility category, in relation to its relatively small area, approximately 27% of its area exhibited landslides. The frequency of landslides shown in Table 1 demonstrated that more landslides were in moderate, high, and very high susceptibility, compared to low and very low, despite the fact that these last categories present much more area, which also supports the model effectivity to identify prone areas.

4.5 Landslides susceptibility zonation

Figure 10 shows three zones of high, moderate, and low hazard to coseismic landslides. These were defined by selecting the classification that best encompasses the modeled scenarios in a consistent manner. Low hazard zones (48%) are associated with slopes less than 20°, where the occurrence of significant landslides is not expected. Moderate hazard zones (35%) are characterized by slopes that generally exceed 20°, in sectors of the upper and middle Río Toro catchment. The high hazard zones (17%) are distributed near the upper catchments areas, in sectors of higher altitude and with steep slopes (> 30°). The foregoing is associated with pyroclastic pediments, canyoning river, and the slopes of the Poás and Platanar volcanic complexes. These conditions allow the development of landslides, mudflows, and complex movements.

Fig. 10
figure 10

Zoning of coseismic landslide susceptibility for the NW sector of the Poás Volcano, Costa Rica

Full size image

5 Discussion

5.1 Methodological and model implications

The seismic and geomorphological assessment determined that the faults with the greatest seismic potential are the San Miguel, Viejo-Aguas Zarcas, Sabanilla, and Carbonera fault. Some studies proposed the Venecia fault (Barquero 2009); however, in this investigation, according to Montero et al. (2010) and the geomorphological mapping, we interpreted it as a continuation of the El Angel Fault toward the NW (Fig. 2, 4 and 5a). The possible existence of a fault in the canyon of the Río Toro would also be important, nevertheless, there is not enough evidence to confirm its existence. We also identified and mapped a Caldera (V1) (Fig. 4) which had been previously reported in the Tectonic Atlas of Costa Rica (Denyer et al. 2003). Bergoeing and Brenes (2017) suggested that the Platanar-Porvenir Volcano could be part of a large collapsed caldera of up to 36 km in diameter, which highlights the importance of the hypothetical epicenter in the Viejo-Agua Zarcas fault.

The Mmax estimate was a 6.8 Mw, similar to the seismic potential associated to the C6 seismic zone and its faults in other studies (6.7–7.0 Mw) (Climent et al. 2006, 2008; Alvarado et al. 2017; Arroyo 2019; Arroyo-Solórzano and Linkimer 2021). The earthquake scenarios modeled have been defined based on the fault trace strike, modeling it as a simple fault, assuming a centroid moment tensor in the epicenter selected as a trigger point. This exhaustive selection of scenarios and the rupture modeled reflect the most probable natural behavior, with maximum PGA values along the fault traces. We consider that this source proximity effect has a very significant influence in the obtained values, which are very important to be considered in the seismic hazard approaches. However, this effect is not always determinant because large-scale landslides could occur in areas with low PGA, an example was observed during the 2016 Kumamoto earthquake in Japan (Chen et al. 2021).

The highest probability of occurrence for the modeled earthquakes is ~ 95% over 500 years. This model covers the most probable scenarios such as the observed Mmax of 6.4 Mw (Fig. 6). The used approximation method assumes that no earthquakes are occurring at the corresponding time. Furthermore, there is no consideration of other processes that contribute to the released energy in an aseismic way such as post-seismic or creep aseismic deformation. Despite the possibility that the method overestimates Mmax, it is widely used because it provides the upper limit of the seismic potential based on observed magnitudes.

The susceptibility maps have some key factors that improve the estimation from the previous MVM method. First, the Sh modification allows to differentiate humidity content susceptibility in different parts of the basin. This factor varies according to biotemperature and evapotranspiration. In high-altitude and mountainous areas, it is usually lower than low-altitude and flat areas. This condition favors a higher humidity content in the high-altitude areas, assigning greater Sh weight to the susceptibility of landslides in these sectors. It is important to indicate that this factor changes during the year due to the heavy raining season, but it works as an average annual value of humidity. Second, the morphometric parameters (M) respond to the presence of faults and fractures as well as to lithology, specifically the Paleo Poás Units (800–200 ka), which exhibit strong weathering, favoring the dissection depth and density. This allows to determine morphometric parameters that could substitute in a good way the geological maps in zones with scarce geological information. Finally, in the seismic analyses, we incorporated the seismic catalog recurrence evaluation, which was not applied in previous studies. Additionally, for the PGA expected values, the local effects of seismic amplification were also incorporated, based on Vs30 topographic values from Heath et al. (2020). It is very important to consider seismic amplification as a key factor for landslides susceptibility. We propose to improve this site effect map (Fig. 7a) with field geotechnical measures for a better characterization of the soils.

The susceptibility model shows an accurate and effective simulation (Fig. 9b, c). Despite the large areas of low and very low susceptibility, the slides zones were mostly in very high, high, and moderate susceptibility areas. According to Mora et al. (1992), landslides could occur in areas of moderate, high, and very high susceptibility. Consequently, 52% of the study area is prone to present landslides (moderate and high susceptibility zones, Fig. 10) and only the low-angle slopes would be safe. However, low susceptibility sectors with deep river valleys could develop mudflows/lahars due to the presence of weathered lavas and juvenile or recent pyroclastic deposits. If earthquakes such as those modeled occur, it could create hazardous situations due to landslides triggered, even in areas located more than 15 km from the epicenter. For both hypothetical events, the geological units of the Paleo-Poás or Platanar Volcano temporal phase show greater susceptibility, which agrees with the observations of Ruiz et al. (2019a).

5.2 Geomorphic and seismological approach for coseismic landslides hazard zonation

The baseline of the methodology developed follows the MVM method and its modifications over time (Mora and Vahrson 1993; Ruiz et al. 2019a). This method is similar to the logic applied in other studies worldwide (e.g., Tian et al. 2019; Shao et al. 2019; Chen et al. 2020). Other studies have carried out an inverse process, based on the evaluation of landslide catalogs, the higher incidence of these events and determining the possible trace of the fault responsible for the earthquake (Xu et al. 2013; Fan et al. 2018; Li et al. 2020, 2021). This methodology is more useful for areas with little knowledge of local tectonics or to complement seismic source analyses. In volcanic and tropical environments such as the Poás volcano, both methods are feasible and effective overcoming the difficulty of tracking surface fault rupture in the field due to dense vegetation. This limitation in tropical volcanic environments could be solved using unmanned aerial vehicles, even low-cost ones (Granados-Bolaños et al. 2021), radar satellite images, and remote sensing new technologies (Maharjan et al. 2021; Smail et al. 2021).

Many studies (e.g., van Westen et al. 2006; Fell et al. 2008; Gorum et al. 2013; Morell et al. 2020; Chunga et al. 2019; Chen et al. 2020) concluded that the integration of geomorphic and seismic data of sufficient quantity and quality into seismic hazard approaches greatly improves coseismic landslides hazard models. Our model improves previous MVM and morphometric models since it integrates and complements both methodology types and assigns a weight categorization according to the original MVM model. This is because geomorphological and seismic data have been exhaustively analyzed for the definition of the trigger, and some control factors have been adapted to obtain better detail (as with Sh) and substituting other parameters to make the evaluation more accessible (as with M). Multiple studies also indicated that a precise geomorphic assessment, morphometric data, humidity, soils, and detailed seismological data improve the prediction of events in tropical-volcanic environments (e.g., Dai et al. 2011; Lin et al. 2017; Otto et al. 2018; Fan et al. 2018; Ruiz et al. 2019a; Li et al. 2020; Arabameri et al. 2019; Tian et al. 2019; Chen et al. 2020; Maharjan et al. 2021; Smail et al. 2021; Zhang et al. 2021).

The main contribution of the proposed methodology is the detailed evaluation of the seismicity. This is accomplished through a seismological probabilistic analysis from the seismic catalog to determine the Mmax and complemented by a deterministic analysis based on earthquake scenarios obtained from the seismotectonic analyses and seismic potential assessment. This aspect has demonstrated its high relevancy in studies such as the ones of Xu et al. (2013), Morell et al. (2020), Shao et al. (2019), Chen et al. (2021) and Zhang et al. (2021) who evaluated the effects in areas with respect to the seismogenic fault and its intensity as a crucial parameter for spatial prediction of intensities and coseismic landslides. Our results show that the method is accurate for landslide susceptibility zoning and applicable to land use planning. We additionally consider that the greatest susceptibility for settlements and road infrastructure would be the transformation of these landslides into lahars, which was also mentioned by Ruiz et al. (2014).

Finally, we considered that the tropical context and the topography of volcanic environments such as in the NW Poás volcano provide a significant opportunity to keep developing methods to characterize the behavior of earthquake-trigger landslides in these sectors in a better way and to compare them with other contexts. An example of this is shown by Ruiz et al. (2019a), where the magnitude of landslides from the Cinchona earthquake was compared with another earthquake with similar characteristics in California, USA. They found a great difference in the number of landslides and a clear difference in the aforementioned factors, showing more coseismic landslides in tropical volcanic environments.

6 Conclusions

Based on the integration of the modified Mora Vahrson Mora (MVM) and morphometric methods, as well as the geomorphological and seismicity analysis, landslide susceptibility maps were determined for the NW sector of the Poás volcano in Costa Rica. Control factors were adapted, two of these being derivatives of the MVM (slope angle (Sp) and humidity content (Sh)), and two (dissection density and dissection depth) taken from the methodologies that evaluate the morphometric parameter, M. Moreover, the earthquake trigger (DT) was added based on the geomorphic and seismic analysis. Therefore, worst-case earthquake scenarios are proposed as two hypothetical earthquakes. The 2009 Cinchona earthquake was also modeled as an earthquake scenario and their susceptibility categories were validated through the landslide catalog of this earthquake, showing a very good fit between landslides and the more susceptible categories. Then, the two scenarios were integrated and a susceptibility zoning map to coseismic landslides was obtained.

Our results improved the local seismic and landslides knowledge due to the control factor detail implemented in the model. More precise and detailed information was obtained using the Sp and M factors by the digital elevation model based on the LiDAR images, as well as the Sh, by means of a water index. In addition, 14 local geomorphological units have been mapped to show differences between endogenic and exogenic landforms and considered for the seismotectonic analyses and the selection of the location of the epicenters in the earthquake scenarios. Furthermore, the seismic potential (Mmax) results in a 6.8 Mw and a depth of 5 km for the hypothetical earthquakes. These earthquakes have been modeled: one with an epicenter associated with the Viejo-Aguas Zarcas fault (scenario 1) and the other with the San Miguel fault (scenario 2). The validation of the method with the 2009 Cinchona earthquake showed an adequate correspondence between categories and real landslides with most of the landslides found in the moderate, high, and very high susceptibility categories predicted by our model. The susceptibility zoning map to coseismic landslides indicated that 52% of the area has a high landslides probability. Conical volcanic landforms and those on slopes greater than 15° are the most prone morphologies to landslides, showing that the slope is the most relevant control factor.

The results also serve as the basis for generating similar studies in other parts not only of Costa Rica, but of the world, in volcanic and seismically active regions with limited geological data. The main virtues of the proposed method are the improved detail in some control factors and the use of easily obtainable morphometric parameters to replace lithology, as well as the integration of a geomorphic and seismic analysis. The methodology is recommended as an input for risk mitigation and territorial planning and based on these models and a computational development, early warning systems (EWS) can be implemented in real time to prioritize higher-risk areas of attention to disasters. It is also recommended to integrate these results with probabilistic landslides hazard models in the future and complement them with the assessment of landslides triggered by heavy rains. Coseismic landslide susceptibility zonation is the first step to integrate the different elements at risk such as the hazard, exposition, and vulnerability.