SWE Modelling of Debris-flow body-sedimentation and Tail-flow Remobilization in a Check-dam Controlled Gully (Unzen Volcano, Japan) using UAV LiDAR, SfM-MVS

On 13th August 2021 at Unzen Volcano (Japan), an 81 mm.hr-1 peak-rainfall (1486 mm in 2 weeks) triggered series of erosion and deposition features in the Tansandani and the Gokurakudani gullies, all adding up to 57,800 m3 of erosion and 39,600 m3 of deposition. Upstream of the Sabo dam located at the exit of the Tansandani Gully, a large deposit has been visually identified as a potential debris-flow front candidate. However, the absence of direct observation leaves some uncertainty on the process that deposited the material and the magnitude of the flow. In the present contribution, the authors are investigating (1) the role of the debris-flow body in constructing the deposit and (2) the role of the tail of the debris-flow in eroding the fresh deposit. To reach this objective, the authors have combined a field investigation with direct observation, UAV (Unmanned Aerial Vehicle) LiDAR (Light Detection And Ranging) and Photogrammetry, and GIS (Geographic Information System) analysis of the field data. From this data, a 2D hydrodynamic simulation and sediment transport models were applied. The results show that the debris-flow ran beyond the first Sabo dam, with part of the material trapped on both sides of the gully. Afterwards, a central channel connected to sub-channels conveyed the diluted flow, most probably < 25 cm with maximum velocities between 4 to 5 m.s-1 at peak flow. Only the debris-flow phase went over the internal shoulders of the Sabo dam. A lobe occupies the top half of the study area and its deposition has been discussed to be related to the sudden widening of the gully, while the lower half is connected to the base-level created by the check dam.


