Determination of Primary and Secondary Lahar Flow Paths of the Fuego Volcano (Guatemala) Using Morphometric Parameters

Abstract

On 3 June 2018, a strong eruption of the Fuego volcano in Guatemala produced a dense cloud of 10-km-high volcanic ash and destructive pyroclastic flows that caused nearly 200 deaths and huge economic losses in the region. Subsequently, due to heavy rains, destructive secondary lahars were produced, which were not plotted on the hazard maps using the LAHAR Z software. In this work we propose to complement the mapping of this type of lahars using remote-sensing (Differential Interferometry, DINSAR) in Sentinel images 1A and 2A, to locate areas of deformation of the relief on the flanks of the volcano, areas that are possibly origin of these lahars. To determine the trajectory of the lahars, parameters and morphological indices were analyzed with the software System for Automated Geoscientific Analysis (SAGA). The parameters and morphological indices used were the accumulation of flow (FCC), the topographic wetness index (TWI), the length-magnitude factor of the slope (LS). Finally, a slope stability analysis was performed using the Shallow Landslide Susceptibility software (SHALSTAB) based on the Mohr–Coulomb theory and its parameters: internal soil saturation degree and effective precipitation, parameters required to destabilize a hillside. In this case, the application of this complementary methodology provided a more accurate response of the areas destroyed by primary and secondary lahars in the vicinity of the volcano.

1. Introduction

Primary lahars are flows which are formed as a direct result of a volcanic eruption. They tend to be bulky (107–109 m³) and record high speeds (>20 m/s). Their maximum flows are commonly between 103–105 m³/s. These features provide the ability to flow long distances, even hundreds of kilometers downstream. They occur primarily when during an eruptive event incandescent material causes the fast melting of large volumes of ice and snow of the glaciers that cover some volcanic edifices and flows by descent gullies. Secondary lahars mainly include lahars caused by the rains. The unbound, by previous eruptions, pyroclastic material can be easily removed by the rains. In general, these are at lower speed, volume and they travel shorter distances as compared to primary lahars, however, they are most frequent during periods of rain [1].

Guatemala is the only country in the Central American region that has trained local observers to generate reports on volcanic activity, who then transmit the information via radio and/or telephone three times a day. The information is transmitted to the National Institute of Seismology, Volcanology, Meteorology and Hydrology (INSIVUMEH) plant when the following characteristics are present: changes in the release of energy, increase in the number of explosions, increase in the expulsion of ash, increase in seismic activity, rumblings, shock waves, and the manifestation of block avalanches that descend through the ravines in the volcanic perimeter.

On 3 June 2018, a strong eruption of the Fuego volcano in Guatemala produced a dense volcanic ash cloud rising to 10 km high. Following the collapse of the volcanic column, pyroclastic flows and descending lahars caused significant damage and a high number of fatalities around the mountain [2]. After the eruption, INSIVUMEH, with support from the Volcano Disaster Assistance Program (VDAP) of the United States Geological Survey (USGS), the University of Edinburgh, and Michigan Technological University, elaborated 2 scenarios of lahars for medium and heavy rains based on the numerical models obtained from the program LAHAR Z (Figure 1). These scenarios can be seen on the MapAction page [3].

On 19 November and for the fifth time during 2018, the Fuego volcano erupted, with a 1200-m long lava flow that descended into the Ceniza canyon. This volcanic flow with fragments of rock producws by erosion on the slopes of the volcano moved downhill, incorporating enough water, so that they formed mud flows and volcanic debris or lahars that descended the slopes of the volcano, affecting sectors such as Las Lajas, El Jute, besides the Ceniza and Seca ravines, according to INSIVUMEH in its last special volcanological bulletins (Nos. 157-2018 and 158-2018) [4,5].

