UAV-Based Quantification of Dynamic Lahar Channel Morphology at Volcán de Fuego, Guatemala

Abstract

This study quantified erosional and depositional processes for secondary lahars in Las Lajas drainage at Volcán de Fuego, Guatemala, during the rainy season from May to October 2021. Abundant pyroclastic material from ongoing eruptive activity is remobilized seasonally during heavy precipitation, which can impact infrastructure and populations living near Fuego. Our region of focus was in an agricultural zone 6 to 10 km from the summit, surveyed with an unoccupied aerial vehicle (UAV) quadcopter at monthly intervals. Imagery was processed into overlapping time-lapse structure from motion digital elevation models (DEMs). DEMs were differenced to find volumetric changes as a function of the channel flow path distance (quantified in 500 m sections) to track channel morphology changes over time. The largest measured volume changes were a 490 m³/day loss in the upper section (~6 km from summit) and a 440 m³/day gain in the lower sections (~10 km from summit). We discussed how the natural channel’s constriction and widening of Las Lajas in more distal sections control the behavior and stability of the stream evolution. Above the constriction, the channel is primarily downcutting and meandering within an old flood plain, which had been filled in by pyroclastic materials deposited by the June 2018 paroxysm.

1. Introduction

Volcán de Fuego (14.473°N, 90.88°W) is a 3800 m open vent stratovolcano with 3000 m of local relief located in central Guatemala. Fuego actively produces pyroclastic material during small Strombolian eruptions with varying frequency up to tens of times per hour [1] depositing material ranging from ash-sized to bombs and blocks on the upper slopes. The current style of activity, including occasional larger paroxysms producing pyroclastic currents (PCs) has been continuous since May of 1999, but Fuego has been intermittently erupting for centuries. Fuego’s usual mode of Strombolian activity is punctuated with larger paroxysmal events including 50 eruptions with a Volcanic Explosivity Index (VEI) ≥ 2 since 1524 [2]. These larger events produce PCs [3], which fill drainages radiating from the volcano with hot debris easily mobilized during the rainy season. Our study site, Las Lajas drainage, was filled in during a series of PCs in June of 2018 including one flow, which jumped the bank of the channel and devastated the village of San Miguel Los Lotes, killing many and burying the structures in hot rock [3]. The estimated volume deposited in this series of PCs was 15.1 ± 4.2 × 10⁶ m³ [4]. The majority of the material originated from the bulking and subsequent erosion of pyroclastic material at the top of Las Lajas near the vent. This led to increased mobility, allowing the PC to travel much farther than earlier events in the series, reaching 12 km from the summit [4].

Lahars are slurries of mud and debris common at Fuego and other volcanoes. They are analogous to debris flows except that they occur on the flanks of volcanoes and transport primarily volcanic material. Lahars account for about 25% of the volcano-related deaths since 1500 AD, second only to pyroclastic currents (PCs) [5]. Most infamously, a primary lahar, induced by an eruption in 1985 at Nevado del Ruiz, Colombia, resulted in 28,000 deaths in the City of Armero [6]. Smaller secondary lahars, induced by rainfall rather than eruption, are a much more common and persistent hazard at many volcanoes. Lahars are particularly dangerous because they do not require an eruption trigger. Many years after eruptive activity has ended, secondary lahars still present hazards controlled by the local weather. In October of 1998, for instance, Casita volcano in Nicaragua, which has no historic record of eruptions, was hit by a hurricane that caused an avalanche and lahars which killed an estimated 1600 people [7].

