Twenty-Five Years of Geomorphological Evolution in the Gokurakudani Gully (Unzen Volcano): Topography, Subsurface Geophysics and Sediment Analysis

: In the aftermath of pyroclastic density current-dominated eruptions, lahars are the main geomorphic agent, but at the decadal scale, different sets of processes take place in the volcanic sediment cascade. At Unzen volcano, in the Gokurakudani gully, we investigated the geomorphologic evolution and how the topographic change and the sediment change over time is controlling this transition. For this purpose, a combination of LiDAR data, aerial photography and photogrammetry, ground penetrating radar and sediment grain size analysis was done. The results show choking zones and zones of enlargement of the gully, partly controlled by pre-eruption topography, but also by the overlapping patterns of the pyroclastic ﬂow deposits of 1990–1995. The ground penetrating radar revealed that on top of the typical lahar structure at the bottom of the gully, side wall collapses were trapping ﬁner sandy sediments formed in a relatively low-energy deposition environment. This shows that secondary processes are taking place in the sediment transport process, on top of lahar activity, but also that these temporary dams may be a source of sudden sediment and water release, leading to lahars. Finally, the sediments from the gully walls are being preferentially oozed out of the pyroclastic ﬂow deposit, meaning that over longer period of time, there may be a lack of ﬁnes, increasing permeability and reducing internal pore pressure needed for lahar triggering. It also poses the important question of how much of a past event one can understand from outcrops in coarse heterometric material, as the deposit structure can remain, even after losing part of its ﬁne material.


Introduction and State of the Art
Quantifying landforms by balancing erosion and uplift on mountain and volcanic slopes is a century-old quest: e.g., [1]. Although long-term evolution models have evolved to include the complexity of soil generation and slope angle-related transport [2], scientists have argued that the models would benefit from further direct measurements (e.g., at Merapi Volcano [3]), which unfortunately are often too scarce on volcanoes. Furthermore, these measurements also need to be extended over long periods of time to account for the variability in the sediment cascade at the decadal scale (e.g., at Unzen Volcano [4]) and at the Holocene-and Quaternary-scale period. Although there is a large body of research on syn-and post-eruptive research as well as long-term (Holocene, Quaternary and older Geosciences 2021, 11, 457 2 of 24 period) research, work at the decadal scale is still in need of attention, including at Unzen Volcano, the location of the present research [4].
At Unzen Volcano (Japan), despite a large amount of unconsolidated material remaining on the apron of the volcano, the frequency of lahars and other forms of gully erosion has decreased dramatically in recent years [4], with a change in the relation between eroded material and rainfall intensity and duration, as well as a change in the watershed geometric configurations [4]. Tsunetaka et al. [4] investigated the role of the geomorphology and the rainfall on the erosion, and in the present contribution, the goal is to qualify and quantify the role of material composition change over time and the role of the sedimentary structure on the geomorphologic evolution of the Gokurakudani Gully, and so at the decadal scale after the eruption.

Syn-Eruptive to Decadal-Scale Erosion on Volcanoes
On stratovolcanoes, sediment fluxes increase sharply in the aftermath of eruptions. At Merapi Volcano (Indonesia), out of all the 851 lahars produced during the 44 years of existing records, 11.2%, 9.2% and 15.2% occurred immediately after respectively the eruptions of 1969, 1975 and 2011 [5]. For the millennial eruption of 2010, 0.03 to 0.06 km 3 pyroclastic flow deposits were generated [6], resulting in more than 250 lahars in the Kali-Putih valley alone, in less than half a year [5]. In the month after the eruption stopped, more than 50 lahars were counted in the Kali-Putih valley, around 25 in the Kali Boyong valley and 25 in the Kali-Gendol valley [7]. At Mount Pinatubo, an approximated 6 km 3 of loose volcanic sediments mantled the volcanic slopes due to the 1991 eruption [8], and 450 lahars were reported for the Pasig and Sacobia catchments-two of the eight catchments impacted by the eruption-between 1991 and 1997 [9]. At Mount St. Helens, an approximated 2.3 km 3 of pyroclastic material was deposited in the upper Toutle River in 1980, feeding syn-eruptive lahars that transported about 140 × 10 6 Mg (megagrams) of sediment in 1980 alone [10].
The syn-eruptive and immediately post-eruptive sediment flux decreases over time. In the aftermath of the 1984 eruption at Merapi Volcano, the number of lahars per year reduced from 35 in 1985-1986, to 27 in 1986-1987, to 29 in 1987-1988 and to 0 in [1988][1989][1990] in the Kali Putih [11]. At Mt. St. Helens, the recorded sediment fluxes in different valleys ranged from 1 to 100 Mg/km 2 /year (mass per catchment surface per unit time) in 1980, and they reduced to only 0.2 to 0.3 Mg/km 2 /year in 1994 [12].
In controlling the flux, both rainfall and the sediment availability play a crucial role. At Merapi Volcano, the largest number of monthly lahars was triggered just after the 2010 eruption, with 50 to 60 lahars for a month when the cumulated rainfall events were 1500 mm to 2000 mm, while the highest rainfall in the subsequent months (>3000 mm) triggered fewer lahars: 40-50 per month [7]. The timing of the rainfall and the availability of removable material is therefore crucial. This pattern could also be observed at Mt. Pinatubo, where broadband seismometers show an increase in high sediment concentration events between 1991 to 1994, corresponding to increased-intensity rainfalls over the years. However, from 1994 to the end of 1997 the number of debris flows and hyperconcentrated flows decrease every year, leaving space for lower-density flows, regardless of the rainfall intensity. At Mt. St. Helens, because the period 1980-1994 was relatively dry, the decreasing trend in sediment fluxes in the first 14 years after the volcanic collapse suddenly ramped up in 1995-1996, starting a new decreasing trend [11].
From these events, it appears that at the decade and multi-decade scale, the sediment fluxes mostly carried by lahars are first controlled by the amount of material available, before being then controlled more predominantly by rainfall trends. Despite a relatively rapid (considering the life of a volcanic structure) decrease in lahar frequency, their deposits dominate the material volcanic aprons are made of [8,13,14].