Conceptual background
Debris-flows and their volcanic counterparts' lahars have been defined as mixtures of debris (sediments and bio-debris) and water [1,2] flowing rapidly from mountain and volcanic slopes.Their behavior is still poorly understood with different models balancing equilibrium slopes and sediment concentrations under ideal conditions [3], but the role material plays in the slurry can vary depending on the dynamic, the position, etc. [4,5].As these challenges persist in controlled environments, simulating these flows in the field is still restricted to hazard purposes rather than understanding the underlying mechanics.
In the field, observations of lahars and debris-flow events show further complexities, with debrisflow phases (sediment concentration >80% mass), hyperconcentrated phases (sediment concentration between 60-80%) and a tail with lower sediment concentration and Newtonian-flow behavior (stressstrain linear relation with an intersect at 0 for both variables).Linking deposits and landforms to debrisflow models is therefore a difficult and often impossible task.
These difficulties extend to water flows on volcanoes.Indeed, measuring the flow of water in rills and gullies on stratovolcanoes such as Merapi or Semeru Volcanoes in Indonesia, or Unzen Volcano in Japan (just to cite a few) is particularly challenging, because part of the precipitation gets routed to the subsurface and even when water flows at the surface, scientists can only observe a small portion of it.The porous nature of the bulk pyroclasts also means that the amount lost is important and only during heavy-rainfall events one can start observing concentrated surface water.The issue of measuring and calculating the rainfall/discharge relation is made further complex over time, because (1) the shape and size of water catchments can change rapidly due to erosion [6], and because the grainsize distribution of the material is increasing the further fractions, with the fines being predominantly exported, making the media more porous [7].In arid and semi-arid arheic hydrological systems, similar issues exist, and especially in subhorizontal floodplains where the natural levees controlling the flood change from one year to another, and where there is no base-flow nor predictable water pathways, so that water discharge monitoring and flood preventions are particularly challenging.To curb these issues, the combination of LiDAR data over large-scales and hydrologic simulations have proven to bridge some of the gaps left by the lack of field measurements [8].
Similarly, in volcanic gullies where no base-flow exists and thus no instrumentation has been emplaced (notably due to risk of losing the measurement station due to lahars and other sediment hazards), a similar approach can be used to simulate the water discharge and its relation to water depth at a given location in a gully.
In the present contribution, the authors are hence (1) to define the role of the debris-flow main body in constructing the deposit and (2) to define the role of the tail of the flow (i.e. the end of the debrisflow that has usually a lower sediment concentration. Calibrated against sedimentological evidences from the field of the maximum height reached by an unobserved maximum flow, HRT (High-Resolution Topography) -based hydrodynamic simulation should thus provide estimates of the surface water discharge.Moreover, even if a burst of sediments exist in such system, it is also difficult to prove whether it was deposited by a debris-flow or a flow with lower sediment concentration, even if sedimentologists and geomorphologists provide explanations based on the shape of a deposit.Using the hydrodynamic approach and sediment incipient motion equations, it is also possible to check whether the sediments must have been transported by a lowconcentration flow or by a debris-flow, and it is this approach In order to limit the effect of losses due to groundwater movement, the present research is investigating a gully reach that is bounded by a check-dam, in such a way that the groundwater flow loss can be capped and also at the dam, there will be a net loss due to seepage movement that can be considered to be null, or limited compared to other upstream reaches.

Survey Location in the Tansandani at Mount Unzen in Japan
The Unzen is a stratovolcano located in Southwest Japan on the island of Kyushu where it forms a peninsula in the Ariake-sea.It last erupted between 1991-1995 after almost 300 years of dormancy.The eruption was made famous due to the repeated dome-collapse block-and-ash flows intercalated with lahars, with the largest block-and-ash flow that took the lives of more than 40 journalists and two French scientists who underestimated the ability of pyroclastic surges to decouple from the main flow.Although the post-eruption was characterized by numerous lahars that lead to the construction of the world's most comprehensive system of Sabo dams every developed on a volcano, and where new technologies like remote excavation tools were first tested and used, often decades before it started to be used in other sites overseas (e.g.under the Sumner cliff in Christchurch (New Zealand) in 2011).The frequency of debris-flow has been decreasing over time [9], and the triggering rainfall have increased in duration and intensity.

The August 2023 "debris-flow" in the Tansandani and Gokurakudani
Centered on a peak hourly rainfall of 81 mm.hr-1 on 13 August 2023 (Fig. 3), the period August 7 to 20th recorded cumulated rainfalls in excess of 1486 mm at the Unzendake station [10].The event resulted in extensive sediment movement in both the Tansandani and the Gokurakudani Gully, with deposition upstream the first Sabo structure (where the present research occurred) and erosion in lower Sabo structures.This event contributes to the increase in rainfall threshold necessary for triggering debris-flow at Unzen Volcano, with thresholds inferior to 38 mm before 1999, 39 mm between 2011 and 2016, 60 mm in 2018-2020, and 81 mm in 2021 [10].For hazards and disaster risk purposes, as well as calibration of Sabo-structure, understanding those one-time rogue events, basic data of water discharge, velocity and depths are essential.In order to better understand these events, often considered as secondary, the present contribution aims to investigate upstream of the top Sabo dam what part of the deposit can be attributed to the debris-flow and what part of the landform may have been modified by the tail of the debris-flow and subsequent water flow.Fig. 3 The August triggered debris flow with maximum rainfall of 81 mm,hr -1 from the recording station of Mt Unzendake near the Unzen Fugendake summit (See the summit on Fig. 1-C).

High Resolution Topography acquisition and processing
The High-Resolution Topography was generated using a DJI Mavic-Pro 2® and the original camera, flying at an altitude of 30 to 40 m above the ground, and combining a total of 376 photographs (out of a total of 800 taken and filtered out by location and based on whether the photographs were sharp or blur).The model was not tied to the ground using the traditional tie-points from ground targets, but instead using a LiDAR pointcloud collected by UAV from a DJI Matrice® 300 RTK with a Zenmuse® LiDAR sensor combined with a camera sensor.For the present research a high-density pointcloud was generated and compared to the LiDAR data at 38.4 points per square meter.The 3D model was constructed using the Agisoft Metashape-Pro V2.0 photogrammetry software with the geographic framework in Tokyo Rectangular Coordinate I.The processing steps were (a) the generation of the sparse pointcloud, (b) meshing of the sparse pointcloud, (c) generation of the densepoint cloud.Then the RMSE (Root Mean Square Error) between the SfM-MVS and the LiDAR pointclouds was calculated (eq.1): , and the result [11] is 0.022 m, value that is negligible for the present contribution.

GIS (Geographical Information System) Analysis
The data generated from a combination of UAV LiDAR and Photogrammetry were then converted into a gridded dataset for GIS analysis in the open-source QGIS environment.The gridded data includes an orthophotographs at 3 cm horizontal resolution and a DEM at 10 cm resolution.From this set of raster dataset, the authors have digitized manually the visible blocks at the surface of the deposit.The surface area and the disc-equivalent diameter was then calculated and at the barycentre of each polygon, the local slope was calculated.This process resulted in the construction of a database of 3625 points with the block diameter and the local slope, so that the process necessary for incipient motion could be calculated.

Hydraulic Simulation using the 2D shallow-water equation
The hydraulic simulation was conducted using the 2D shallow-water equation that can be formulated using the conservation of mass and the conservation of momentum with x and y being the two horizontal directions.For an inviscid and incompressible fluid where the horizontal dimensions are dominant against the vertical dimension, the St-Venant equation can be formulated as follows (eq.2,3,4): Where ℎ is the depth, ‫ݑ‬ and ‫ݒ‬ are the depth integrated velocities in x and y (horizontal), g is the gravitational acceleration, z is the topographic surface, ‫ݎ‬ ‫ݔ‬ , ‫ݎ‬ ‫ݕ‬ are the bed friction (in x and y) , ρ is the water density and m is the mass source/drain.‫ݎ‬ ‫ݔݔ‬ and ‫ݎ‬ ‫ݕݔ‬ are the viscous stress and turbulent stress (Reynolds stress) bundled together.The viscous stress is calculated from the fluid kinematic viscosity and the turbulence with a k-epsilon model [12].The bed friction is estimated using the flow density, a coefficient of friction and the depth-averaged equation in x and y (eq. 5 and 6): Finally, C_f is the bed friction that is calculated using the Manning formula formulated with the water height h instead of the typical hydraulic radius for each element of the grid (in our present case an unstructured mesh): ݊ 2 ‫ܥ‬ ݂ = g ℎ 1/3  (Eq.6) In the present experiment, the typical Manning coefficient of 0.06 for rough irregular gravel bed was used [13], and the calculation was performed using finite volume method [14].

The Brownlie critical shear stress and the Peter-Meyer Muller Equation
In order to confirm whether the results of the simulation are pertinent with the fieldwork, and know whether the different clasts have been rather deposited by the debris-flow or the tail of the event (closer to a water flow or banjir), the authors have applied two commonly used equations in gravel transport: (a) the Peter-Meyer Muller sediment transport equation and (b) the Brownlie critical shear stress compared to the shear stress generated at the base of the flow.The different parameters to calculate these equations were obtained from the 2D simulation and gravel size extraction from the orthophotographs.The reader will note that the two equations were generated from laboratory data and they should be viewed as semi-quantitative indicators rather than "final" values.

The Brownlie critical shear stress equation
For his PhD thesis at CalTech, Brownlie defined a critical shear stress equation for gravels sediment transport (eq.7) as follows [15]: , where ߩ ‫ݏ‬ is the sediment and ߩ is the fluid density and ‫ݒ‬ is the kinematic viscosity.This criterion was then compared to the bed shear stress exerted by the simulated fluid and for a given grain-size d, the two shear stresses were compared in order to define whether movement was indeed possible or not.The basal shear stress was defined using Manning's equation (although this equation does not take into account the initial bore).

The Peter-Meyer Muller Equation
One of the commonly used sediment transport equation is the Meyer-Peter and Muller (MPM) equation.The MPM equation ( 1948) is a transformation of the Meyer-Peter formula (eq.8), which gives: where: ρ is the specific mass of water [in metric ton sec -2 /m] ‫ܭ(‬ ܵ ‫ܭ/‬ ‫ݎ‬ ) is an adjustment for the slope S that is adjusted so that only a portion of the total energy loss is due to the grain resistance is used for bed-load motion.
In the present case-study, the local slope varies in a complex manner due to local variation of the bed and the flow.Therefore, averaged values of S at the meter scale were used (by averaging the DEM calculated slope data).Finally, the adjustment slope parameter was formulated from the empirical relation as obtained from laboratory experiments by Peter-Meyer Muller (eq. 9 and 10), also replacing hydraulic radius with flow depth: where: (Eq. 9) And the coefficient Kr was determined by Muller to be 26/d 1/6 where d90 is the size of sediment for which 90% of the material is finer.In the present case, the authors have considered that the large fractions that were digitized consists in the top end of the material and was assimilated to the d90 (this is a necessary simplification and an indicator of whether the largest fraction would move).

Results
The results present successively the field data in a GIS format, then the simulation calibrated against the field data and finally the results of the entrainment equations for the digitized material.

GIS analysis of the survey area: clastic blocks and micro-landforms
The geospatial data shows over a reach of length < 100 m, two main units with a fan-like formation with a general convex (usually concave) surface geometry upstream, and in the lower half a straight channel in the center with a width that approximates 10 m, matching the lower base-level of the sabo dam.On the sides of this channel, two slightly topographically higher benches are characterized with blocks (d=0.16m to d=0.6 m).The upstream fan occurs as the valley width suddenly increase at the location where the Gokurakudani used to meet with the Tansandani gully.The former is now an abandoned channel filled with a debris-flow deposit.The fanning structure shows that a main central lobe with smaller rills, while on both the left and the right side of the fan, larger channels filled with sand are present.The true left-side of the channel is topographically higher than the one on the true right side (Fig. 4).The distribution of the largest blocks (size > d=0.35) shows a concentration on both the right and the left side, of the fan (Fig. 4 A), although the largest blocks were not moved by the flow but fell in the channel due to lateral erosion (from field investigation).However, the largest blocks are concentrated in the area identified as the limit of the lobated fan (Fig. 4 B and C).The Heat map of blocks > 25 cm per 10 m radius shows this dominant distribution (Fig. 4 C).In the central part of the lobe, fewer large blocks were digitized as a concentration of large gravels and other small blocks create a large sheet, that is the apex of the fan.
The lower half of the study area has a fine to coarse-grained central channel and on the sides alignments of blocks creating local ridges following a slightly slanted orientation in the direction of the flow can also be observed.

Hydrodynamic simulation of a flow "matching" the topography and sediments
The geometry of the channel is displaying a complex set of braiding channels before converging back to a single-channel geometry with arms reaching it can be confirmed to be related to the flow dynamics (Fig. 5).The different depths of the streamflow provide indications on the level of connectivity of the different active reach, with the lowest depth segments becoming disconnected first, and eventually connected last.There is a central channel that is doing two full meandering convolution, before straightening in the lower half of the study area (Fig. 5).
Based on the level of connectivity, depth and, the lower-connectivity channels are showing lower depth and velocity and more of the gravels remain while the finer material is not covering the stream area, as it is the case for the main stream.Indeed, the main stream, indeed, although showing the highest velocities (5 m.s -1 ) are characterized post-events by sand and coarse sand deposits.
From the set of experiments and by visual comparison with the orthophotographs and the field photographs and observations, it is a flow between 1 and 2.5 m 3 .s - that best fits the geomorphology of the deposit observed in the surveyed area (please note that this is only the water flowing over the surface).These observations and data, which consider the ground to be fully impermeable only provide an estimation and the values may be underestimated (due to infiltration losses).This complexity and the different uncertainties in dealing with the different connectivity level and height above the main channel (meaning that infiltration rates would vary spatially) can be fortunately alleviated at the Sabo dam, where the flow is consolidated into one and there is virtually no subsurface flow (the dam section does not have any water evacuation gate pierced in the wall), and the relation between the depth of the flow and the velocity at an opening 10 m wide provides a relation where at 4 m.s-1 the depth is about 15 cm, and at 4.5 m.s-1 the depth is <20 cm (Fig. 6).From the combination of the flow simulation that reproduces an inundation extent that fits or resemble the one observed in the field (Fig. 5) and the deposits found at the Sabo dam structure (Fig. 6) it can be then inferred that the tail of the flow or the water flow deposited the material located in the lower 10 m wide section of the dam, but the material on the shoulders have certainly been carried in and sedimented by another process, or a phase of the flow for which the existing sediments and material have been erased by subsequent activity.The presence of fine sand and clastic blocks suggest that the material was transported up on the 50 cm ridge.From Fig. 4-B however, the true left side of the gully shows an accumulation of material and series of subparallel ridges created by clasts on top of the deposit suggesting (a) that either a Newtonian flow with finer material eroded or not deposited or (b) a debrisflow type of flow that would explain how the large clasts have been transported and deposited.On the true right side of the gully, the same explanation is difficult to bring forward, because the banks have been eroded, especially around the location where the flow is creating a reach at the base of the wall (Fig. 5), and it is most likely that some of the large blocks that are on the bed are the results of the bank collapse.Actual visual observations over time have confirmed this idea, and some of the blocks that were on top of the banks are not at the bottom of the gully.

Comparing the critical shear stress, the shear stress and the transport equation
Based on the simulation that matches the last phase of the flow and the position of blocks and the height of deposits, one can state that the clastic blocks on the lobe that is not inundated by the water and the true left side and the shoulders of the Sabo dam are not the results of the diluted tail of the flow.However, for the material that is on the inundated bed, were the blocks transported by the dilute flow, or was it transported by another phase in the process.To solve this issue, the bed shear stress in a diluted water flow was compared to the critical shear stress (Fig. 7 A,B,C,D) and then the critical shear stress for each clastic block was compared to a bedload transport equation, to estimate whether sediments could still be "moving" despite of the fact that incipient motion was not probable (Fig. 7 E).
Using the general slope and the local slope where the material deposited, the critical shear stress is overwhelmingly higher than the shear stress induced by the simulated tail flow, and it confirms that the motion of the blocks is most unlikely under the maximum depth and velocities calculated.Although these equations are used at the reach scale or using "coarse" DEM data, the use of local slope at the meter scale is further emphasizing the results of the blocks not being able to move.Furthermore, the few blocks that may be experiencing incipient motion are not of a given size or located in areas where velocities would have been higher (Fig. 7 A to D).The MPM equation further shows that for all the blocks the potential discharge of material is extremely low, further confirming that the material may not be set in motion by a typical water flow (Fig. 7 E).This pattern is then further confirmed when comparing the ratio of shear stress per critical shear stress for each block and the size 1313 (2024) 012024 IOP Publishing doi:10.1088/1755-1315/1313/1/01202412 of the blocks.Whether the digitized blocks is of a few centimetres to >10s of centimetres, the material remaining was not set in movement by the flow that left the sediment marks and micro-landforms in the channel.

Discussion
The hydrodynamic results performed for the present study reveal the complexity of the flow at the exit of the Tansandani-gully as the flows fans out, because the other gully (the Gokurakudani gully) has been plugged by a debris-flow deposit, and no flow is occurring anymore (another confluence has opened upstream and this arm of the Gokurakudani is 'dead').This fanning area is then constricted back towards the 10 m wide lower-opening of the Sabo dam.However, the side of the gully are in 2022 higher in elevation than the main central gully and it results in complex flows that can be compared to a braided river system.Based on the hydrodynamic simulations (which are not perfect, nor including numerous elements such as erosion and water infiltration -to be added in the next step of this research), the relation between the velocity and the depth of the flow shows that for all flows that match the sediment heights in the channel, heights at the Sabo dam < 0.25 m are most probable, for velocities < 5.0 m.s-1.However, field observations are also confirming the presence of coarse and fine sands on top of the lateral rebar that is 0.5 m high, showing that the flow is reaching such height.There are arguable two possibilities for those.First of all, the material found 50 cm above the lowest point of the Sabo dam is only composed of sand and no gravels nor other size have been found, suggesting that only low energy flow have deposited those sands, and that (1) the flow may have experienced a time-limited surge when large blocks are rolling downstream, as it has been observed at Semeru Volcano at an even larger scale [16]; or (2) the flow being braided upstream, the simulations also show that the flow comes and hit the side of the Sabo dam, so that the material can either splash or accumulate at a higher level than the main central drainage channel, handing up in having a sediment deposit on the rebar, while the flow itself is < 0.5 m deep.The shear stress and critical shear stress data suggest that the material was probably deposited by the debris-flow phase or at least a phase that is rich in clastic blocks, and the tail would have washed away the finer fraction, only leaving the coarse material that it cannot transport.
It is thus most likely that a part of the debris flow stopped where the gully becomes larger, fanning out, but leaving a fan that has a lobated shape, to the difference of alluvial fans.This fanning lobe was then incised by the tail of the flow, which did not carry out the large clastic blocks, but did remove the finer material and transported it downstream the Sabo dam.It is also most likely that a portion of the debris-flow did flow downstream before the formation of the fan, as clast rich deposits are found on the true left side of the gully, but previous observations had already reported the presence of a true-left side bench on the side of the gully, and the concentration of large clasts, may just be the result of finer material erosion.Balancing these results, it is however important to remember that (1) the shear stress and the MCM equations have not been developed for arheic environments and it is difficult to decide precisely how accurate the estimates are.This being stated the comfortable difference between the critical shear stress and the basal shear stress (Fig. 7 B and D) allows for a good error margin.
The present results show that subsequent remobilization can occur at the end of a debris-flow, and it is a call for caution for scientists working from outcrops and deposits, as it has been previously shown at Semeru Volcano (Indonesia), where the relation between the flow and the deposit is actually multiple [16] and outcrops are also not representative of the full deposits [17], and field observations need to be further combined with physical simulations to test whether inferences are indeed in the range of what is "physical" plausible and possible.

Conclusion
The present research has demonstrated that the lobated part of the deposit was most certainly the result of a debris-flow, because the blocks were too big to be transported by a tail or more dilute flow, but important reworking still occurred at the end of the debris-flow, notably with the opening of breaches in the deposit and the fanning of finer sediments.

Fig. 1
Fig. 1 Research Location in South-Japan (A), on Kyushu Island in Nagasaki Prefecture (B) at Unzen Fugendake on Shimabara Peninsula (C), just upstream of the first check-dam structure at the confluence between the Tansandani Gully (North arm on D) and the Gokurakudani Gully (South arm on D).

Fig. 2 (
Fig.2(A) Mount Unzen Fugendake with the dome "hanging" at the summit and on the right side with the Gokurakudani gully (left) and the Tansandani Gully (right) just upstream the research area, where Using airborne LiDAR data of December 2020 and September 2021, Himano et al. (2021) defined erosion volumes of 42,400 m 3 and 15,400 m 3 in the Tansandani and the Gokurakudani respectively, while deposits ranged from 25,900 m 3 and 13,700 m 3 for a single event.It is however unclear what material has been remobilized by the debris-flow only, and what part was remobilized by the tail and/or following water flows.
the bedload rate in underwater weight per unit time and width [metric ton/s/m] ߛ is the specific weight of water [metric ton/m 3 ] ߛ ܵ is the specific weight of sediments [2.45 to 2.65 metric ton/m 3 ] ℎ is the flow depth, replacing the original hydraulic radius [m] S is the slope of energy [m/m] mean particle diameter [m]

Fig. 4
Fig. 4 Block distribution, density and simplified geomorphologic features in the study area (A) imagery and digitized blocks (red color inside and yellow outline); (B) simple geomorphologic interpretation of the surface; (C) Blocks density map of > 0.25 m diameter blocks.

Fig. 5
Fig. 5 Spatial distribution of the surface flow depth (at simulation time 2 minutes, after output and input reached a steady balance).All depth below 0.05 m were visually cut off as they mostly correspond to velocities close to 0, and infiltration and DEM error makes the use of such data aphazard (the small values were kept in figure 04 for comparison purposes on the lower reach of the studied segment).

Fig. 6
Fig. 6 Relation between the velocity of the flow and the flow depth at the Sabo dam, and the geometry of the top of the Sabo dam.(IMWH: Inferred Maximum Water Height most certainly during the "debrisflow phase").

Fig. 7
Fig. 7 Relations between shear stress, critical shear stress with general slope and local slopes (A to D) and the relation between critical shear stress and qb from MPM (E); and (F) the ratio of shear stress to critical shear stress and the block size (d).