Secondary lahars are generally characterized by a bouldery snout followed by a watery tail [5,8]. These events may be very common due to erosional activity caused by rainfall remobilizing pyroclastic materials [9,10]. Lahars at Volcán de Fuego, featured in this study, commonly incise deep channels that can be over 100 m deep in places with steep and unstable walls. Lahars can quickly erode such walls and bulk up creating hazards that are extremely dangerous and destructive near to a drainage as well as in fluvial fans far from the volcano [10]. They can be an especially significant threat to populated communities, local infrastructure, and property situated adjacent to these channels. Lahars are capable of forcing boulders and logs against bridges or dams with impact forces of 10 to 1000 tons per square meter [7]. Channel paths are often meandering and dynamic and can erode into farmland and create steep cliffs. Flooding is another concern when lahars pass beyond their confining channels. Structures may be inundated and destroyed due to a shifting or reactivated channel. Road destruction and temporary closures are a common occurrence at Fuego when the activity is high because most vehicular transport crosses channels without bridges. Even small lahars are highly erosive and can undercut banks and bridges. Lahars can also be harmful to infrastructure through deposition. They “freeze” when the flow energy is dissipated enough that the head cannot be pushed forward anymore and the water seeps away, leaving a deposit of condensed sediment and boulders [5]. If that deposit is on or around infrastructure or agricultural lands, removal can be costly. Lahar deposition may divert streams, widen flood plains, and shift meanders forcing land use adaptation [3,11].

Unoccupied aerial vehicles (UAVs) are growing in popularity for use in natural hazard disaster management due to their ever-increasing accessibility and ease of use [12]. Although they are limited in areal coverage and surveys may be impacted by weather conditions, they are versatile, low-cost tools for collecting high-resolution data with short notice [13]. Trends in research concerning rapid onset natural hazards utilizing UAVs under 25 kg have been on the rise since 2014, when the number of journal articles went from 9 or less per year to over 180 by 2020 [12]. UAV use in the geosciences is expanding outside of natural disaster research as well. Many areas of research that have classically involved heavy field mapping with boots on the ground can be streamlined with these tools [12,13,14,15,16,17,18,19,20,21].

Given the rapid morphological changes present in active lahar channels UAV-acquired structure from motion (SfM) repeat 3D mapping of the terrain is a useful and practical tool for the study of both lahar erosion and deposition. SfM has been used to map geomorphic change over time in many studies ranging from volcanic craters to river morphology [14,15,16,19,21]. UAV-acquired data enhanced utility over satellite data for a focused study area because the cost of operation for producing high-resolution maps (spatially and temporally) is relatively low. Such focused mapping can complement higher spatial coverage and lower-resolution satellite-derived topography. Inexpensive, consumer-grade UAVs have been proven to produce quality SfM products at sub-meter resolution [14,15,16,19,20,21]. This permits high accuracy DEM of Difference (DoDs) models. Our study showed the capabilities of SfM DoDs for quantifying the erosional and depositional regimes of a lahar channel over the course of a rainy season.

2. Materials and Methods

2.1. Study Area

This study focused on a particularly active lahar drainage at the active Volcán de Fuego in Guatemala. The region-of-interest (ROI) was a ~3.5 km long zone, with a total surveyed area of approximately 1.6 km², of the Las Lajas lahar channel on the southeastern flanks of the volcano (Figure 1). This portion of the channel was chosen because it is relatively accessible for repeat field surveys and it regularly experiences lahars that can be monitored by time-lapse cameras and other ground-based geophysical sensors. The ROI was also part of the path impacted by a series of PCs in June of 2018, which lasted 16 h and filled the lahar channel with tens of meters thick deposits. The main phase of the eruption lasted 5 h and destroyed the village of San Miguel Los Lotes and a section of La Reunion Golf Resort [3,10]. Monitoring instrumentation belonging to the Instituto Nacional de Sismología, Vulcanologia, Meteorología e Hidrología (INSIVUMEH) was situated at various locations along the ROI and this permitted timely public notices to be published when a secondary lahar was in progress. Erosion deeply cut the sediments in the drainage and provided a dynamic environment for study. Our four-month study coincided with the Guatemalan rainy season (May through October) when secondary lahar activity was common, occurring, on average, about once per week in Las Lajas.

2.2. Data Collection

SfM made use of a series of overlapping photographs taken in each survey zone of the ROI to build high-resolution DEMs and orthomosaics. Our surveys were made in four zones (labeled A–D; Figure 1) and repeated at semi-regular monthly intervals (36–48 days apart) throughout the 2021 rainy season.

Table 1. Details of the flights made during the 2021 rainy season from May to October. Each survey is roughly separated by a month. Each flight record contains a 6-digit date, YYMMDD, and the number of photos taken during the flight. Red cells indicate no data.

	Survey 1	Survey 2	Survey 3	Survey 4