Quaternary-Scale Erosion Rates on Volcanoes
At the Quaternary scale it becomes more difficult to identify the influence of one process over another, but for stratovolcanoes, mechanical erosion is the main driver compared to shield and island volcanoes where chemical weathering dominates the erosion process-e.g., 49 to 306 t/km 2 /year from river transport and 113 to 1382 t/km 2 /year as dissolved material in groundwater, at Mayotte French Volcanic Island [15]. At the stratovolcanoes Cayambe, Cubilche, Cushnirumi, Cusin, FuyaFuya, Imbabura and Mojanda, on the Ecuadorian Arc, the estimated erosion rates from isotopic analysis are respectively 0.02, 0.02, 0.12, 0.02, 0.1 and 0.14 km 3 /1000 years, using calculations over the periods of 30,000 years to 1,050,000 years [16]. These rates are in line with erosion rates provided in other studies, such as 0.2 km 3 /year at the Tungurahua volcano and the Huisla volcano [16], and at the Cayambe volcano [17]. On top of the geochemical approach, scientists have used high-resolution topography to measure erosion, following the movement of knickpoints at waterfalls for instance. In this way, Hayakawa et al. [18] measured recession rates of 0.013 to 0.068 m/year at 11 waterfalls located in the Aso volcanic structure (Japan), and 0.08 to 0.15 m/year at Shomyo falls, cut into ignimbrites [19]. These calculations and measures are to be taken carefully, because they encompass climate and landcover change, etc. For instance, at Terceira Island in the Azores, the change from pasture to bare soil has shown a threefold increase of soils erosion based on measurements from 1996 to 1998 [20], showing that minute changes can have important impact on long-term erosion rates of volcanoes. Furthermore, the breath of processes occurring at different timescales after an eruption, makes the interpretation of these calculations even more arduous.
Between the first syn-and immediately post-eruptive scale and Quaternary scales, there is thus a relative knowledge gap at the tens to hundreds of years scale.

Lahar: The Major Observed-Agent of Change
From field surveys and event observations, lahars appear to be an important terragenic driver of stratovolcanoes. Since its first use in the scientific literature-almost 100 years ago for a mudflow in Alaska [21]-lahar research has flourished, especially due to its imperatives for associated hazards and disaster risk. In the scientific literature, a lahar has been defined as an Indonesian term most commonly defined as a rapidly flowing, gravity-driven mixture of rock, debris and water from a volcano [22][23][24].
Large events will typically have several flow peaks and phase changes [25], and the lahar will also vary with time and distance downstream (dilution, tributary intakes and material deposition). The lahar may comprise one or more flow types, which include debris flow, transitional or hyperconcentrated flow and muddy streamflow or flood flow [26], all differentiated based on the sediment concentration [27]. Typically, the debris flow phase has a sediment concentration (SC) of 60% of volume or 80% of mass [28], and the hyperconcentrated flow phase has a SC between 20 and 60% of volume or 40-80% of the mass [24,29]. As the lahar often dilutes in hyperconcentrated flow and eventually streamflow [30], the deposited material at a given site only represents a portion of the originally mobilized material, which makes it often haphazard to survey post event, notably using traditional outcrop analysis [31].
Because of the wide range of environments where lahar can occur [32], and because lahars can happen during the timespan of an eruption and up to several decades after it has stopped (e.g., hot and cold lahars [33]), the processes affecting lahar triggering are assumed to vary with the environment and with the local geological and geomorphologic setting. Notably, the proportion of fines delimitates cohesive and non-cohesive lahars. Non-cohesive lahars are most common when Dacitic (e.g., Unzen Volcano in Japan, with poor grain size distribution in clay size material [34]) and Andesitic material (e.g., Merapi and Semeru in Indonesia) are present. Such lahars usually originate from dome collapse pyroclastic flows and other forms of pyroclastic flow material, where the number of fines is relatively limited. On the contrary, the development of cohesive lahars requires a larger amount of clay, which can be produced via hydrothermal eruptions [35] or from Geosciences 2021, 11, 457 4 of 24 flank-collapse debris avalanche deposits for instance [27]. Further variability in lahar and the related deposits appear over time, especially when the timespan from an eruption increases [36]. As the progressive sediment depletion after an eruption modifies the relation between rainfall and lahar triggering, and as the geometry of valleys changes [37], so are the erosion rates and processes (e.g., at Merapi and Semeru Volcano [14]). Long-term monitoring is therefore essential to grasp the variability in lahars over time and their interplay with the valley geomorphology [37].

From Pyroclastic Density Currents' Deposits to Lahar Deposits
Consequently, sediment flux and lahar research are often the link between pyroclastic flow deposits and their transformation into lahar deposits, with the finer portion of the material being exported further downstream, generating other forms of deposits. Furthermore, the irregular nature of those events often requires scientists to work from the deposits, where pyroclastic flow deposits mingle with lahar and other deposits generated by the remobilization of the original material.
Pyroclastic density current deposits (PDCD) typically display tongue-shaped topography with a front lobe, and they are often trapped in the valleys, while surge deposits are a thinner deposit spreading around it. The internal structure of PDCDs is typically made of sets of bedded layers with or without grain-size gradients within the layers [38][39][40]. From their remobilization, lahar deposits are typically clasts rich and supported by a sandy to gravelly matrix with sets of crude horizontal layers [41], which alternate with foresets and backsets depending on the valley configuration [42]. Local variations exist along the valley gradient and at the sub-meter scale [31]. Those intricacies have been related to multiple peak flows and eventually complex deposition patterns during single events [25]. Over several years to decadal scales, the topography of lahar deposits and the gullies they shape experience a complex evolution with an alternation of erosion and deposition, notably under the influence of ravine wall collapses [37]. The work of Vasquez et al. [37] in the Montegrande ravine of Volcan de Colima shows that after a single event, erosion can occur in portions of the channel, while wall collapses and floor filling can temporarily raise the valley floor by more than 2 m in other portions of the channel. In detail, it is therefore not erosion or depositional processes, it is a mix of the two, with dominant trends at a broader scale. As the floor slope of the ravine tends toward a general equilibrium slope, upstream sections can still be dominated by erosion (8 to 10 degree slopes) while the lower slopes (2 to 4 degrees) are more prone to deposition. At Unzen Volcano, it is this type of PFD-dominated environment from which lahars are emerging.

Topographic and Geological Setting
The Unzen volcano is located on the peninsula of Shimabara, on the Island of Kyushu (South Japan), and the actual summit is at 1486 m asl (Figures 1 and 2). The basement rock of the volcano is Paleogene to middle Pleistocene sedimentary sequences. Japanese scientists have recognized two phases in the edifice construction: the Older Unzen volcano between 500 kBP and 100 kBP and, after 100 kBP, the Younger Unzen volcano. During the last period, four volcanic structures have topographically emerged: the Mayuyama volcano (Figure 3), the Nodake volcano (Figure 3), the Myokendake volcano and the Fugendake volcano ( Figure 3).        [43]. The Volcanic apron and fan that starts at the present summit of Mt. Unzen and extends towards the east is trapped in the volcanic graben, in which the newest structure is sinking.