The theory of LAHAR Z’s operation will not be studied in this article, since it is not the objective of this study. Its operation is reviewed briefly. In this article, reference is made to the LAHAR Z program and its results as a basis for superimposing the results obtained by the Differential Interferometry Synthetic Aperture Radar Differential (DINSAR) analysis, the morphometric indices and the unstable areas near the volcano.

LAHAR Z was written to delimit areas of potential lahar inundation from one or more user-specified lahar volumes. LAHAR Z is a code that is executed within a geographic information system (GIS), which was created by the USGS [4] and is based on a semi-empirical model that delimits flood risk zones by lahar (that is, areas drawn to represent the probable floods during a lahar event) in a digital elevation model (DEM).

The program uses two semi-empirical equations calibrated by statistical analysis in the cross section of an area flooded by a lahar (A) and the flooded planimetric area (B) measured in 27 lahars deposits of 9 volcanoes in the United States of America, Mexico, Colombia, Canada and Philippines [6] (Figure 2).

The equations are:

$$A={\alpha }_1{V}^{2/3}$$

(1)

$$B={\alpha }_2{V}^{2/3}$$

(2)

where A is the maximum section area flooded, B is the total area flooded and V is the volume of the lahar, ${\alpha }_1=0.05$ and ${\alpha }_2=$ 200, are constant values.

The LAHAR Z model works under assumptions such as:

• Floods from past lahars can provide the information basis to predict possible future floods.
• Distant hazard zones are confined to the valleys and are directed to the slopes of the volcano.
• The volume of the proximal lahar controls the extent of the distal flood forming downstream.
• Very voluminous lahars occur less frequently.
• No one can predict the size of the next lahar that will descend by a river.
• The volume selection for the lahar and the origin site thereof control the flow range.

As a complementary mechanism to the results obtained by DINSAR and Lahar-Z, the morphometric model has been used to determine the trajectories of the destructive lahars with greater precision by combining morphometric parameters obtained from the digital elevation model of 12 m of spatial resolution such as the flow accumulation (FCC), the topographic wetness index (TWI), the length-magnitude factor of slope (LS), the degree of internal saturation of the soil (h/z), where z is soil depth, h is water level above the failure plane and the effective precipitation (q/T) required to destabilize a slope, where h is the height of the water table and z is the thickness of the colluvium that slides above the failure plane. All these parameters are obtained from the shallow landslide slope stability model (SHALSTAB) to complement the mapping of this type of lahars using remote sensing (differential interferometry, DINSAR) in Sentinel images 1A and 2A, to locate areas of deformation of the relief on the flanks of the volcano, areas that are possibly origin of these lahars. To determine the trajectory of the lahars, parameters and morphological indices were analyzed with the software System for Automated Geoscientific Analysis (SAGA). The parameters and morphological indices used were the accumulation of flow (FCC), the TWI, the length-magnitude factor of the slope (LS). Finally, a slope stability analysis was performed using the Shallow Landslide Susceptibility software (SHALSTAB) based on the Mohr–Coulomb theory and its parameters: internal soil saturation degree (h/z) and effective precipitation (q/T), parameters required to destabilize a hillside. In this case, the application of this complementary methodology provided a more accurate response of the areas destroyed by primary and secondary lahars in the vicinity of the volcano [6] (Figure 3). These parameters are described in more detail in the Materials and Methods chapter.

Although the limits of the areas delimited by the numerical modeling of lahars must be taken with caution and are considered as references and not as absolutes, the comparison between the limits of hazards by flows of lahars obtained from modeling with LAHAR Z and the morphometric model is mainly differentiated by their theoretical conceptualization.

Based on the damage caused by the pyroclastic flows, it was found that the volcanic hazard maps and physical vulnerability contained errors for two main reasons: a lack of field observation data to adjust the limits of the hazard map simulated with LAHAR Z, and the poor tracing of possible secondary lahars and the analysis of their distal trajectories originating in proximal sources [5].

2. Materials and Methods