Zone A	210607 590	210720 258	210902 258	211008
Zone B	210602 279	210720 265	210902 265	211008 270
Zone C	210602 286	210720 283	210831 282	211008 283
Zone D	210603 382	210720	210901 213	211012 279

summarizes the 14 flights; surveys 2D and 4A were not successful and were not used in the analysis. Preprogrammed flight paths were flown with a consumer-grade UAV, the DJI Air 2, which has a ½ inch Complementary Metal Oxide Semiconductor (CMOS) RGB sensor with a 24 mm equivalent focal length capable of an effective 48 MP photo [22]. The flight planning software, Dronelink version 2.41, used a fixed altitude set at a constant 100 m height above the take-off point. This software did not allow for terrain tracking, and so, surveyed regions with a range of topography resulted in maps with variable resolution. The flight pattern was set to the normal ‘lawnmower’ style of acquisition, with an average of 260 photos taken per survey zone and a forward and side photo overlap of 80% at 100 m above the takeoff point (Figure 2, Table 1). The typical aerial extent of each survey zone was 0.5 km². The area covered by each survey was primarily controlled by balancing the battery limitations of the UAV against acquiring the best data possible for reconstruction in the available work time before the afternoon rain began.

The blue lines show the planned path for the UAV. The white lettered points A-E are used to adjust the size and shape of the polygon in the app and the white + signs are to add more points to the polygon. The green marker is the start point at 100 m above the takeoff point and the red marker is the end of the planned flight.

2.3. Data Processing

2.3.1. Agisoft Metashape Pro

Initial data processing was carried out with Agisoft Metashape Pro version 1.7.3 on the highest resolution settings to produce the DEMs and orthomosaics at a sub-decimeter resolution with an average ground sample interval of 5.22 cm/pix. The workflow used in Metashape was adapted from the USGS guidelines (Figure A1) in Over et al. [20]. Parallel processing was carried out across 5 Alienware computers with Intel(R) Core (TM) i9-9900K processors and NVIDIA GeForce GTX 1660 Ti graphics cards. Iterative gradual selection and camera optimization were used to improve the quality of the tie points before the dense point cloud was populated [20]. Definition of a local coordinate system was accomplished using synthetic ground control points (SGCPs) chosen by picking stationary features that exist in all the surveys. Real ground control points (GCPs) were not possible for this study due to access and equipment limitations; however, SGCPs allow registration of the DEMs to each other spatially and can correct any warp present from inaccuracies in the GPS elevation metadata collected by the UAV. The first and last surveys’ orthomosaics were processed for direct comparison to establish static features. Ten to fifteen SGCPs dispersed through each zone were chosen to guarantee good coverage as recommended by Agisoft [23]. We recognized that while SGCPs are not ideal for absolute georegistration, they are appropriate for intra-image registration and DEM differencing. They rely upon one set of coordinates held as a control from the non-georegistered survey, which is determined by the onboard UAV GPS. The Mavic Air 2 accuracy is claimed to be ±0.5 m horizontally and ±1.5 m vertically [22]. In Agisoft, the camera accuracy was set to the default 10 m, and the SGCP accuracy was set to the default 0.005 m. Our primary goal was to identify changes in channel morphology over time and this does not require absolute georegistration relative to a global coordinate system but rather a local coordinate system. Toward this objective, DEMs and orthomosaics were exported in the TIFF file format with 10 cm/pixel resolution and fixed, consistent geographic bounds permitting efficient comparison and differencing. In total, 28 products (14 DEMs and 14 orthophotos) were produced and are provided in a linked database.

2.3.2. QGIS-Based Processing

Co-registration of images is impacted by the growth of the plants in the study area and required consideration during processing. To address this influence, the orthomosaics for the final survey, when the foliage would presumably be at its largest, were imported into QGIS 3.22.11, and the semi-automatic classification plugin (SCP) was used to mask the foliage. The program was trained by selecting polygons containing the target classifications, namely, rock versus foliage. The resultant classified regions were then sieved to reduce the miscategorized rock shadows posing as plants. Edges were still visible on the plant classification, so they were dilated by 2 pixels to add a buffer around the plants. These products were then exported as a TIFF file to create a foliage mask for the DEMs.