The Unzen Volcano and the 1990-1995 Eruption
After a period of 198 years of dormancy, its eruption began in November 1990 and became infamous after taking the lives of 44 victims including French and American scientists and Japanese mass media professionals in June 1991. The eruption lasted until 1995. The eruptive sequence started with a swarm of earthquakes in November 1989. The first tremor underneath the summit was then recorded in July 1990. The eruption officially ended in February 1995 with the halt of the effusion and a slow caving of the dome [44]. The eruption resulted in a series of dome growth, fracturation and subsequent gravity collapse, sweeping mostly over the slopes located on the east side towards the Mizunashigawa valley basin. However, during the latest stage of the eruption, pyroclastic density currents also flowed northbound. On the southern flank, the pyroclastic flow deposits filled the valley with material >50 m thick [45], greatly modifying the pre-eruption topography and by then the subsurface hydrological mechanisms and processes. From the PDCDs (probably mixed with other secondary deposits), lahars have been generated, sweeping the flanks of the volcano, with the majority occurring in the the Mizunashigawa valley at first and in recent years in the smaller-scale gullies of the Gokurakudani and the Tansandani [4], where the lahars have been scouring the 1990-1995 eruption material down to the pre-eruption material [46][47][48].

Methodology
To reach the objective, the research combines (1) fieldwork, (2) laboratory analysis of the samples and collected empirical data and, (3) geophysical exploration based on ground penetrating radar (GPR).  [43]. The Volcanic apron and fan that starts at the present summit of Mt. Unzen and extends towards the east is trapped in the volcanic graben, in which the newest structure is sinking.

The Unzen Volcano and the 1990-1995 Eruption
After a period of 198 years of dormancy, its eruption began in November 1990 and became infamous after taking the lives of 44 victims including French and American scientists and Japanese mass media professionals in June 1991. The eruption lasted until 1995. The eruptive sequence started with a swarm of earthquakes in November 1989. The first tremor underneath the summit was then recorded in July 1990. The eruption officially ended in February 1995 with the halt of the effusion and a slow caving of the dome [44]. The eruption resulted in a series of dome growth, fracturation and subsequent gravity collapse, sweeping mostly over the slopes located on the east side towards the Mizunashigawa valley basin. However, during the latest stage of the eruption, pyroclastic density currents also flowed northbound. On the southern flank, the pyroclastic flow deposits filled the valley with material >50 m thick [45], greatly modifying the pre-eruption topography and by then the subsurface hydrological mechanisms and processes. From the PDCDs (probably mixed with other secondary deposits), lahars have been generated, sweeping the flanks of the volcano, with the majority occurring in the the Mizunashigawa valley at first and in recent years in the smaller-scale gullies of the Gokurakudani and the Tansandani [4], where the lahars have been scouring the 1990-1995 eruption material down to the pre-eruption material [46][47][48].

Methodology
To reach the objective, the research combines (1) fieldwork, (2) laboratory analysis of the samples and collected empirical data and, (3) geophysical exploration based on ground penetrating radar (GPR).

Data Acquisition
The topographic data used in the present research are LiDAR (light detection and ranging) data for the years 2004 to 2016, completed with UAV-based data using a DJI Mavic Pro UAV (DJI, Shenzhen, China) for the years 2017 and 2018. The UAV-based DEM was only used to rectify aerial imagery in the present work. The LiDAR data was commissioned by the Unzen Restoration Work Office directed by the Ministry of Land, Infrastructure, Transport and Tourism, while the UAV data were collected by the authors. Details on the dataset for both the UAV and the LiDAR can be found in Tsunetaka et al. [4].

UAV Data Processing
Structure from motion (SfM) is the photogrammetric method used with the UAV and the implementation was performed using Agisoft PhotoscanPro ® software, from which the sparse point cloud, dense point cloud, the meshed surface, the DEM and the orthophotographs were successively generated following methods described in the literature [49][50][51] and used previously in various volcanic and non-volcanic settings by the authors [45]. To calibrate and scale the DEMs, eight targets in the valleys were measured using GNSS PPK with a Trimble GeoExplorer (Geo-7X, Trimble, Sunnyvale, CA, USA) hooked to an external antenna. To ensure the best data generation, tie points were added through a trial-and-error approach. For further details on the data please see Tsunetaka et al. [4].

GIS Analysis
Planform analysis-From aerial orthophotographs of the years 2005, 2006, 2015 and 2016, a total of 50 gully sections of 30 m length were digitized and calculated in UTM 52, using the open-source GIS environment QGIS. Each section was digitized using polygonshapefiles and the limits of each cross section were kept constant for the 4 different years, so that for a specific section, widening or size reduction could be retrieved. The planform analysis is thus providing a dataset of the planform stability/instability of the gully.
Topographic analysis-Imagery and topographic data was incorporated in the QGIS environment in UTM 52. The raster dataset was sampled using vector files (shapefiles). A set of 116 cross sections was acquired along the Gokurakudani gully with a sampling interval of 10 m.
On each cross section, 100 points were acquired from the LiDAR data of 2005, 2015 and 2016, creating a dataset of 348,000 points to compare at 116,000 locations along the gully (data available upon reasonable request). The data of each cross section was used to calculate the cross-sectional area, in order to detect vertical and horizontal adjustments in the geometry of the gully using the following formula for each cross section: where A is the cross section area, n is the number of sampling points along the cross section, i is the iterator over n, d i is the distance from the first point of the cross section to the ith point, and the subscript i + 1 is the iteration to the next point; Z max is the maximum altitude at one of the n points on a given transect, and Z i is the averaged altitude given by: The calculation was made in the R language, using the shapefile id (each transect has a distinct id number) to run the above calculations on each transect separately.

Characterization of the Gully Floor and the Walls' Material
After characterizing the change and rates of change in the gullies, it is then important to understand the role of the material and its structure (next section). In Tsunetaka et al. [4], the questions of rainfall and modifications of the geometry of the catchment have been examined, and to complement this research, the authors look at how the material change may also influence the evolution of the gullies. Indeed, as erosion occurs, the different size fractions travel at differential speed, creating a natural horizontal sieving between the walls (dominated by PFDs), the floor of the gully and its exit. A total of 18 samples were collected from the floor of the valley with a shovel or from the walls. For the walls, the outcrop was refreshed to remove potential contamination from falling material, except in the case of the acquisition of the wall drape (see the results for explanation), composed of the fine fraction that had oozed out of the wall. The sediments were sealed and double bagged in 800 mL plastic bags.