The Sentinel 1 and 2 images obtained from the Alaska Satellite Facility [7], were used to determine the deformations of the terrain using the DINSAR technique [8] and those from Sentinel 2 were used to determine the changes in the morphology of the relief, which in this case was the delimitation of primary and secondary lahars, and compare them with the results obtained from the morphological model. The digital elevation model of 12 m of spatial resolution was obtained from the Copernico Project [9]. The planimetric information such as the results of the LAHAR Z software application was provided by the National Coordinator for Disaster Reduction (CONRED).

To complement the results obtained with the DINSAR analysis and the morphological method, in this study, the spatial distribution of potential sources of primary and secondary lahars during eruptions and torrential rains was analyzed using the slope stability model of shallow landslides. SHALSTAB is a plugin of the SAGA GIS.

The general flow of processes used in this study can be seen in Figure 4.

The processes A, B and C are described below:

Process (A): DINSAR analysis. The Sentinel 1 and 2 images were obtained from the Alaska Satellite Facility [9]. Sentinel 1 images were used to determine the deformations of the terrain using the differential SAR interferometry (DINSAR) technique [10] and those from Sentinel 2 were used to determine the changes in the morphology of the relief, which in this case was the delimitation of primary and secondary lahars and compare them with the results obtained from the morphological model. Process B: the planimetric information such as the results of the LAHAR Z software application was provided by CONRED. Process C: digital elevation model of 12 m of spatial resolution was obtained from the Copernico Project. The field observations carried out by CONRED technicians and sent in the format of points in shp format, helped to calibrate and verify the results obtained from the presented methodology, as well as to understand which of the three morphological parameters used had greater weight according to the scenario in the one that developed the lahars.

For the study of the primary and secondary trajectories of lahars using the DINSAR technique, the morphometric method and the stability analysis, the following parameters were analyzed:

The relief deformation analysis by the Sentinel Application Platform (SNAP) [11] followed the workflow indicated in Figure 4 from the Sentinel 1B images taken on May 21 and June 14, 2018.

This analysis was carried out to determine the possible zones of deformation of the ground where the possible secondary lahars start.

Brief synthesis of the processes to obtain the Interferogram:

• Base images: Sentinel 1B
  • S1B_IW_SLC__1SDV_20180521T001317_20180521T001345_011012_0142C5_A797
  • S1B_IW_SLC__1SDV_20180614T001319_20180614T001347_011362_014DD5_D06F;
• Mission: SENTINEL-1B
• Pass: ASCENDING
• Polarization: VV

Process diagram in ESA-SNAP [10] (Figure 5).

Step 1

Read-Topsar-Split-Orbit file: Open the products (Figure 5A).

• S1B_IW_SLC__1SDV_20180521T001317_20180521T001345_011012_0142C5_A797
• S1B_IW_SLC__1SDV_20180614T001319_20180614T001347_011362_014DD5_D06F;

View the products. (Figure 5B): In the Sentinel-1 IW SLC products, there are 3 subsets IW1, IW2, IW3. Each subset is for an adjacent strip of satellite data collection by TOPS mode. The Fuego volcano can be seen in IW3 (Figure 5C).

Step 2

View a band geocoding-enhanced spectral diversity. To view the data, Intensity_IW3_VV band.

Coregister the images into a stack: For interferometric processing, two or more images must be coregistered into a stack. One image is selected as the master and the other images are the slaves. The pixels in slave images will be moved to align with the master image to sub-pixel accuracy (Figure 5D).

Step 3

Form the Interferogram-Topsar Deburst-TopoPhase Removal-Filtering-Sarsim Terrain Corection-Write. (Figure 5E): Interferogram formation from the InSAR products menu.

The current location of the primary and secondary tabs after the eruption of 3 June 2018, were obtained with a Sentinel 2 image of 10 m of spatial resolution where bands 4, 3 and 2 were combined to obtain the current position of the lahars (Figure 6). A brief synthesis of the processes to obtain the real position of lahars is presented below:

• Base Image: Sentinel-2B: 2018-10-22T16: 23: 19.024Z
• Pass: DESCENDING
• Bands: 4-3-2

Process diagram in ESA-SNAP

Step 1

Read-Open the products. Granule-img-bands (Figure 6A).

View the products. (Figure 6B): In the Sentine-2B-Granule-Img products, there are 13 spectral bands in the visible and near-infrared (VNIR) and short-wavelength infrared (SWIR) spectrum, as show in the below table (Figure 6C).

Step 2

Combine bands 4-3-2 RGB natural colors. Location of primary and secondary lahars (Figure 6D).

Process (D): mass movements analysis. To evaluate the susceptibility for landslides and define potential areas of sediment production through mass removal processes, the SHALSTAB program was used, which is a model to map potential shallow landslides. SHALSTAB combines the characteristics of the stationary subsurface flow (water flow that moves at shallow depths below the surface of the land with a constant velocity), considering the morphology-hydrology relationship [12]. SHALSTAB uses two types of numerical models: the slope stability model and the hydrologic model.

The slope stability model is based on an infinite slope form of the Mohr–Coulomb failure law in which the downslope component of the weight of the soil just at failure, τ, is equal to the strength of resistance caused by cohesion (soil cohesion and/or root strength), C, and by frictional resistance due to the effective normal stress on the failure plane Equation (3).

$$\tau =C+\left(\sigma -\mu \right)\tan \phi$$

(3)

where $\sigma$ is the normal stress, $\mu$ is the pore pressure opposing the normal load and $\tan \phi$ is the angle of internal friction of the soil mass at the failure plane. This model assumes, therefore, that the resistances to movement along the sides and ends of the landslide are not significant. As SHALSTAB’s goal is the mapping of landslide hazards, then establishing zero cohesion maximizes the degree of instability possible at the study site. By eliminating cohesion, Equation (4) can be written as:

$${\rho }_{s}gz\cos \theta \sin \theta =\left[{\rho }_{s}gzco{s}^2\theta -{\rho }_{w}ghco{s}^2\theta \right]\tan \phi$$

(4)

where z is soil depth, h is water level above the failure plane, cρs and ${\rho }_w$ is the soil and water bulk density, respectively; g is gravitational acceleration and $\theta$ is slope angle (Figure 7A). This Equation (5) can then be solved for h/z which is the proportion of the soil column that is saturated at instability:

$$\frac{h}{z}=\frac{{\rho }_s}{{\rho }_w}\left[1-\frac{\tan \theta }{\tan \phi }\right]$$

(5)