2.3.3. CloudCompare Co-Registration

The DEMs in several survey zones possessed long-period slope artifacts that also required correction during co-registration. The point clouds from Metashape were, thus, imported into CloudCompare version 2.13. alpha and aligned in all directions using manual alignment and the Align tool with an average root mean squared error (RMSE) of 1.64 cm. Each point cloud was then exported as a raster TIFF file with a 10 cm pixel resolution. These DEMs were then imported into MATLAB version 9.6 R2019a to match the pixel dimensions of the other DEMs of the same zone to permit vertical feature alignment for static points.

2.3.4. MATLAB DEM Calculations

The DEMs were imported to MATLAB and processing was performed on areas outside of the vegetation mask made in QGIS. The individual surveys of a zone were then subtracted from one another to find the DEM of Difference (DoD)(Equation (1)) [15].

ΔZ₂₁ = Z(x,y)_t2 − Z(x,y)_t1,

(1)

Z(x,y)_t2 in Equation (1) indicates the vertical relief in each x and y coordinate in time 2 (t2). From the DoD, we could then calculate volume change (ΔV_t21) in Equation (2) as the integrated elevation difference ΔZ multiplied by 0.01 m² to convert from pixels to decimeters. Since the DEMs were exported with 10 cm/pix resolution, the conversion from pixels to meters is a simple unit conversion.

$${\mathsf{\Delta }\mathrm{V}}_{\mathrm{t}21}={\sum }_{ }^{Section}{\mathsf{\Delta }\mathrm{Z}}_{\mathrm{t}21} \times 0.01 {\mathrm{m}}^2$$

(2)

Since the Mavic Air 2 logs distance above the take-off point (rather than absolute elevation), the z direction was unconstrained. To limit the z offset error, the zones were adjusted to a local coordinate system by manually zeroing out known static areas between DEMs, such as roads and roofs. The DoDs in the static areas then showed a zero value which persisted everywhere there was no change detected. The DEM and differenced DEMs were used to create profiles showing erosion (or aggradation) over time in each zone. Each DoD was cropped to only include the drainage and a small portion of the surrounding banks in order to quantify the volume of material moved in each interval.

Volume calculation errors are interpolated across the bridge, which obscures visibility to the ground in survey zones B and C. To rectify this, in the preprocessing, the tie points for the bridge structure have been removed, effectively cropping the structure [21]. Any remaining artifacts were masked in a later step, leaving a void that will be filled with NaNs in MATLAB for DoD drainage volume calculations.

The survey zones were purposefully flown with large areas of overlapping terrain that needed to be accounted for in the volume sum. These overlaps were merged by averaging DEM values. The full mosaic of 4 zones (Figure 1) was split into 500 m longitudinal sections according to the stream path distance within the channel (Figure 3). Segment breaks were oriented approximately perpendicular to the channel flow direction and statistics from each section were calculated separately for the DoD. Missing surveys in Zones A and D resulted in some of the segment comparisons being absent (Table 1).

2.3.5. Channel Morphology Calculations

SfM-derived DoDs are used to calculate changing geomorphology and quantify the sinuosity index (SI) for each section and survey interval. The sinuosity index is a measure of the streamflow path divided by the straight line length of the stream [24] and is important as a key metric used in the Rosgen natural river classification method. The streamflow path was calculated using a lowest-neighbor approach from a starting point selected inside the channel.

SI = (Streamflow path length) / (straight line path).

(3)

SI values are used to categorize a stream channel as straight (<1.1), sinuous (1.1–1.5), or meandering (>1.5) [21]. SI variations over time may be used to quantify the geomorphic dynamism of a rapidly changing fluvial system. The amount of change occurring in the sinuosity index in each section was calculated using a relative percent range (Equation (4)).

ΔSI (relative percent range)_i = ((SI_i − SI_min)/SI_mean) × 100.

(4)