Dry Sieving and Grain Size Analysis
At the laboratory, each sample was weighed, and then the organic debris were floated out of the sediments using water washing and drying. The amount of organic debris found was close to none in most samples, if not for a few leaves and branch bark trapped in the gully floor samples. Their amount was judged negligible. After drying the remaining sediments in a dish on a hot plate at 60 degrees, the samples were dry-sieved, and the size distribution was measured per mass for a column of 8 mesh sizes (8 mm, 4 mm, 2 mm, 1 mm, 0.5 mm, 0.25 mm, 0.125 mm, 0.0625 mm, <0.0625 mm). The grain size distribution was calculated and so was the grain sorting using the grain size variation using [52] inclusive standard deviation: Finally, the kurtosis, or measure of the strength of the peak was derived as:

Subsurface Geophysical Data Using Ground Penetrating Radar (GPR)
To complete the material characterization, the gully floor subsurface was imaged with ground penetrating radar (GPR). GPR was used in order to investigate the internal structure and find the presence (or not) of bedrock that may not be visible from the surface following the methods used on PDCDs [39,40] and lahar deposits [31,42,44]. The method was also aimed at better understanding the evolution of the gullies: if the channel is going through a purely erosive phase, one is expected to find solely PDCDs or the pre-eruptive material; however, if a mix of lahar deposits, wall collapse material, etc. is present, it shows a non-linear evolution. This is essential, because even yearly topographic surveys are blind to topographic variations that occur at the minute to hour scale.

The GPR Method
Subsurface data was gathered using GPR from the floor of the gully. GPR is a geophysical technique based on the emission and reception of electromagnetic waves through the ground. As the wave travels in the substratum and encounters irregularities, part of the energy is refracted and recorded by the receiving antenna [53]. Electromagnetic waves travel at the speed of light in a vacuum and decrease in substratum, mostly as a function of the electricity resistivity (ξ r ) and the relative magnetic permeability (µ r ), so that for a wave travelling back and forth from the emitting antenna (E t ) and a punctual object or a layer (O), the time (t) can be calculated using the slope of diffraction hyperbola as: where x is the horizontal distance and h is the vertical distance underneath the antenna, and V 1 is the velocity until the object that creates the hyperbola; and because the velocity (V) can also be estimated to be: where c is the speed of light (considered to be a constant in a vacuum), information on the material invisible at the surface can be inferred. This calculation is made using the relation of electric resistivity and velocity between layers using the reflection coefficient R c : where the indices '1' and '2' indicate two different layers of contrasting dielectric properties.

Field Data Acquisition with the Mala-ProEx GPR
The GPR data acquisition was performed at two locations, one along a set of transects of~140 m at the downstream end of the Gokurakudani gully using a 500 MhZ and a 800 MhZ antenna (centered on UTM 52 623163/3624407) and upstream in a reach of the Gokurakudani gully, where material had accumulated UTM 52 623353/3624242, and from which 4 GPR transects were collected: one along the gully direction and three perpendicularly across (cf. Figure 1 with the locations).
For the present study, the GPR is a Ramac ® Pro-Ex (Mala GeoScience, Mala, Sweden) mounted with shielded antennas 500 MhZ and 800 MhZ. The distance was measured using a calibrated coding wheel and a measuring tape to confirm the wheel's data, as well as a hand-held GNSS to mark the track position in the Gokurakudani gully. The data used in the present research follows semi-horizontal topographies for which no complex topography inversion was necessary (cf. Figure 1 for location).

GPR Data Processing
The data produced by the Mala system are .RDS and .RD3 binary files that were converted into ascii data using ReflexW (Sandmeier Software, Karlsruhe, Germany). The processing steps used in the present contribution are standard for GPR [41]: -Static time correction to remove the signal time before entry into the ground; -DEWOW filtering (low-cut filter) to remove the air-ground reflection (this step was removed from one of the transects as it was erasing some of the data associated with the removed signal); -AGC gain amplification with a 40 ns window and a 1.12 multiplier; -Removal of the background noise with an average moving window every 100 traces; -Energy propagation loss compensation using an energy decay multiplier of 1.44; -Measure of the velocity through the ground. As the common mid-point method is not possible with the shielded antenna [42], the authors used the slopes of the hyperbola curves to calculate the velocity profiles; -Finally, the velocities are then converted into depths using the calculated velocities.
Because the data were acquired on semi-horizontal surfaces, no topographic data migration was necessary.
The method therefore combines imagery with topographic analysis with sediment grain size characterization and GPR-based sediment architecture measurements.

Planform Changes of the Gokurakudani Gully
Derived from aerial and UAV photographs for the period 2005-2016, the planform evolution of the Gokurakudani gully shows surfaces between 0 m 2 to 1500 m 2 for the surveyed series of 30 m length sections of the gully (Figure 4).

-
Finally, the velocities are then converted into depths using the calculated velocities.
Because the data were acquired on semi-horizontal surfaces, no topographic data migration was necessary.
The method therefore combines imagery with topographic analysis with sediment grain size characterization and GPR-based sediment architecture measurements.

Planform Changes of the Gokurakudani Gully
Derived from aerial and UAV photographs for the period 2005-2016, the planform evolution of the Gokurakudani gully shows surfaces between 0 m 2 to 1500 m 2 for the surveyed series of 30 m length sections of the gully (Figure 4). During the 11 years, the gully is irregular in its width with broad and narrow sections ( Figure 5). Although one could expect a homogenization of the gully over time under the effect of lahars and other processes, the gully has however kept the same spatial distribution of wide sections and narrower sections. The distribution of the gully widths shows a narrow section in the first upstream 800 m, then an abrupt widening for almost 500 m, and then a narrowing again, with both the downstream and upstream narrow sections experiencing widening. The two upstream sections widen from a <190 m 2 /30 m (6.3 m width) to a 750-910 m 2 /30 m length (25 to 30 width), representing a rapid average lateral growth of 2.5 m/year for the section. The section between ~400 and ~700 m from the upstream ends also widens dramatically in the 11-year period, passing from <190 m 2 /30 m to 910-1500 m 2 /30 m, or an average growth of 3 m/year ( Figure 5). By comparison, the central section remains mostly unchanged. The planform evolution of the gully is therefore irreg- During the 11 years, the gully is irregular in its width with broad and narrow sections ( Figure 5). Although one could expect a homogenization of the gully over time under the effect of lahars and other processes, the gully has however kept the same spatial distribution of wide sections and narrower sections. The distribution of the gully widths shows a narrow section in the first upstream 800 m, then an abrupt widening for almost 500 m, and then a narrowing again, with both the downstream and upstream narrow sections experiencing widening. The two upstream sections widen from a <190 m 2 /30 m (6.3 m width) to a 750-910 m 2 /30 m length (25 to 30 width), representing a rapid average lateral growth of 2.5 m/year for the section. The section between~400 and~700 m from the upstream ends also widens dramatically in the 11-year period, passing from <190 m 2 /30 m to 910-1500 m 2 /30 m, or an average growth of 3 m/year ( Figure 5). By comparison, the central section remains mostly unchanged. The planform evolution of the gully is therefore irregular in space, and the gully is not homogenized nor does it become calibrated in its horizontal shapes. Temporally, the gullies have not evolved in a homogeneous manner either. The originally disconnected section upstream of the gully fails to grow after connecting, and between 2008 and 2016, there is a stabilization and a narrowing of the gully, with more vegetation encroachment. and between 2008 and 2016, there is a stabilization and a narrowing of the gully, with more vegetation encroachment.
In detail, the evolution of these sections is not a linear process and it is notably driven by the erosion of the side walls and the creation of erosion notches, corresponding to water preferential pathways that have eroded the sediments (Figures 4 and 5). Comparing more recent data obtained from UAV imagery in 2017 and 2018, one can observe the erosion of side gullies connected to the main Gokurakudani gully ( Figure 5).