h/z could vary from zero (when the slope is as steep as the friction angle) to $\frac{{\rho }_s}{{\rho }_w}$ when the slope is flat (tan θ = 0). An important assumption, however, will be used below which sets a limit on what h/z can be. It is assumed that the failure plane and the shallow subsurface flow is parallel to hillslope, in which case h/z can only be less than or equal to 1.0 and any site requiring h/z greater than 1 is unconditionally stable—no extreme rain can cause it to fail. Figure 7B illustrates the relationship between h/z and $\tan \theta$ for an angle of internal friction φ of 45° and a bulk density ratio of 1.6. There may be 4 stability states: any slopes equal to or greater than the friction angle will cause the right-hand side of (3) go to zero, hence the site is unstable even if the site is dry (h/z = 0). This scenario is “unconditionally unstable”, commonly corresponds to sites of bedrock outcrop. Because h/z cannot exceed 1.0 in this model, if $\tan \theta$ is less than or equal to tanφ(1-(ρs-ρw) then the slope is “unconditionally stable”. The two other stability states are “stable” and “unstable”, with the former corresponding to the condition in which h/z is greater than or equal to that needed to cause instability (given by the right hand side of Equation (5)) and the latter corresponding to the case in which h/z is less than that needed to cause instability.

Hydrologic model

To model the hydrologic controls on h/z. SHALSTAB use a steady state shallow subsurface flow based on the work by [14,15]. Assume that the steady state hydrologic response model mimics what the relative spatial pattern of wetness (h/z) would be during an intense natural storm which is not in steady state. This assumption would break down if precipitation events are sufficiently intense that thin soils on non-convergent sites can quickly reach destabilizing values of h/z before shallow subsurface flow can converge on unchanneled valleys.

Figure 8 illustrates the geometry and routing of water off the landscape used in the SHALSTAB hydrologic model. If we assume that there is no overland flow, no significant deep drainage, and no significant flow in the bedrock, then q, the effective precipitation (rainfall minus evapotranspiration) times the upslope drainage area, a, must be the amount of runoff that occurs through a grid cell of width b under steady state conditions. Using Darcy’s law, we can write that (6):

$$qa=k_{s}h\cos \theta \sin \theta b$$

(6)

where $\sin \theta b$ is the head gradient, $k_s$ is saturated hydraulic conductivity. At saturation the shallow subsurface flow will equal the transmissivity, T, (the vertical integral of the saturated conductivity) times the head gradient, sinq and the width of the outflow boundary, b and this we can approximate as follows (7):

$$\frac{h}{z}=\frac{q}{T} WI;WI=\frac{a}{b sen \theta }$$

(7)

where $h$ represents the height of the water table on the sliding surface (m), $z$ the depth of the soil (m), $q$ effective precipitation [mm], $T$ soil transmissivity [m²·day-1], $WI$ (Wetness Index) factor describes the effect of morphology on the subsurface flow (this is calculated in SAGA as the topogrraphic wetness index (TWI)).

In order to delimit areas of danger, new studies have applied a methodology with a geographical information system (GIS) base, identifying zones where flows start along the hydrographic network [12].

Morphologically, the triggering of a detritic or lahar flow depends on two factors: the slope of the channel and the critical flow. The speed of an elementary gravity wave is called the critical velocity. When the velocity of the surface water in a channel equals the critical velocity, the elementary gravity waves with upstream motion remain at the point where they were generated and the flow in the channel is called critical.

When the superficial velocity of the water in the channel is equal to the critical velocity, the elementary waves of gravity with movement upstream remain at the point where they were generated and the flow in the channel is called critical. Both factors can be obtained from the morphological characteristics of the basin. The first one is evaluated directly from the digital terrain model, while the liquid flow is evaluated indirectly through the drained area [12].

Studies conducted in the Swiss Alps found a relationship for the critical slope (S) that triggered detrital flows and the drained area of the basin (A) (Equation (8)).

$$S=c A^{-n}$$

(8)

The coefficients c and n have values of 0.32 and 0.20, respectively. The inverse relationship between S and A indicates that a greater slope of the channel will require a smaller liquid flow to trigger a debris flow (a smaller amount of liquid within the flow area).

Other authors have proposed classifying lahar flows into four categories: triggering, propagation, deceleration, and deposit [13].

Triggering occurs on slopes greater than the threshold proposed by other authors; 25° is known as a slope of minimum probable slope for the initiation of debris flow [13,14,15], but with a lower angle than the internal friction of the ground and a drained area smaller than 10 km². The maximum slope corresponding to the internal friction angle of the deposit is such that, in areas with steeper slopes, the amount of detritus is generally modest or negligible, while for drained areas exceeding 10 km², the type of transport becomes fluvial in place of detrital. The model calculates the degree of the internal saturation of the soil with the relationships of the height h/soil depth of the groundwater level z; (h/z) (Equation (9)) is the necessary factors to destabilize the slopes.

$$\frac{h}{z}=\frac{\rho s}{\rho w}(1-\frac{tan\theta }{tan\varnothing })$$

(9)

This formula allows two limit states to be defined. When h/z is negative, the slope is unstable for any degree of saturation and is called “unconditionally unstable”.

On the other hand, when h/z is greater than one, it is called “unconditionally stable” because even in cases of saturation, the slope is stable. To define the intermediate situations, it is necessary to add the hydrological model.

The hydrological model considers a subsurface flow in a permanent regime.

Applying the law of Darcy, on the characteristics of the movement of water through a porous medium, it is possible to arrive at the following expression that links the degree of saturation (h/z) with effective precipitation (q), which is the rainwater that enters the soil and is assimilated by the plants (Equation (10)).

$$\frac{h}{z}=\frac{q}{T} WI;WI=\frac{a}{b sen \theta }$$

(10)

With

• $q$ = effective precipitation [mm];
• T = soil transmissivity [m²·day⁻¹];
• θ = slope [°];
• a = contribution area [m²]; and
• b = is the length of the contour through which the flow transits [m].

where T is the transmissivity of the soil; a is the drained area; and b is the dimension of the raster cell. The relation q/T represents the magnitude of the precipitation relative to the capacity of the soil to conduct the water. The WI (wetness index) factor describes the effect of ground morphology on the subsurface flow.

The analysis of the stability on the slopes of the volcano with SHALSTAB through the hydrological relationship q/T, the magnitude of precipitation, q, in relation to the capacity of the subsoil to transmit water along the T (transmissivity) slopes determined that the greater q was in relation to T, the more likely the soil would be saturated, and the greater the number of sites with unstable slopes. This relationship used the classification of pixels with stable positive values to quasi stable and chronically unstable values with negative values according to the stability classes (

Table 1. Stability classes defined by SHALSTAB and modified from Montgomery and Dietrich [15].

Classes SHALSTAB	Interpretation of Class
Chronic instability	Unconditionally unstable and unsaturated
Log q·T -1 < -3,1	Unconditionally unstable and unsaturated
3,1 < Log q·T -1 < -2,8	Unstable and saturated
-2,8< Log q·T-1 < -2,5	Unstable and unsaturated
-2,5< Log q·T-1 < -2,2	Stable and unsaturated
Log q·T -1 > -2,2	Unconditionally stable and unsaturated
Stable	Unconditionally stable and saturated

).

Process E: The morphometric method with SAGA. Once the areas of surface deformation of the relief were found with DINSAR (Sentinel 1 images), the current location of the lahars (Sentinel 2 images), the superficial flow accumulators of the runoff and instability zones (SHALSTAB) were determined, the morphological method was applied to complement the location of trajectories of lahar flows with the analysis of the following morphological index.

2.1. LS-Factor (Slope Length and Steepness Factor)

This is a combination of the slope and the length of the slope as when combined, it is a useful attribute to predict erosion potential (Equation (11)). The length factor of the slope, L, calculates the effect of the length of the slope on the erosion, and the slope factor of the slope, S, calculates the effect of the slope on the erosion [16,17]. The higher values of brown color represent a greater susceptibility to erosion (Figure 9).

LS = L^(1/2)/100*(1.36 + 0.97 + 0.1385*S^2)

(11)

where L is the length of the slope in meters, and S is the gradient of the slope in %.

2.2. Topographic Wetness Index (TWI)

The topographic wetness index models the dynamics of surface and subsurface flows based on the topographic control of runoff, which offers a better perspective in relation to the prediction of sites where the saturation and high concentration of runoff can act as initial flow paths to processes greater than flood (Figure 10).

The TWI combines the contribution to runoff from a local area drained and the pending of it, and it is commonly used to quantify topographic control on hydrological processes and is defined as [6,7,8,9,10,18,19,20] (Equation (12):

TWI = Ln (af/tan β)

(12)

where af is the local area drained for a calculation point, and tan β is the directional slope of the cell of interest (and of the eight neighbors in the case of using an algorithm D8. Ln is natural logarithm.

2.3. Flow Accumulation (FCC)

Flow accumulation is calculated as the upslope contributing (catchment) area using the D8 model or multiple flow direction approach [17,18,19]. The result is an accumulated flow raster for each cell, determined by the accumulation of the weight of all of the cells that flow into each cell of the descending slope. The accumulated flow is based on the total number of cells, or a fraction of them, flowing to each cell in the output raster.

The output cells with a high flow accumulation are areas of concentrated flow and can be used to identify watercourse channels (Figure 11).

3. Results

The differential synthetic aperture radar interferometry (DINSAR) analysis, based on the Sentinel 1 images, showed the lack of magmatic deformation before the eruption, but there were surface subsidences on the associated slopes due to the passage and deposits of lahars (Figure 12). This lack of magmatic deformation corroborates the studies undertaken by several researchers in the volcanoes of the area [20,21,22], which attribute the absence of deformation to differences in the storage of magma in relation to other continental arcs. The absence of magmatic deformation before the lahars can mean that the intrusion of magma is not the trigger of the lahars. Most likely, the rain is the cause.

The analysis of stability on the slopes of the volcano with SHALSTAB determined that the zones near the volcano were the most unstable with a qualification of quasi stable to chronic instability through the erosion of the channels where the lahars had flowed and where there had been greater destruction as is the case of El Rodeo, San Miguel, Monte María, and sectors near the Guacalate River (Figure 12B,C).

The current locations of the lahars produced by the eruption of 3 June 2018 and later were delimited with the Sentinel 2 images and their destructive power was verified with field data mainly from La Réunion, San José Las Lajas, Santa Clara Las Lajas, San Miguel, and El Rodeo, amongst other towns (Figure 12C). Figure 13A shows the sectors before the eruption. The implementation of lead-colored lahars on that date. In line of red color, the limits of the model of lahars obtained with LAHAR Z. In blue are the ravines and canals where the primary and secondary lahars flowed that flooded the sectors, and in green-yellow-ocher color the areas of ridges and crests. The field control points of the destroyed sites are represented by red crosses. In (B) Sectors close to La Réunion, in (C) sectors near San Miguel and in (D) close to Monte María.

From the analysis of the morphometric parameters LS, TWI, and FlowAcc, the areas with greater susceptibility to erosion were obtained as well as those with the highest concentration and accumulation of superficial flow of rainwater that were directly related to the primary and secondary lahars (Figure 14). These areas of blue color were observed to have been directed by the lahars that destroyed the populated areas. These three morphometric indices [23,24,25,26] are valid parameters to determine the flow paths of primary and secondary lahars. The areas with high values of these indices are areas of moisture concentration in the slopes, flow paths and concentration of sediments that can move from the high parts by intense rains and gravity generating secondary lahars, generally on the slopes.

The field control was carried out by technicians of CONRED.

4. Discussion

Unlike the LAHARZ model, the presented model uses morphometry as a value that changes the trajectory of the primary and secondary lahars, which allowed us to obtain the real flow paths of these phenomena and their effects. A beneficial coincidence of the application of both models is that they kept a central axis path of the main lahars, due to the fact that they were drains or main rivers. However, LAHARZ does not separate secondary lahars, due to this action; the morphological method can determine the flow paths of the mentioned phenomenon.

The lahar danger map, created for the Fuego volcano, based on the LAHARZ program, published in 2018 by CONRED, established high hazards zones (marked in red color) that can be affected by lahar flow paths. Torrential rains have generated secondary lahars that have affected the surrounding areas of the volcano, even those ones of low hazard. The medium hazard zones were affected by lahar flow paths from the La Reunión site to the Alsacia site, as can be evidenced in Figure 13.

In the current paper, with the proposed methodology, a study of the lahars flow paths has been carried out, being the main study areas the valleys of these Barranca Seca, Ceniza, Guacalate, Playa Trinidad, Las Cañas, El Jute and Las Lajas rivers. Lahars originated due to the erosion of the cause ways of detritus-ash-water deposited during the volcanic eruption and torrential rains. The axis of the primary lahar paths studied with the proposed methodology are like the ones obtained to delimit the high hazards zones in the mentioned hazard map.

On the other side, is has been proved that the modeled secondary lahar deposits with the proposed methodology destroyed civilian structures in La Reunión, San José de Las Lajas, Santa Clara de Las Lajas, San Miguel and other areas had location limits in accordance with the data obtained with LAHARZ and the location obtained with Sentinel 2A images, but with an added element that was not obtained with LAHARZ which is the outline of the secondary lahars flow paths, some as contributors to primary lahars and others as independent elements with no connection to the main ones.

The range of the secondary lahars and its respective modeled central axis is similar and bigger to the June and November lahars of 2018, observed in the Sentinel images. This is probably due to the resolution of the digital elevation model and the data given by the analyzed morphometric index.

Finally, a comparison of the affected areas has been carried out and its reach on transit zones, obtained with the proposed methodology and LAHARZ, observing that the entire area affected by lahars is within the high danger zone. The lahars corresponding to the Las Lajas y Guacalate rivers, take the biggest amount of area of the affected areas, with most of the civilian infrastructure destroyed as a result, turning the mentioned areas to maximum danger zones.

The surface terrain deformations, obtained from Sentinel 1B, located in the high mass erosion and removal areas can be part of secondary lahars which were generated from high ground. This is a contribution when the DINSAR analysis is used.

In addition to this work, further research studies can improve the location of the lahars flow paths using the deformed areas of the relief obtained with the DINSAR analysis, the stability criterion of slopes and analysis of morphometric parameters in digital models of high elevation spatial resolution obtained from Lidar images and integrating them into specialized programs in lahars and floods analysis such as LAHAR Z, Titan 2D, Ash3D, Flo-2D, IRIC, Iber, among others. The final objective of these investigations is to put in the hands of government authorities and the community, vulnerability maps to this type of threats to reduce them.

5. Conclusions

As a contribution of the application of this methodology, based on a flow modeler such as LAHAR Z or another program, to know the Flow Paths of the primary and secondary lahars that have caused the destruction in the vicinity of the volcano, this methodology deepened the study of these trajectories based on the analysis of the morphometric parameters, which allowed us to have a better knowledge of the possible paths before a new volcanic eruption. This methodology is currently being applied in the Fuego volcano by INSIVUMEH technicians with the objective of putting it in the hands of the government authorities, as well as the institutions linked to the subject, which allows these conditions to be ascertained in advance thereby implementing structural measures and non-structural measures to reduce this danger.

To complement the results obtained with the DINSAR analysis and the morphological method, in this study, the spatial distribution of the sources and potential routes of primary and secondary lahars during the eruptions and torrential rains that can erode and produce mass movements in areas of critical instability was obtained from the stability model of shallow landslide slopes with SHALSTAB.

The integration of the DINSAR methodologies to determine the surface deformations of the ground, where the SHALSTAB model determines the stability-erosion of the slopes of the volcano, and the morphometric method to analyze the primary and secondary lahar currents, is a procedure that can help to prevent the destructive effects that they cause during these events as demonstrated in the same destroyed localities that were not mapped in the maps of existing hazards by CONRED.

The field observations carried out by CONRED technicians and sent in the format of points in shp format, helped to calibrate and verify the results obtained from the presented methodology, as well as in understanding which of the three morphological parameters used had greater weight according to the scenario in the one that developed the lahars. This parameter corresponds to the TWI, the same one that determined that the lower areas of the streams are the ones that get wet and fill quickly overflowing and flooding nearby areas.

The morphological model could obtain the primary and secondary lahar paths in all sectors near the volcano in a very short time, unlike the numerical models such as the LAHAR Z, which are applied in a single axis of a current or drainage. In this case, the paths of the lahars obtained are sufficient to produce maps of the vulnerability reduction of civil structures and socio-economic systems including the preservation of human lives.

This combination of procedures can be applied in any volcanic landscape based on digital models of high-resolution elevation, obtaining very good agreement with the paths of the real flows.