We proposed a novel metric, the morphological activity factor (MAF) (Equation (5)). We used MAF to show the overall changes in the channel when both erosion and deposition are taking place in the same section, often directly adjacent to one another, because of the relationship between the point bar and cutback in fluvial systems. In order to quantify the dynamism of the morphological changes in volume, we computed the absolute volume transfer as a MAF,

MAF = Gain + abs (loss).

(5)

The overall methodology for this study is outlined in the flowchart in Figure A2.

3. Results

The average ground sampling interval of the initial photogrammetric modeling prior to export from Agisoft was 5.22 cm/pix. The average RMSE for SGCPs in the xy direction was 3.44 cm and the average check point RMSE in the xy direction was 12.66 cm. The estimates for the z direction through Agisoft Metashape Pro were not included in the error estimation because the z was manually adjusted to a tolerance of less than 1 cm. The ground sample interval during export from Agisoft was adjusted to 10 cm/pix to reduce the computational work and make a uniform format for DEM merging and DoD calculation.

SfM reconstructions were used to generate accurate DEMs and orthomosaics for four different time steps to use in spatial and temporal analysis of landform change (Table 1). These reconstructions span approximately 3.5 km of the Las Lajas drainage with a total surveyed area of 1.6 km². Volume change approximations can be made by subtracting these high-resolution DEMs at different time steps. The sample DoD in Figure 4, for instance, shows the aggradation and erosion occurring for the entire season between June and October (Figure 4). DoDs are also calculated for time steps spanning 36 to 48 days each, during which as much as 10 m of vertical erosion was evident in some areas. The observed erosion followed normal river meander patterns, though it was remarkably rapid given the energetics of the lahars [24].

The DEM differences allowed high-resolution profile changes to be extracted in key areas throughout the channel reach. For example, the cross-section A–A’ in Figure 5 shows a deep cut in the apex of a meander with 7 m of vertical erosion and 12 m of lateral erosion over the course of the entire rainy season. The steep and unstable walls were often undercut and slough off in this area. In cross-section B–B’, lateral migration and vertical erosion were evident without the extensive volume loss from collapsing channel walls compared to cross-section A--A’. Cross-section B–B’ still represented a strong erosional regime, but with less of an impact due to the relatively minor lateral erosion into the outer bank. Cross section C–C’ showed a 10 m wide constriction formed from an old, indurated deposit. Further downstream, cross-section D–D’ showed greatly reduced erosional features and some aggradation as the stream path freely meandered in the wider unconfined channel.

The volume change in each zone fluctuated throughout the season (Figure 6) and was evident primarily as either losses or gains, whose entire volume tended to cancel out when summed across the ROI (Figure 6a,b,d). In order to quantify the dynamism of the morphological change we computed the absolute volume transfer as the MAF (Equation (5)), which is shown in Figure 6c.

Fuego’s common secondary lahars were triggered by precipitation remobilizing material on the flanks of the volcano. Correlation between rainfall totals from May through October and lahar events was evident, but the amount of rainfall, measured at the FGLJ rain gauge, did not necessarily correspond to the qualitative intensity of a lahar (Figure 7). Since FGLJ was a point location, it might not represent the total rainfall distribution over an area, which fed channels in a variety of watersheds. The relative lahar intensity was a qualitative intensity of the lahar assessed by INSIVUMEH scientists based largely on seismic amplitudes. Values ranged from 1 to 3 on the graph, with 3 being the strongest lahar typically seen in this area. Every lahar event had a corresponding rainfall event, but the converse was not true. Out of the 24 precipitation events with an intensity greater than 2 mm/h, 19 resulted in lahars.

The sinuosity of each 500 m section was calculated separately and suggests that the upper sections (a–e) possessed much more variability in SI than the lower sections (f–i) (Figure 8). The average sinuosity of the entire channel was between 1.4 and 1.5, which, according to Rosgen [24], is at the boundary between sinuous and meandering. Meandering was common in sections a–f, as opposed to the strongly sinuous classification below the indurated deposit knickpoint in section f at km 2.6. The tortuosity of the channel was evident when the stream path was mapped onto the DEM (Figure 9a). The SIs for each section and each period compiled into one figure showed significant variability over time (Figure 9b–d).