Vertical Changes of the Gokurakudani Gully
The comparison of topographic data between 2008 and 2016 shows that the gully is stable overall, and its floor has not been eroding nor accumulating deposits significantly, as the majority of the topographical change is within 50 cm ( Figure 6). There are three unstable zones between 250 and 300 m and between 550 and 600 m and at 1100 m. These features appear as a topographical step, with the top of the step being stable, while the In detail, the evolution of these sections is not a linear process and it is notably driven by the erosion of the side walls and the creation of erosion notches, corresponding to water preferential pathways that have eroded the sediments (Figures 4 and 5). Comparing more recent data obtained from UAV imagery in 2017 and 2018, one can observe the erosion of side gullies connected to the main Gokurakudani gully ( Figure 5).

Vertical Changes of the Gokurakudani Gully
The comparison of topographic data between 2008 and 2016 shows that the gully is stable overall, and its floor has not been eroding nor accumulating deposits significantly, as the majority of the topographical change is within 50 cm ( Figure 6). There are three unstable zones between 250 and 300 m and between 550 and 600 m and at 1100 m. These features appear as a topographical step, with the top of the step being stable, while the lower part erodes. The recorded change has been 1.5 m to 2 m per year. There are three zones of such erosion, creating drops of several meters in the gully floor ( Figure 6) for the studied period. On top of the three steps with erosion >1.5 m, there are smaller secondary steps that also exist with erosion values between 0.5 m and 1 m.  The gully is therefore changing over time in an irregular manner, both horizontally and vertically, with the alternation of zones of relative stability against zones of increased erosion. Completing the longitudinal profile evolution, cross-profile analysis using the angles of the wall face as a measure of stability ( Figure 7) show that the sub-vertical walls are all comprised between 26.7 and 71 degrees. They present a sloping angle comprised between values below 35.4 degrees (the angle of repose or friction angle) and up to subvertical angles (71 degrees). The distribution of these angles is controlled by the status of the wall: whether it has collapsed or not and also its position on the right or left side of the gully (Figure 7). Dividing the walls that have shown signs of collapse from those that have not, the data show that the true left walls can sustain an angle of 56.7 degrees, while the true right walls can withstand angles of up to 68 degrees. The collapsing mechanisms are arguably different on the true right and the true left of the gully. This is further evidenced by the distribution of collapsed or partially collapsed wall profiles that can be straight up to 71 degrees on the true right, while the value only reaches 45 degrees for the true left (Figure 7). The "not collapsed" walls show the same pattern with maxima at 68 and 56.7 degrees for the true left and the true right respectively. The gully is therefore changing over time in an irregular manner, both horizontally and vertically, with the alternation of zones of relative stability against zones of increased erosion. Completing the longitudinal profile evolution, cross-profile analysis using the angles of the wall face as a measure of stability ( Figure 7) show that the sub-vertical walls are all comprised between 26.7 and 71 degrees. They present a sloping angle comprised between values below 35.4 degrees (the angle of repose or friction angle) and up to subvertical angles (71 degrees). The distribution of these angles is controlled by the status of the wall: whether it has collapsed or not and also its position on the right or left side of the gully (Figure 7). Dividing the walls that have shown signs of collapse from those that have not, the data show that the true left walls can sustain an angle of 56.7 degrees, while the true right walls can withstand angles of up to 68 degrees. The collapsing mechanisms are arguably different on the true right and the true left of the gully. This is further evidenced by the distribution of collapsed or partially collapsed wall profiles that can be straight up to 71 degrees on the true right, while the value only reaches 45 degrees for the true left (Figure 7). The "not collapsed" walls show the same pattern with maxima at 68 and 56.7 degrees for the true left and the true right respectively.

Cross Section Areas' Evolution from 2005 to 2016
Using the sampled transects along the Gokurakudani gully (Figure 4), the statistical distribution of cross sectional areas of the gully was calculated. It shows that the distribution of the gully areas is bimodal in 2005, becoming increasingly unimodal in 2015 and 2016 ( Figure 8).

Cross Section Areas' Evolution from 2005 to 2016
Using the sampled transects along the Gokurakudani gully (Figure 4), the statistical distribution of cross sectional areas of the gully was calculated. It shows that the distribution of the gully areas is bimodal in 2005, becoming increasingly unimodal in 2015 and 2016 ( Figure 8).

Cross Section Areas' Evolution from 2005 to 2016
Using the sampled transects along the Gokurakudani gully (Figure 4), the statistical distribution of cross sectional areas of the gully was calculated. It shows that the distribution of the gully areas is bimodal in 2005, becoming increasingly unimodal in 2015 and 2016 ( Figure 8).  In 2005, the two modes of the cross sections are~200 m 2 and~400 m 2 with respective density of 0.003 and 0.004. In 2015 and 2016, the series tends to converge with one single mode~380 m 2 , but this mode is a plateau rather than a peak and the plateau ranges from 300 m 2 to 500 m 2 . There is thus an increase of the sizes of the cross section areas and a thickening of the positive tail, while the density of small cross sections <200 m 2 reduces (Figure 8). Overall, one can observe a shift of all the series towards larger values. The cumulated density representation of this dataset shows this pattern with a shift to high values. In 2005, a cross section <200 m 2 represents >20% of the cross sections, when in 2015 and 2016, it drops to~10% (or 0.1 on Figure 8).

Topographic Results' Interpretation
The geometry of the Gokurakudani gully displays stretches where the width alternates between narrow and wide reaches. It also shows a relative stability of the bed, with localized points of erosion, and the lateral erosion and stability of the walls of the gully show a dissymmetry between the left and the right wall. To attempt to understand those differences, the interpretation of these results will rely on the hypothesis that antecedent conditions must have a control on the geometry and geomorphological evolution of the gully. The upper section of the Gokurakudani gully cuts into the 1990-1995 PDCDs before wrapping around the pre-eruptive topography (PeT on Figure 9), so that the left is against it, while the right is limited by the PDCD. For the lower half of the gully, it cuts through the 1990-1995 PDC, which is present on both sides of the gully. Cross sections show that the PDC topography is slanting away from the gully right wall, while at transect C-D ( Figure 9) the transversal slopes converge towards the gully. At transects E-F and G-H, the transversal slope diverges from the gully right wall, and converges towards the left wall. Finally, at transect I-J, the transverse slope tends towards 0 (flat topography). The Tansandani gully is also much lower than the Gokurakudani gully. For most of the gully, the direction of the transverse topography is opposite for the left and right walls of the gully. From the topographic surface of the PDCDs (Figure 9), this characteristic is related to the location of the gully compared to the PDCD and the pre-eruptive topography. This dissymmetry between the left and the right side of the gully can be related to the difference between the stability of the left and the right wall, as precipitated water tends to either diverge or converge towards the gully. From the aerial imagery ( Figure 5), the direction of the topography and the location of the side gullies are related. In places where the transversal slope converges towards the gully, side gullies are present, otherwise not. Most side gullies are therefore located on the left side of the Gokurakudani, except in the mid-upper section, where both transversal slopes converge towards the Gokurakudani, and then side gullies are denser on the right side. One will also note the lack of gullies in the PeT material, as it is not loose material generated by the latest eruption.