SI variation over time is another useful metric to quantify the geomorphic dynamism of a rapidly changing fluvial system. The amount of change occurring in the sinuosity index in each section was calculated using a relative percent range (Equation (4)). The data showed a large shift in sinuosity (ΔSI) from section e to f. Although all ΔSI values were positive, the greatest relative percent range was in section d with 20% and the smallest changes in SI were in sections f and g (Figure 10).

The slope of the channel was measured and compared for each of the 500 m channel sections using the vertical drop from DEMs and horizontal distance from the tortuous flow path. Slopes changed very little over time (less than 1 degree), and this was largely attributable to increasing flow distance due to increasing sinuosity. Slope values for the ROI ranged between −4.6 and 9.0%, as measured in the June survey (Figure 11)

4. Discussion

Quantified erosional and depositional volume rates in the Las Lajas ROI provide insights into processes that occur over both spatial and temporal domains. During the three intra-survey measurements, we note that both the amount of rainfall and number of lahars were relatively constant (Figure 7). Assuming geomorphic change is principally related to the lahar activity, we expected similar volumetric transport when averaged over approximately one-month intervals in the rainy season.

Notably, ΔV_t21 and ΔV_t32 showed opposite behaviors in terms of the net volume change (Figure 6a; sections a–f). We considered the possibility that a single significant event, such as a large wall collapse, might be considerably impactful and could have dammed up the channel following the July survey (t2), resulting in sections b, c, and d showing exceptional deposition. (Figure 6b,d). Section d showed a high MAF with a substantial gain and an even greater loss (for a slight net loss). The DoD illuminated that the loss was from a large portion of the channel wall collapsing and a gain from the subsequent damming and upstream and capture of sediments. Another plausible explanation is inconsistency in sample intervals; data below section e was collected prior to a September 1 small lahar, which preceded the surveys in Zones A and B. That small lahar might, thus, have caused the geomorphic change that was only recorded below Zone B (section e and below) (Figure 6a). It is reasonable that a small lahar with shorter runout may have deposited more material in upper sections of the ROI (Figure 12). The net volume changes calculated for DoD 3 showed a steady increase with increasing flow distance out through section i, suggesting that during this time interval (July through August), some lahars retained enough energy to carry and drop large quantities of material at the lowest extent of the ROI (Figure 6a).

Although the meanders in the stream path appeared relatively dynamic throughout the 2021 rainy season, it is notable that the meandering flow was confined to the gross channel topography that existed prior to 2018. Google Earth time-lapse imagery from 2006 (prior to the construction of the La Reunión golf course), 2018 (prior to the PC), and 2021 each showed features associated with a terraced floodplain and in-filling with potential PC and lahar flow deposits. Las Lajas drainage was relatively stable within the wider channel for some time prior to 1969 according to Google Earth and these early images showed no evidence of recent migration in the form of scoured vegetation and slide scarps (Figure 13A). The floodplain was effectively unimpacted until the catastrophic 3 June 2018 PC event, which overtopped the Las Lajas channel and flowed towards San Miguel Los Lotes [3]. The addition of material into the system filled and dammed the channel temporarily, which diverted the flow (Figure 13C).

Analysis of imagery since 2018 showed that loose PC material within Las Lajas was downcut quickly, resulting in the development of a meandering path over time without any influence from the surrounding or buried topography apart from the outer banks of the old channel. The sinuosity indices of sections a through e were higher and more widely varying over time compared to those from sections f to i (Figure 8 and Figure 10) and we suggest this was because tall unstable walls are easily undercut and collapse more frequently, which diverts the flow within the channel and adds to the material available for bulking.

Observations of how the meanders were confined to the pre-2018 terrace suggest an opportunity to evaluate risk when developing infrastructure and planning land use. The golf course was built partially within the previous channel infill and it is unsurprising that these deposits were impacted by subsequent lahars. In the case of San Miguel Los Lotes, the PC did jump the channel and devastated the community, which was located outside of the Las Lajas drainage proper (Figure 1).