The Material: Grain Size Analysis and Material Differentiation
Although the geomorphology of the gully and its evolution over time and space show variations that can be interpreted as topography-related, the constituents of the gully must have changed from the original PDCDs and the secondary processes (such as lahars and waterflows as evidenced by [4]. The grain size distribution of the collected samples from the valley floor ( Figure 10) shows a material dominated by poorly sorted material (>1 on the sorting scale) with a platykurtic distribution (i.e., the tails at both ends of the distribution are not important). For instance, there is no material of fine sand size and smaller in sample D02 taken at the center of the gully floor, and the gravel fraction is also limited. Sample D01, (interpreted as a lahar debris flow phase facies, based on the material and the rest of the lahar deposit in the Gokurakudani gully) has even less fine material but a larger fraction of coarse material and gravel. The other samples from the gully floor, D03 to D06, were deposited by Newtonian flows and trapped behind a side gully debris cone.

The Material: Grain Size Analysis and Material Differentiation
Although the geomorphology of the gully and its evolution over time and space show variations that can be interpreted as topography-related, the constituents of the gully must have changed from the original PDCDs and the secondary processes (such as lahars and waterflows as evidenced by [4]. The grain size distribution of the collected samples from the valley floor ( Figure 10) shows a material dominated by poorly sorted material (>1 on the sorting scale) with a platykurtic distribution (i.e., the tails at both ends of the distribution are not important). For instance, there is no material of fine sand size and smaller in sample D02 taken at the center of the gully floor, and the gravel fraction is also limited. Sample D01, (interpreted as a lahar debris flow phase facies, based on the material and the rest of the lahar deposit in the Gokurakudani gully) has even less fine material but a larger fraction of coarse material and gravel. The other samples from the gully floor, D03 to D06, were deposited by Newtonian flows and trapped behind a side gully debris cone. Figure 10. Gokurakudani gully floor samples (UTM 52S 623353/3624242). Sample D01 was retrieved from a pseudo-terrace on the true right, made of a mixture of material without fabric, which can be attributed to a debris flow deposit. Sample D02 was retrieved from the valley floor, near the gully wall. Samples D03 and D07 were retrieved from small side gully debris cones with a gentle slope of 7.02 degrees. Samples D04 to D06 were all retrieved in the central part of the gully where the GPR transects 7, 8 and 9 were taken. Erosion features and the GPR transects show that the material has been deposited by Newtonian flow with a sub-horizontal layering.
In the lower section of the gully, lahar deposits S09 to S12 are poorly to moderately sorted, and the distributions are very platykurtic (Figure 11). This material contrasts with the wall PDCDs that are at the limit between platykurtic and mesokurtic. If the sampled PDCD material did not display a distinctive difference with the lahar deposits, a drape of fine and wet material was oozing out of the outcrop and its sample showed a significantly finer distribution of material (S08 on Figure 11). This shows that although the wall is stable and not changing from a topographic perspective, its constituents are being selectively washed away. The diminution of fine material from this process may explain why the PDCD material in sample S13 does not present a higher proportion of fine material and is very close to the material found in the lahar deposits. Figure 10. Gokurakudani gully floor samples (UTM 52S 623353/3624242). Sample D01 was retrieved from a pseudo-terrace on the true right, made of a mixture of material without fabric, which can be attributed to a debris flow deposit. Sample D02 was retrieved from the valley floor, near the gully wall. Samples D03 and D07 were retrieved from small side gully debris cones with a gentle slope of 7.02 degrees. Samples D04 to D06 were all retrieved in the central part of the gully where the GPR transects 7, 8 and 9 were taken. Erosion features and the GPR transects show that the material has been deposited by Newtonian flow with a sub-horizontal layering.
In the lower section of the gully, lahar deposits S09 to S12 are poorly to moderately sorted, and the distributions are very platykurtic (Figure 11). This material contrasts with the wall PDCDs that are at the limit between platykurtic and mesokurtic. If the sampled PDCD material did not display a distinctive difference with the lahar deposits, a drape of fine and wet material was oozing out of the outcrop and its sample showed a significantly finer distribution of material (S08 on Figure 11). This shows that although the wall is stable and not changing from a topographic perspective, its constituents are being selectively washed away. The diminution of fine material from this process may explain why the PDCD material in sample S13 does not present a higher proportion of fine material and is very close to the material found in the lahar deposits. Figure 11. Sample S08: drape over the ignimbrite deposit wall; sample S09 to S11: debris flow deposits near the first confluence with the Tansandani gully; sample S12: potential tail of the debris flow deposits sampled with S09 to S11; sample; S13: ignimbrites from the wall. This set of samples was collected in the vicinity of UTM52 623163/3624407 in the downstream area of the gully.

The Subsurface Structure of the Gully Floor from GPR
The geometry of the gully is changing and so is the material that makes it. The topographical analysis showed areas of erosion, but an overall stability of the gully. We thus investigated with GPR three different locations, two where the topography was stable during the 2004-2018 period and one where a temporary dam created by a wall collapse has raised the gully floor level locally, creating a topographic "step", where finer sediments have accumulated.
The first "stable" location is located at the lower end of the studied gully. The deposit structure ( Figure 12) near the check dam shows a set of units with a weak structure and delineated by a strong reflectivity layer, signaling a difference between the material on the left side of the wedge (from 0 to 20 m on the first radargram of Figure 12) and the rest of the deposited material. Underneath 2 m depth, sub-horizontal layers rich in blocks dominate the lower part, while a set of layers with a series of overlapping layers is located between 0 and 2 m depth. The layers have been numbered based on the deposition sequence (1 to 10 on Figure 12). Layer 1 is a bulk of blocks without structure that lies against the wedge with what was interpreted to be of anthropogenic origin. Upstream, layers 2, 3 and 4 are the connection with the lower unit. Then, located on top of these units, layers 8 to 10 show a downstream prograding pattern, which has the typical morphology of a foreset. On top of this last set of units, a number of short sub-horizontal layers can be found. Figure 11. Sample S08: drape over the ignimbrite deposit wall; sample S09 to S11: debris flow deposits near the first confluence with the Tansandani gully; sample S12: potential tail of the debris flow deposits sampled with S09 to S11; sample; S13: ignimbrites from the wall. This set of samples was collected in the vicinity of UTM52 623163/3624407 in the downstream area of the gully.

The Subsurface Structure of the Gully Floor from GPR
The geometry of the gully is changing and so is the material that makes it. The topographical analysis showed areas of erosion, but an overall stability of the gully. We thus investigated with GPR three different locations, two where the topography was stable during the 2004-2018 period and one where a temporary dam created by a wall collapse has raised the gully floor level locally, creating a topographic "step", where finer sediments have accumulated.
The first "stable" location is located at the lower end of the studied gully. The deposit structure ( Figure 12) near the check dam shows a set of units with a weak structure and delineated by a strong reflectivity layer, signaling a difference between the material on the left side of the wedge (from 0 to 20 m on the first radargram of Figure 12) and the rest of the deposited material. Underneath 2 m depth, sub-horizontal layers rich in blocks dominate the lower part, while a set of layers with a series of overlapping layers is located between 0 and 2 m depth. The layers have been numbered based on the deposition sequence (1 to 10 on Figure 12). Layer 1 is a bulk of blocks without structure that lies against the wedge with what was interpreted to be of anthropogenic origin. Upstream, layers 2, 3 and 4 are the connection with the lower unit. Then, located on top of these units, layers 8 to 10 show a downstream prograding pattern, which has the typical morphology of a foreset. On top of this last set of units, a number of short sub-horizontal layers can be found.  Figures  12 and 13) is preponderant. Large blocks have been marked on the interpretative side of Figure 13, but the vast majority can be seen from the radargrams themselves. These gully floor features are toped in the first few 10 s of centimeters by finer units with blocks less than >1 m depth. To these unstructured deposits, the local sandy material deposited upstream of temporary dams generated by the side gullies collapsing in the gully shows an important contrast. The internal structure shows a set of thin sub-horizontal layers made of sand and coarse sands (samples D04, D05 and D06 in Figure 11). The velocity in the upper part of the material is approximately 0.065 m/ns and 0.035 m/ns in the lower part of the material, richer in blocks.
At the third location, the authors collected a series of transects (transects 7 to 10 in Figure 14), from a "pool" of sandy sediments, which were trapped upstream of a temporary dam created by a small wall collapse, which barred the gully. The radargram shows two sets of layers, one with velocities of 0.065 m/ns and one at 0.035 m/ns ( Figure 14). The unit near the surface is deeper in the center of the gully and it thins outwards towards the gully wall. The limit between the two layers connects to the topography of the side wall collapses, showing that the side wall collapse surface continues underneath the horizontal layer deposited over it. Although the top layer is characterized by a set of horizontal layers, with little hyperbolae, showing only a few blocks, the layer underneath is more chaotic, with numerous hyperbolae showing blocks and sets of less elongated units.  Figures 12 and 13) is preponderant. Large blocks have been marked on the interpretative side of Figure 13, but the vast majority can be seen from the radargrams themselves. These gully floor features are toped in the first few 10 s of centimeters by finer units with blocks less than >1 m depth. To these unstructured deposits, the local sandy material deposited upstream of temporary dams generated by the side gullies collapsing in the gully shows an important contrast. The internal structure shows a set of thin sub-horizontal layers made of sand and coarse sands (samples D04, D05 and D06 in Figure 11). The velocity in the upper part of the material is approximately 0.065 m/ns and 0.035 m/ns in the lower part of the material, richer in blocks.
At the third location, the authors collected a series of transects (transects 7 to 10 in Figure 14), from a "pool" of sandy sediments, which were trapped upstream of a temporary dam created by a small wall collapse, which barred the gully. The radargram shows two sets of layers, one with velocities of 0.065 m/ns and one at 0.035 m/ns ( Figure 14). The unit near the surface is deeper in the center of the gully and it thins outwards towards the gully wall. The limit between the two layers connects to the topography of the side wall collapses, showing that the side wall collapse surface continues underneath the horizontal layer deposited over it. Although the top layer is characterized by a set of horizontal layers, with little hyperbolae, showing only a few blocks, the layer underneath is more chaotic, with numerous hyperbolae showing blocks and sets of less elongated units. Figure 13. GPR transects 7, 8, 9 and 10 of the finer sediments trapped behind a side gully talus, creating a local zone of accumulation and a topographical step. On the two 3D diagrams, one can identify the top layers made of sands and coarse sands with subhorizontal internal structures, covering the talus slope and buried talus slope, in such a way that we know that the valley floor was at a lower level before being raised again. This is further evidence that the erosion of the gully is not a linear process.  . GPR transects 7, 8, 9 and 10 of the finer sediments trapped behind a side gully talus, creating a local zone of accumulation and a topographical step. On the two 3D diagrams, one can identify the top layers made of sands and coarse sands with subhorizontal internal structures, covering the talus slope and buried talus slope, in such a way that we know that the valley floor was at a lower level before being raised again. This is further evidence that the erosion of the gully is not a linear process.
Geosciences 2021, 11,457 19 of 24 Figure 13. GPR transects 7, 8, 9 and 10 of the finer sediments trapped behind a side gully talus, creating a local zone of accumulation and a topographical step. On the two 3D diagrams, one can identify the top layers made of sands and coarse sands with subhorizontal internal structures, covering the talus slope and buried talus slope, in such a way that we know that the valley floor was at a lower level before being raised again. This is further evidence that the erosion of the gully is not a linear process. Figure 14. Transects in the lower section of the Gokurakudani gully, just upstream of the confluence with the Tansandani. Figure 14. Transects in the lower section of the Gokurakudani gully, just upstream of the confluence with the Tansandani.

Topographic Change over Time
In the aftermath of the 1990-1995 eruption of Unzen volcano, material accumulated several tens of meters thick [45], and in this material, the Gokurakudani gully has been eroding vertically, developing an incision mostly 10-20 m deep with steep side walls (in the upper part of the gully, the wall was measured to be 40 m high locally in 2021). During the period 2005-2008, the connection between the upper part of the proto-gully (linked to the dome) and the main gully occurred via regressive erosion. While increasing its length and by then the surface of its watershed, the gully has also been carving its bed vertically and the gully is at present meeting with the pre-eruptive topography at a few locations along the gully [4,46]. It has also been shown in [4] that the tributary gullies were capturing varying portions of the watershed over time, further modifying the amount of water being transferred to the gully and also the location where the water is being conveyed to the main gully. As the gully is becoming deeper, reaching pre-1990-1995 eruption material, the pre-eruptive material located underneath the gully floor has emerged at a few locations, and they provide some general stability, although immediately downstream of the erosion-resistant pre-1990-1995 eruption material, an erosion step can often be found.
In comparison with other volcanoes, however, it is difficult to determine whether this evolution will be representative or marginal. Indeed, compared with lahars reaching more than 10 km at Casita volcano [54] or 20 to 30 km at Merapi Volcano [55], the lahars and the place where they are supposed to occur is more modest (Tsunetaka et al., 2021), and the observed processes and the role of the secondary processes may only be typical of "end of life" lahar gullies.

Grain Size of Different Features in the Gully
The grain size analysis has shown that the gully is not only changing in shape, but also in composition. As the walls are oozing fine material due to slow differentiation of the material, the longer the timespans, the coarser the material becomes. This means that the same deposit, which may seem homogeneous, may actually evolve very differently depending on the topography it is associated with.
Similar lahar deposit variability within one single lahar unit has also been shown by Vazquez et al. [56] (cf. Figure 12). The authors have demonstrated that sediment statistical parameters can vary within one single deposit, and the representativity of sediment grainsize analysis is essential to ponder. Grain size provides an appraisal of a deposit, but not a "full picture", and this spatial variability may further evolve over time, with a change in the composition, which must be typical of coarse unconsolidated deposits. Indeed, as the skeleton of the sediment structure is holding the general shape of the deposit, fine fractions in between may not-for a given concentration-play any role in this skeleton. They can then leave the deposit, without any change in the topography. It is just changing the density of the deposit.

Subsurface Structure: Multiple Processes
The ground penetrating radar data show, at the two lower locations, typical subhorizontal units with terminal overlapping units. A set of units below 100-120 cm deep shows higher relative dielectric permittivity data and reflectors, suggesting that finer material or material with higher moisture contents (and thus with finer grain size) must be present, as has been demonstrated from GPR at other similar volcanoes, in lahar deposits of Semeru volcano [31] and at Merapi volcano in pyroclastic flow deposits [39,40]. At the upstream location, however, a set of horizontal fine-bedded units with fewer blocks, trapped upstream of a local natural dam generated by a wall collapse, show that the general pattern can be broken. More importantly, the finer material thin units, which are semi-horizontal, cannot be directly related to the wall collapse nor to lahars, and this confirms the fact that water is transporting sediments in between lahars, being, with wall collapses and side gully connections, another set of processes occurring in the gully.
Put in perspective, due to the oozing of drapes from the pyroclastic walls, the fine material is thus being exported from the gully, and only coarser material remains, and this secondary process is in turn making the triggering of the lahar more difficult (it is difficult to raise the pore pressure leading to material liquefaction if the pores are very coarse). This oozing process is also a reason for the low popularity of this theory among traditional geologists and sedimentologists. Indeed, as those processes occur, one can start to wonder what is being measured and calculated from outcrops several thousands of years later. Unless the material is welded or fully consolidated during the deposition process, it is most likely that any interpretation (as I have done myself) is biased.

Conclusions
The present volcanic geomorphology study at Unzen volcano has demonstrated that: (1) The left and right walls of the Gokurakudani gully evolve at different rates and that the collapsed profiles and maximum sustained slope angles vary by more than 20 degrees. This difference corresponds to layers in the pyroclastic density current deposits either tilting towards the gully or away from it. (2) The longitudinal profile is overall close to an equilibrium profile, but there is an alternation of stable and of erosive and deposit-dominated zones. These patterns have appeared over time, when the erosion of the material from the last eruption let pre-eruptive rocks emerge. (3) Underneath the surface, the ground penetrating radar confirmed the fact that the gully is undergoing periods of erosion and eventually deposition in the same location. Underneath the present valley floor, there is evidence of gully wall collapses in the valley, which were then filled up. Even if the gully has been eroded during the last 20 years following the eruption, it however does not mean that it is a purely erosive process. Episodes of deposition also occurred. Methodologically, this has implications on the representativity of two sets of LiDAR data taken over time and compared against one another. Any erosion rate calculated from such datasets would not represent the full extent of erosion, and it might be necessary to complete the topographic data with subsurface geophysical investigation. (4) On top of the spatial variability (vertical and horizontal), there is a temporal variability in the material as well. Indeed, as the pyroclastic density current deposits are oozing the fine fraction out, and so without changing the surface topography, it means that (1) over time, there is less fine material available in the gully (notably to trigger lahars) and (2) "precise" grain-size analysis, which is often a pillar in geomorphology or sedimentology, is not a reliable parameter. In other words, the same pyroclastic density current may have a very different grain size distribution at an outcrop, depending on the time at which scientists go and sample it. Even if the deposit does not show any sign of erosion, in the case of unconsolidated coarse and mixed grain deposits, outcrop analysis cannot be performed as scientists would do for fluviatile or estuary data for instance.
Overall, volcanic eruption simplifies the processes occurring on volcanic structures for a period of a decade or so, eventually starting from pyroclastic density currents and then turning towards lahar processes, before the rhythm (frequency + amplitude + coupling with other processes + "non-frequenciable" processes) progressively changes, leaving some space for secondary sedimentary processes to come and play a role in shaping volcanic morphology (Newtonian flows, wall erosion and grain size distribution transformation). Those processes, in turn, are not linear either, and accounting for their complexity is essential to provide an exact picture of the geomorphologic processes occurring on volcanoes. This issue is particularly true because scientists have developed increasingly precise instruments to measure topography notably, but the paradigms that drive them are still the same as the previous generation of researchers, and the assumptions made at a broad scale do not actually stand anymore at a finer scale.