While the upper sections of our study area appeared particularly dynamic, the lower sections were more static. We observed a transition from the upper sections to the lower in section f at a natural constriction, which then became wider before going over a waterfall in section g. This waterfall, just below the bridge, was at an old, indurated deposit with resistant rock serving as a knickpoint. The local base level at the knickpoint accounted for some of the stability of the channel from the constriction above it. The slope of the channel between sections e and f decreased with distances downstream from −7.6% to −5.9%. The channel width, on the other hand, increased from ~10 m at the constriction to between 40 and 60 m before the waterfall. We assumed that much of the energy carrying the type 2 lahar material was lost here. The stability of the geomorphology was reflected in low sinuosity and the cross section (Figure 5D) indicated an evenly mixed regime of deposition and erosion. This was also evident when comparing the high MAF to the low net volume difference. The findings in this study may inform identification of higher risk zones at a course scale in areas with high erosion and lateral channel movement. The placement of the bridge, for example, appeared to be in a low-risk zone for typical type 2 lahars. The location was advantageous, given the combination of the constriction upstream and the indurated layer as a local base level. This was the most stable section of the channel to have major infrastructure.

There were several limitations to this work, motivating future study and improvements. UAV surveys such as the featured examples were limited by weather, battery life, and altitude restrictions. Only a relatively small areal area can be effectively surveyed for analysis and it could be desirable to increase the ROI area for broader volume calculations. The lack of ground control points makes accurate geo-registration impossible. The temporal resolution, at approximate monthly intervals, allows for the generalized impact of many lahar events, but does not allow for study of individual lahar depositional and erosional statistics.

Overall, the Las Lajas channel quickly eroded the PC material at the top of the study area and deposited much of the material at the bottom or beyond. The MAF indicates where major erosion and deposition were simultaneously happening with a net volume change relatively close to zero. The relative percent range of the sinuosity index decreased dramatically where the channel was restricted and changes from erosion-dominated to deposition-dominated behavior. Future work could focus on collecting DoDs with higher temporal resolution in a channel that does not have a constriction. This could be used to investigate if the MAF and sinuosity index can robustly indicate where the change from erosion to deposition occurs.

5. Conclusions

Secondary lahars at Volcán de Fuego are frequent and capable of moving large volumes of material over the course of a rainy season. The constant addition of new material near the vent through volcanic activity, including PCs every few years, combined with heavy precipitation makes Volcán de Fuego a natural laboratory to study the impacts of lahars on the geomorphology of the channels they occupy. Analyses of DoDs in 500 m channel flow sections elucidate the movement of material in both erosional and deposition sense. Using UAV SfM, we quantified volume changes of 490 m³/day loss in the upper sections (~6 km from summit) and a 440 m³/day gain in the lower sections (~10 km from summit).

Erosion and deposition along a lahar path are not solely dependent on the slope of the channel. This was evident using time-lapse volume estimates from four distinct aerial surveys. ΔV_t21, corresponding to volume changes between June and July, showed deposition in several sections that only had erosion in the other periods. The slope in these sections changed negligibly during that span and the number of lahars and amount of rainfall differed very little from the other two periods.

We propose that the MAF is a useful measure of the dynamism of a channel as a function of erosion and aggradation. The net volume change, gross erosion, and gross deposition do not tell a clear enough story of the total activity occurring in the channel. The MAF allows for a quantitative measure of the total activity occurring within the channel. The area where the regime switches from erosional to depositional is the area where the MAF has the highest value due to the shifting channel eroding and depositing at near-equal rates.

The shape of the channel and some underlying indurated deposits significantly influenced the behavior of the channel in sections e through i. In section e, there was a constriction 12 m wide made of an indurated deposit, directly after a bend in the channel, which slowed the flow and presumably shed some of the energy. The constriction opened into a 60 m wide channel that flowed over a knickpoint formed from another indurated deposit, creating a local base level and stabilizing the stream behavior concerning sinuosity index and downcutting. Below the knickpoint, lahars ran out into a wide channel and gradually deposited material as they lost momentum.

The sinuosity index showed a distinct change in character from the upper sections, which were tortuous up high and less so further down Las Lajas, starting at section f. An increase in sinuosity over time was observed in the upper sections and can be a useful geomorphological statistic because it shows how lahars are actively changing flow paths.