How Can the Morphometric Characteristics and Failure Conditions of a Historic Gully Caused by Intense Rainfall Be Reconstructed?

In January 1980, during exceptional cyclonic rainfall, an atypical landslide, called déboulé, rapidly generated the permanent 700 m-long gully of the Ravine de l’Eglise on an inhabited plateau in Reunion Island (Indian Ocean). Retrieving the initial conditions that led to this historical process is both challenging and necessary for understanding the mechanism of gully incision and providing pointers for improving risk mitigation in relation to this phenomenon. In this study, we reconstruct the preand post-failure topographies using SFM (structure from motion) applied on archive aerial photographs. Based on the comparison of these digital elevation models, we estimate the volume of material eroded to be ca. 0.63 Mm3. Groundwater level increase, part of the triggering mechanism, is hindcast in the catchment of the gully using a lumped hydrological model. This model shows that in only a fortnight the groundwater level probably rose by 36 m, which could have caused a progressive increase in pore pressure and triggered formation of the gully by retrogressive landslides. We test this hypothesis by considering the pre-failure topography and the hindcast groundwater level in a deterministic model based on limit equilibrium equations to explore ground stability. The evolution of ground stability with a rise in the water table shows that the gully may have extended in a headward direction by retrogressive landslides. This is the first quantitative reconstruction of an exceptional historical event affecting the territory of Reunion Island. The methods used to investigate the Ravine de L’Eglise incision thus offer new complementary insights and challenges for understanding the mechanism and the temporality of gully formation.


Introduction
The reconstruction of historical hazardous landslides and erosion processes is a milestone to understand their formation and perform the appropriate hazard assessment [1][2][3][4][5].
Cyclone Hyacinthe lasted from 15 January to 27 January 1980, hitting Reunion Island twice. It drenched the island with rainfall of world-record levels: 6083 mm in 15 days in Commerson and 1144 mm in 12 h in Foc-Foc (c.f. http://www.meteofrance.re/climat/ pluies-annuelles (accessed on 12 December 2021)). This heavy rainfall caused thousands of slope instabilities throughout the island [6,7]. It also led to the formation of gullies on steep slopes, but also on gentle slopes [6][7][8]. The Ravine de l'Eglise, formed during Cyclone Hyacinthe, is the most impressive gully that has formed over the last 40 years. This gully, 700 m long and 50 m wide, crosses the inhabited Grand Ilet plateau. The process resulting in the rapid formation of this gully is locally called déboulé [6,8]. For 40 years, no cyclonic conditions similar to Cyclone Hyacinthe have been recorded and no other déboulé has been Earth 2022, 3 observed on the island. So far, the lack of information has made it difficult to understand the mechanisms involved in its formation. In order to mitigate the risk associated with déboulé in Reunion Island, we need to reconstruct this historical event quantitatively.
A gully can form or extend by hundreds of meters in a single storm by different processes such as overland flow, subsurface piping and/or repeated mass movements [9,10]. The materials eroded during the extension are rapidly removed by streamflow within the gully. Gully initiation is controlled by multiple factors, including morphology, lithology (e.g., geotechnical characteristics), land use and rainfall [9,[11][12][13][14]. Under intense rainfall, surficial and subsurface flows play a significant role in gully formation and expansion. An increase in rainfall intensity causes higher runoff that may intensify gully erosion [15,16]. Groundwater levels, flows and seepage have an impact on gully incision by mass movements and/or internal erosion [17,18].
As this fast and hazardous process cannot be observed in situ and because most of the evidence of its origin has been removed, resorting to physically-based dynamic models has to be considered in order to test the scenarios of its formation. The key factors to simulate such processes are rainfall intensity and duration, morphology and the hydrogeological response of the terrain, and the geotechnical properties of the material [19][20][21][22]. However, the uncertainty of past environmental conditions is high [23], and therefore constraining these factors into a form useful for application in modelling is challenging [15].
Here we tackle this issue, by studying the formation of the Ravine de l'Eglise. The first aim is to retrieve the pre-and post-topography. Since the 20th century, systematic aerial surveys have been conducted in many parts of the island every few years. These archives can be used to build historical ortho-photographs and digital surface model (DSM) and aid detection and quantification of ancient surface changes such as landslides, river incision and beach volumetric changes [24][25][26][27][28][29][30]. It has been demonstrated that a structure from motion and multi-view stereo (SFM-MVS) processing pipeline allows a high-resolution dense point cloud to be automatically and rapidly produced from digitalized archive aerial photographs [24,[29][30][31][32]. As with other photographic sources, some limitations in historical photographs induce errors in the final reconstruction models, such as the quality of the photograph, camera resolution and digitalization, picture overlaps and orientations, and ground control points [24,[33][34][35][36]. Aerial photographs of Reunion Island shot by IGN-France in 1978 and 1984, dates preceding and following Cyclone Hyacinthe (1980), respectively, are available online. These photographs offer excellent potential for retrieving the pre-and post-topography of the Ravine de l'Eglise and for studying the gully's morphology.
The second objective is to assess the hydrological conditions during Cyclone Hyacinthe. Groundwater rise is one of the possible causes of multiple slope failures inducing gully formation, a hypothesis we wish to test in this study. To hindcast historic hydrological conditions, two types of hydrological models can be used: lumped models and physically-based models [35]. Lumped models enable the temporal simulation of a system's various hydrologic processes. The parameters used in a lumped model, such as recorded meteorological data, represent the spatially averaged characteristics of the hydrological system. The advantage of lump models over physically based models is that the conceptual parameterization in the models is simple, and does not require extensive detailed special data to provide satisfying results, in particular for slope stability analysis [37][38][39][40].
Last, considering the reconstructed topography and the hindcast groundwater level, we test the scenario according to which the Ravine de l'Eglise formed by successive headward retreats due to a rise in groundwater. We evaluate the ground stability of the Ravine de l'Eglise during Cyclone Hyacinthe using a deterministic model based on limit equilibrium equations [41]. This simple 2D analysis can help assess the influence of the increase in groundwater on slope stability and verify the shape and depth of the slip surface [42][43][44][45].
The overall aim of this approach is to provide new insights into the understanding of the déboulé phenomenon.

Historical Facts
The gully of the Ravine de l'Eglise incises the northeast zone of the Grand Ilet plateau located in the cirque of Salazie (Reunion Island) over a length of about 700 m (Figures 1 and 2). This gully formed during Cyclone Hyacinthe's heavy rainfalls when 12 days of cumulative rainfall, totalling 5245 mm, was recorded at Grand Ilet (c.f. http://www.meteofrance.re/ climat/pluies-annuelles (accessed on 1 December 2021)). No precise chronology of the event has been described, except that the whole gully formed in the last days of the cyclonic storm. The damage caused by the cyclone is reported in Humbert et al. (1981) [6] and Pinchinot (1984) [7], in which a photo of the déboulé taken just after the cyclone is presented (Figure 1) and the following field observations and interpretations can be read.  Pinchinot (1984) [8], of the déboulé of the Ravine de l'Eglise taken just after Cyclone Hyacinthe. The photo was taken in the field, at a date that is not exactly known. Grand Ilet is located by the red square in the image of Reunion Island to the left. PN, Piton des Neiges; PF, Piton de la Fournaise. Houses on the left side of the picture give the scale. AA' and BB' are scarps of Grand-Ilet landslide.

• Scarp A-A':
Head of the gully: 150 m wide, 50 m deep, head scarp slope > 45 • . Lithology: zeolitized basalt and hard layers, rust-coloured, heavily fractured. Movement: Initiation: possible failure of scarp A-A' due to high pore pressure; Progression: saturated material probably flowed rapidly downslope to the Ravine Roche à Jacquot, digging the terrain over 10 m to 30 m deep.
Springs: possible outlets along scarp A-A' after the first collapse, which may have induced the thinning of the terrain.
• Scarp B-B': Head of the gully: hidden by products of the destabilization of scarp A-A'.

Movement:
Initiation: possible failure of part of scarp B-B' due to high pore pressure; Progression: the material flowed downslope but stayed quite cohesive compared to the material from the failure of scarp A-A'. The photo (Figure 1) shows that the resulting channel was covered by blocs of up to several meters in size, and finer materials. Water was still flowing when the picture was taken. Deposits from intersected flows suggest that several destabilization processes occurred during the rainy episode, in agreement with the description made by Pinchinot (1984) [8].

Geological and Hydrological Setting
The gully formed on the eastern part of the Grand Ilet landslide ( Figure 2). This large compound landslide moves northeast slowly at about 40 to 50 cm·yr −1 in the area of the Ravine de l'Eglise [7,46]. The mean slope of the Grand Ilet landslide is about 10 • . Its surface is characterized by numerous features due to the landslide's movement (e.g., scarps, cracks, depressions). These features are mainly oriented northwest-southeast. The main scarps of the Grand Ilet landslide are at least 20 m high and have a dip that is higher than 40 • .
The gully of the Ravine de l'Eglise crosses ( Figure 2): • A superficial layer of colluvium made of decimetre-to meter-sized blocks of basalt, incorporated in an unconsolidated sandy matrix, less than 10 m thick [47]. • A deeper layer of loose volcaniclastic deposits of the so-called Intermediate Unit of Grand Ilet of about 50 m in thickness, which come from large old failures of the flanks of the Piton des Neiges volcano. These volcaniclastic deposits are made of heterogeneous materials with variable proportions of blocks and mixed materials [48,49]. They are composed of more or less zeolitized basalts. Alteration of the material forms thin layers of grey clays in some areas ( Figure S1). The Intermediate Unit is intensively deformed, in particular due to the movement induced by the slow-moving Grand Ilet landslide, and sensitive to erosion on its steep slopes (>20 • ) (i.e., badlands) [8,43,47]. The Intermediate Unit holds the Grand Ilet aquifer, which has high conductivity and porosity [50].

•
The substratum of the Grand Ilet landslide, an indurated layer of volcaniclastic deposits that crops out in the Ravine Roche à Jacquot, where the river has steep sides [7,8,47].
The gully of the Ravine de l'Eglise drains a part of the Grand Ilet aquifer through the river of the Ravine Roche à Jacquot. Two permanent springs are identified: EG1, located at 922 m NGR, and EG2, located at 929 m NGR. The low water level flow equals 0.8 L/s at EG1 and 0.5 L/s at EG2 [47]. During the rainy season, non-permanent springs appear further upstream the gully.

Method, Data and Modelling Strategy
The first hypothesis put forward by Pinchinot (1984) [8] suggested that the gully resulted from debris-flow initiated by collapses of the Grand Ilet landslide scarps. However, recent observations of different gully formations around the Grand Ilet landslide mainly show regressive destabilizations during rainstorms [6,51]. Without totally rejecting the hypothesis of Pinchinot (1984) [8], it seems more likely that the gully formed by headward retreats [10,18,52]. To test this hypothesis, we first reconstructed the morphometric characteristics and failure conditions and then evaluated the stability of the Grand Ilet plateau under these conditions.

Building a Historic Elevation Model
To retrieve the characteristics of the pre-and post-gully incision topography we applied SFM on archive aerial photography.
Aerial photography campaigns of May 1978 and April 1984 imaged the cirque of Salazie before and after Cyclone Hyacinthe. Digital aerial files are freely available from IGN-France's data portal http://remonterletemp.ing.fr (accessed on 12 February 2021). as JPEG200 files. We selected the 60 pictures of 1978 and the 79 of 1984 that cover the cirque of Salazie (Table 1). The camera calibration parameters of these two series are not known and were solved using SFM software. However, the scanned images were slightly rotated and may cover a wider area than the photo diapositives, inducing disparities in the image coordinates. This biases the estimated camera frame, preventing the finding of a relevant solution for camera calibration parameters [31]. Therefore, the first step for building a 3D point cloud from these photographs consisted of the preparation of the aerial photographs. The images were rotated in order to align the camera's X/Y axis with the line and column of the photo. Then the straightened photographs were cropped to the same origin to remove inscriptions and the black borders of the photographs [31]. Once corrected, the digital photos were loaded in commercial SfM-MVS software developed by Agisoft: Metashape ® .
To build DSM and orthophotography we used the following workflow on Metashape ® : (1) gather the photography of 1978 and 1984 in two respective chunks; (2) on the aerial photos, mask all the moving objects that may induce errors in calculation (e.g., clouds, vehicles); and (3) add the known 11 ground control points (GCP) and 9 check points (CP), spread as regularly as possible over the study area on a stable surface. Because 20% of the cirque is occupied by active landslides, including the Grand Ilet landslide [6,7,46], the controlled points could not be located right beside the déboulé ( Figure 3). Reference points are angular objects (e.g., a rock, a road corner), which can be identified precisely on aerial photographs from 1978 and 1984 as well as on orthophotos from 2015. These points were selected in flat and stable areas (i.e., not on a landslide), and in positions bracketing the 3D domain of interest [53]; i.e., at least pinning all four corners of the area and constraining an elevation span between 266 m and 3066 m above sea level. The coordinates were extracted from a 0.5 m DSM and an orthophotograph of 2015 and encoded in RGR92 UTM40S (EPSG code: 2975). These points were set to provide scale and orientation to the 3D model but also to verify the accuracy of the model built. Then, (4) we aligned photos with high accuracy from both chunks [54,55]; (5) the self-bundle adjustment was computed using four independent parameters (principal distance (f), eccentricity (cx, cy) and radial distortion (K1)); (6) for each separate chunk, a dense point cloud was extracted with ultra-high and aggressive filtering; and (7) we processed the calculation of the DSMs and the orthophotographs at a default resolution of 0.6 m per pixel.
Comparison of the 1978 and 1984 DSMs enabled us to characterize the morphological changes that occurred between 1978 and 1984 and in particular calculate the volume of the déboulé using the DEM of difference (DOD [56]).
Because the two models are not independent, due to the alignment stage, the error estimation of the elevation changes between 1978 and 1984 were obtained by comparing the difference in elevation in a stable area between the two photogrammetric models ( Figures S2 and S3). The accuracy of each model was evaluated separately using the control points and by comparing them to a high-resolution DSM, acquired in an unchanged area by LiDAR in 2015 ( Figure S4).

Hindcast Historic Groundwater Level
No data concerning historic hydrological conditions was available for 1980 except the daily rainfall recorded at the Grand Ilet meteorological station located at ca. 1 km west of the head of the gully ( Figure 2).
To hindcast the water level before and during the 1980 Cyclone Hyacinthe, we use Gardenia ® , an application for lumped hydrological modelling [37,39]. Documentation, technical notes and free software are available at https://www.brgm.fr/en/software/ gardenia-lumped-hydrological-modelling-catchment-basin (accessed on 17 January 2021).
Gardenia ® simulates the main water cycle mechanisms in a basin by applying simple physical laws of flows through reservoirs. These laws are governed by parameters (e.g., transfer time, retention capacity of the soil) determined by calibrations to series of observations. One of the main assumptions considered in this application is that the hydrological behaviour of the basin remains unchanged. We use the rainfall series recorded in Grand Ilet and the measurements of water table level to calibrate the models corresponding to the Ravine de l'Eglise catchment. Rainfall data were recorded at the Grand Ilet meteorological station with a tipping-bucket rain gauge. With Gardenia ® , the daily rainfall series recorded since 2010 was correlated to the daily water table level measured in PZA3, a 100 m-deep piezometer. PZA3 was installed in 2010, 250 m ahead of the main scarp of the gully of the Ravine de l'Eglise ( Figure 2). There, the water level was monitored at 30-min time steps since monitoring started in December 2010. Previous studies [46,50] have proven the presence of a continuous aquifer for the Grand Ilet plateau. Consequently, the piezometer readings are representative of the hydrodynamics of the Ravine de l'Eglise area.
Then, considering this correlation and the rainfall records from 1978, the historic water table levels, especially those under Cyclone Hyacinthe's rainfall, were hindcast.

Modelling Gully Incision by Headward Retreat
Rise in groundwater level increases pore pressure and reduces effective stress and shear stress that may induce slope failure [7]. To test whether or not the rise in groundwater would have triggered retrogressive failures, we assumed that the groundwater level increased gradually and we analysed slope stability accordingly. We evaluate slope stability using limit equilibrium equations. This approach is based on (1) the definition of a failure surface; (2) the division of the failure surface into sections; and (3) the evaluation of the driving and resisting forces and momentum [41]. The ratio between the driving and resisting momentum defines the factor of safety (FoS). A FoS equal to or lower than 1 means that the system is considered to be unstable. Above 1, if the FoS is comprised between 1 and 1.2, the failure is moderately likely to occur and between 1.2 and 1.5 is unlikely to occur. This approach is simple to implement and benefits from considerable experience feedback [41,57].
TALREN ® software (Terrasol) was used to compute slope instabilities. Based on Bishop simplified equations [58][59][60], the tool assumes a rotational slip surface, and divides the rigid body into slices to account for its inhomogeneity and to consider the equilibrium of each slice. Generally, this method yields high factors of safety. It is also less conservative than other methods as it accounts for the effects of interslice forces [59,61].
To achieve computations, three types of data are required: (1) the initial topography and the internal structure of the slope; (2) the unit weight (Yd), the effective cohesion (C') and the effective friction angle (Φ') of each material; and (3) the groundwater level. The topography used is the one retrieved from SFM. The mechanical parameters are of the different geological units crossed by the gully, and are presented in Table 2. These parameters were obtained from laboratory tests carried out in previous studies [62,63]. The groundwater levels are those established with Gardenia ® . Table 2. Mechanical parameters of the Grand Ilet geological units crossed by the gully of the Ravine de l'Eglise. Yd, dry unit weight; C', effective cohesion; Φ', effective friction angle. These parameters were obtained from laboratory tests carried out in previous studies on samples from Roche à Jacquot and Grand Ilet [62,63]. * Laboratory tests made on multiple samples in Reunion for similar materials.

Geological Unit Number of Samples Yd (kN/m 3 ) C (kPa) Φ' ( • )
Detritic cover alluvium (DA) For each groundwater level tested, we performed a stability calculation under static hydraulic conditions. The first simulation consists of testing the stability of the retrieved topography of 1978 by gradually increasing the groundwater level. Once a failure occurs (FoS ≤ 1), the maximum headward retreat was selected and the corresponding topography of 1984 retained.
The longer headward retreat selected was obtained by: 1.
Selecting one of the calculated slip surfaces that has a maximum retreat; 2.
Considering that the failure occurred as far as the selected slip surface, but along the retrieved surface of 1984, which is more realistic than a rotational slip surface; 3.
Verifying afterward the instability of the terrain with this polygonal slip surface with TALREN ® . If the FoS is lower than (or equal to) 1, the topography was modified, and the retrogressive failure assumption validated. On the contrary, if the FoS is higher than 1, we searched for a smaller headward retreat that fulfils the instability conditions.
The following runs were carried out by testing the stability of the successive postfailure topographies obtained after each groundwater level. Downstream from this watershed, headward erosion had already initiated the formation of a small gully (Figures 4 and S6). These features suggest that preferential flow paths above the pre-failure surface of the gully may have existed in 1978. In 1984, the gully head has a badland morphology (Figure 4b). The vegetation has started to grow in some areas only, suggesting that the erosion must still have been active four years after the déboulé.

Pre-and Post-Failure Topographies
Since 1984, no large retreat has been observed in the Ravine de l'Eglise ( Figure S7). Only small landslides and minor erosion processes have occurred along the side slopes ( Figure S8). The gully is therefore considered to be stable.

Morphometry of the Gully
The gully of the Ravine de l'Eglise covers a surface of ca. 7.9 × 10 4 ± 2800 m 2 , is ca. 700 m long and is an average of 50 m wide, from 120 m wide at its head to 40 m at its outlets. These current digital results are perfectly consistent with the observations made after Cyclone Hyacinthe [6,8].
The volume of material evacuated from the Grand Ilet plateau between 1978 and 1984 due to the formation of the gully is estimated from the difference of elevation between the 1 m-resolution digital surface model (DSM) produced of 1984 and 1978 (Figure 5a). The median absolute deviation (MAD) between these two models in a stable area is 1.45 m with a median difference of 0.73 m, which means that the elevation difference is accurate at 0.73 ± 1.45 m ( Figure S2). Furthermore, in the area of the Ravine de l'Eglise, the vegetation was quite short in 1978 and 1984 (Figure 4), which limits errors that are due to vegetation change. The total negative difference (i.e., depletion) in the gully of the Ravine de l'Eglise is ca. 6.3 × 10 5 ± 7.0 × 10 4 m 3 . According to the picture of the déboulé taken just after it occurred (Figure 1 The DoD shows a small area with accumulation on the right side between scarps A-A' and B-B' of the gully where the products of the destabilization of scarp A-A' flowed (Figures 1 and 5a).
The maximum depth of the gully is observed at the toe of the scarps A-A'. Pre-failure and post-failure topographic profiles are presented in Figure 4c. The channel surface (topography of 1984 and 2015) is relatively smooth with a mean slope of 12 • , similar to the initial topography (Figures 5c and S7). The gully has side walls from ca. 30 m to 10 m high, 30 • steep (Figures 4b and S5). Figure 6 shows the hindcast piezometric level of the Grand Ilet aquifer since 1978. The correlation of the measured and the modelled groundwater level from 2010 and 2020 with Gardenia ® is clearly efficient, as attested by a Nash-Sutcliffe model efficiency coefficient (NSE) of 0.959. A model with an estimation error variance equal to zero has an NSE equal to 1 [64]. This good modelling allows the groundwater level to be hindcast at piezometer PZA3 from 1976 to 2010. The groundwater level just before Cyclone Hyacinthe was 1014 m above sea level (a.s.l.); i.e., in the average of those recorded during the rainy season in that piezometer. During Cyclone Hyacinthe, the water table rose progressively in response to the gradual increase in rainfall. The groundwater level may have reached at least 1050 m a.s.l., meaning an increase in groundwater of at least 36 m in only 13 days. Reported on the pre-failure surface (Figure 5c), the maximum piezometric level of Cyclone Hyacinthe crosses the toe of the scarp A-A' that collapsed during Cyclone Hyacinthe. This is consistent with observations from inhabitants during the cyclone, who reported that sources emerged along A-A' after its collapse [6,7]. The hindcast water table north of scarp A-A' was above the ground surface. Such a groundwater level is thought to have led to an increase in pore water pressure and total stress in the whole terrain column between the surface and the less permeable substratum.

The Role of Runoff
Runoff is often considered as one of the main factors influencing gully incision [9,52,[65][66][67]. The rainfall-runoff relationship can be estimated from the following simple relation used from a small drainage basin (i.e., <100 ha) [68]: where Q is the peak runoff rate, I is the average rainfall intensity, A the drainage area of the catchment, and C the runoff coefficient. In our study case, A is the surface of the pre-existing catchment estimated from the DSM of 1978, such as A = 0.067 km 2 . For Reunion and intense rainfall, C is considered equal to 0.9 [69]. Average rainfall intensity is estimated from the daily rainfall, the smallest time step available.
Considering this approximation, we obtained a peak of runoff rate of approximately 0.7 m 3 /s. Such a rate of flow is insufficient to transport the 6.3 × 10 5 m 3 of epiclastic material removed during incision of the Ravine de l'Eglise. The runoff itself is too low to induce gully incision. Groundwater must therefore be strongly involved in the formation of the gully.

Parameter Sensitivity Analyses
A wide range of values are available for the C' effective cohesion and Φ' effective friction angle of the Grand Ilet geological units. The sensitivity analysis aims to constrain the values to be used in the stability analysis.
The sensitivity analysis was conducted for the intermediate unit, the material in which the gully formed. The sensitivity analysis is performed taking into account a dry condition corresponding to the minimum groundwater level, and a wet condition for which groundwater level is at its highest. The factor of safety was calculated for the topography of 1978 using different values of C' and Φ'. In each case, the FoS is calculated by varying a single parameter (C' or Φ'), while keeping the others equal to the values of a basic configuration. The values tested are in the range of values presented in Table 1.
Under dry conditions, the slope (Fos > 1) is stable even with the smaller values of effective cohesion and effective friction angle. Figure 7 shows the results of the sensitivity analysis under wet conditions. Under wet conditions, with a variation of ±50% of C' around its average value (10 kPa), the FoS values vary by 12%. With the same conditions, by varying Φ' by ±28% around its average, the FoS values vary by about 60%. Under wet conditions, the FoS is lower than 1 only if we consider the smaller known values of C' (≤5 kPa) and Φ' (≤23 • ) for the intermediate unit. Table 3 gives the parameters selected for simulation of the gully incision. Considering these parameters, scarp A-A' is at the limit of stability (FoS = 1) and the toe of the slope is unstable (FoS = 0.98) (Figure 8). This result is in agreement with the instability deformations (e.g., cracks, leaning trees) observed in the field along scarp A-A' and at the front of the Ravine de l'Eglise watershed before the gully incision.  Table 3.

Headward Retreat Modelling
Our hypothesis is that the gully of the Ravine de l'Eglise was formed by successive retrogressive collapses caused by the progressive rise in the groundwater level of the Grand Ilet aquifer. Seven successive model iterations were carried out, gradually increasing the water table level from 15 m above the substratum of the Grand Ilet aquifer to ca. 80 m, corresponding, respectively, to dry conditions and saturated conditions that may have occurred during Cyclone Hyacinthe according to the Gardenia ® simulation (Figures 6 and 8).
For scenario t 0 , which is characterized by the initial topography of 1978 and a low groundwater level, the slope toe is stable (FoS > 1.5), but the upstream talus corresponding to scarp A-A' is unstable to conditionally unstable (1 < FoS < 1.2). By increasing the groundwater level by a factor of 2, as has been known to occur in the wet season, a failure takes place downstream from the watershed of the Ravine de L'Eglise and induces a headward retreat of the slope toe of about 100 m (scenario t 1 ). The stability of scarp A-A' remains unchanged. With the same groundwater level, the front of the remaining topography of 1984 is unstable (FoS << 1). These conditions cause the gully to retreat by ca. 150 m (scenario t 2 ). Pursuing the next simulations in this way, with a piezometric level similar to that simulated for the first ten days of Cyclone Hyacinthe, we observe a succession of more or less large retrogressive destabilizations of the slope. The last scenario corresponds to the saturation of almost all of the intermediate unit, as has occurred during the last days of the Cyclone Hyacinthe rainstorm (scenario t 6 ) [8,[42][43][44][45]. For this scenario, we observe that the remaining slope is unstable (FoS << 1) and causes a failure that reaches the current gully head. In these conditions, the post-failure topography is conditionally unstable up to 500 m down the remaining unstable head scarp (FoS = 1.12).
For most of the scenarios, if the failure occurs along a non-circular slip surface instead of a circular one as calculated by Bishop's equations, the FoS is always lower or equal to one, which confirms the headward retreat hypothesis (Table 4). However, for the last scenario (t 6 ), the FoS is equal to 1.15, which means that the slope is conditionally unstable. For this scenario, the remaining slope contains the intermediate unit but also a significant portion of detritic cover. This last unit is quite stable due to a high Φ' and low topographic slope.

Discussion
The results obtained allow us to retrieve the historic conditions of the formation of the gully of the Ravine de l'Eglise and validate the hypothesis that the gully formed with successive headward retreats caused by an increase in groundwater level. Nevertheless, several points require discussion.

Retrieved Topographic Conditions of Pre-and Post-Gully Incision
Although the 1978 and 1984 models are not faultless because of possible errors and distortion, they give indications about the pre-and post-gullying topography.
The construction of ground models based on the co-alignment of sets of archive aerial photographs and ground control points extracted from topographic models yields relatively well-aligned 1 m-resolution models. There is a fair agreement between the elevations estimated for 1978 and 1984 as their difference is lower than 1 m in flat and stable areas (Figures S2 and S3). Elevation discrepancy is sizeable but enables assessment with a 1 m-resolution of the minimum volume of displaced material during the déboulé of 1980, as part of the material may have built up in the gully and vegetation has grown in some places.
Independently, the 1978 and 1984 photogrammetric elevation models are also in fairly good agreement with the LiDAR elevation model (DEM) in 2015, with a difference equal to −0.87 m and 2.77 m, respectively ( Figure S4). A significant part of this difference might be explained by the vegetation [32].
The difference in elevation between the photogrammetric models in a stable and unchanged area can be caused by several factors, such as the quality of the pictures (scratches, dust or stains), camera and digitalization resolution, picture overlaps and ground control points [24,30,35,36,70,71]. Better GCPs will improve the georeferencing and should reduce elevation error. Nonetheless, the procedure used here, exploiting coordinates from orthophotographs and elevation models, can be replicated without fieldwork and offers the possibility of studying historical events in remote areas.
An additional source of limitation in the evaluation of morphological change is the presence of vegetation, which adds noise to photogrammetric models [72]. Part of the uncertainty in the evaluation of morphological change may be caused by changes in vegetation. In the Ravine de l'Eglise area this was not a real limitation, because the vegetation was short and the change in vegetation between 1978 and 1984 was insignificant compared to the size of the gully formed.

Hindcast Historic Groundwater Level
Thanks to observations of the groundwater level over the past ten years, we are able to hindcast a possible groundwater level using a simple and efficient lumped hydrological model. This level might have been different as our simulation was calibrated using observations formed after the gully's incision. Gully incision may have lowered the groundwater levels by providing a shorter drainage path to the outlet [9,12,73]. Hence, a longer drainage path than expected before gully incision should raise the groundwater level. Such conditions increase the probability of observing slope failure in the first days of Cyclone Hyacinthe, but also for the last scenario tested. Belle et al. (2018) [50] showed that a large quantity of rainfall water (up to 250 mm) is adsorbed in shallow soil layers in Grand Ilet plateau. A significant portion of water of this reservoir is evacuated by real evapotranspiration (1500 mm/year) due to the dense tropical broad-leaved vegetation, which lowers aquifer recharge. Although the vegetation has an influence on groundwater recharge, vegetation cover has not changed drastically since 1978, meaning that its effect on estimation of groundwater level can be ignored, especially considering the amount of rainfall and the saturation conditions of the shallow reservoir during Cyclone Hyacinthe.
The peak of piezometric level during Cyclone Hyacinthe could have been reached earlier (i.e., before 29 January) than the model hindcast. Due to its daily time-step, Gardenia ® tends to smooth the results and erase high frequency variations, inducing possible delays to the simulated aquifer response to precipitation. This is all the more plausible in the light of the observations from Belle et al. (2014) [46] showing that the response of the aquifer to rainfall is quite fast.

Headward Retreats and Gully Formation
For Pinchinot (1984) [8] and Humbert et al. (1981) [6], the gully of the Ravine de l'Eglise originated from a forward erosion due to the rapid transportation of the material originating from the successive collapses of scarps A-A' and B-B'; i.e., a debris flow [74]. The stability analysis and field observations confirmed that scarp AA' is at the limit of stability.
From a simple stability model and gradually raising the groundwater level, our results are in favour of the formation of the Ravine de l'Eglise by successive retrogressive landslides. Rise in groundwater may cause a decrease in shear strength due to excess pore water pressure and thus lead to slope failure.
Headward retreats in the formation of gullies is a well-known process, which according to some is related to landslides, especially in the case of intense rainfall, as observed for example in volcanic soil in Ethiopia [16,17,75], in saprolites of various lithologies in Madagascar [13,76,77] or in fluvio-deltaic deposits as in New Zealand [78]. Nevertheless, these processes are often related to internal erosion [13,18,52,76,79,80].
In our study, static hydrological conditions, influence of groundwater flow, material transmissivity or even possible overpressure on the formation of the gully of the Ravine de l'Eglise are not taken into account although they may have played a role by leading to internal erosion, for example, as suggested by Pinchinot (1984) [8].
The scenario we tested is based on many assumptions, and although it is realistic it does not explain why the gully formed in that location and not elsewhere along the Grand Ilet plateau. Spatial variations in hydrogeological properties may control the location of the gully [78]. The epiclastic deposits of Grand Ilet's Intermediate Unit are very heterogeneous poorly consolidated volcaniclastic material. Such heterogeneity may create anisotropic swelling and shrinking with different speeds and intensities related to differences in water content that can promote internal erosion and cause the collapse of the topsoil [76,81,82]. Hence, the location of the déboulé might be explained by local permeability anisotropy of the epiclastic material, which may also have high hydraulic conductivity. Preferential groundwater flows suspected beneath the Ravine de l'Eglise area from the topography of 1978 support this inference, but further analyses are required to confirm it.

Gully Incision and Slope Stability
No significant change in the morphology of the gully has been observed since 1984, whereas other slopes nearby have destabilized following cyclonic rainstorms [51]. By creating a natural drainage of the Grand Ilet aquifer, the gully may have enhanced the stability of the area. The extension of the gully may have stopped because the system has reached the limits of the available catchment of the Ravine de l'Eglise and erosion has become self-limiting [10,81].
The formation of the gully has created a permanent drain of the Grand Ilet aquifer, possibly inducing the general level of the water table to decrease. The water table is the driving force of the Grand Ilet landslide [46]. Paradoxically, the formation of the gully may have slowed down the displacement rate of the Grand Ilet landslide by lowering its water level. Further modelling of the gully effect on the dynamics of the Grand Ilet landslide should be undertaken to study this hypothesis.

Conclusions
The Ravine de L'Eglise is an example of how extremely intense rainfall can give rise to the flash formation of a large permanent gully on a gentle slope through a phenomenon known locally as déboulé. This gully of 700 m in length formed in January 1980 during Cyclone Hyacinthe, a 100-year event.
For the first time, we have built a 1-m 3D model of the gully of the Ravine de l'Eglise based on SFM-MAV methods applied to archive aerial photographs. Not only do preand post-topographic models enable us to assess that material of a volume of 0.63 Mm 3 was eroded during the gully's formation, but they also provide critical information on the geomorphological conditions prior to the failure. The latter is particularly useful for discussing the susceptibility of areas to the formation of déboulé under extreme rainfalls.
Gardenia ® is a simple and efficient hydrological model that enables a realistic historic groundwater level to be hindcast. Thanks to this model, we show that Cyclone Hyacinthe's rainfall caused a gradual rise of at least 36 m in the piezometric level of the Grand Ilet aquifer.
Recovering past-topography and hindcast hydrological conditions of Cyclone Hyacinthe allows us to test the scenario in which the gully formed by successive headward retreats induced by a gradual increase in the groundwater level. This scenario is realistic as it is in agreement with current field observations. Using a simple deterministic model based on limit equilibrium equations, we show that this scenario is probable under hydrostatic conditions. However, field evidence and observations on the retrieved pre-event topography suggest that internal erosion due to groundwater flows and seepage might be implied in the formation of the gully. Further studies will be undertaken to improve understanding of the dynamics and the specific location of déboulé formation.
The methods used for studying the Ravine de L'Eglise formation have created new opportunities for understanding historic events, offering considerable potential in risk assessment.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/earth3010020/s1. Figure S1: Unconsolidated coarse materials of the intermediate unit of Grand Ilet in the gully of the Ravine de l'Eglise. Figure S2: Histogram of elevation difference between 1978 and 1984 in Mare à Vieille Place over a stable area of 7 × 105 m 2 . No major geomorphic changes occurred during Cyclone Hyacinthe in that location. The median (black vertical line) (−0.73 m) and the quantiles at 17% and 83% (grey vertical lines) are respectively −2.64 m and 1.42 m. The RSME is equal to 0.88 m and the median absolute deviation 1.45 m. The error in the estimation of the volume erased during the rainfall of Cyclone Hyacinthe is calculated using the RMSE multiplied by the surface of the area of interest. Figure S3: Estimation of the difference of elevation between the DSM 1984 and the DSM 1978 in Mare à Vieille Place (an area where no major geomorphic change occurred during Cyclone Hyacinthe). a. Distribution of topographic slopes in the area considered. The mean and the mean +/− 1 standard deviations are represented by the black and grey lines respectively. The mean slope is 40 • with a standard deviation of 17 • . b. Mean, standard deviation and RMSE of the elevation difference between 1984 and 1978 as a function of the slope bins. No major differences are observed as a function of the slope, meaning that the error in elevation difference can be uniformly evaluated without considering slope effect, at least for the slopes with angles below 50 • . Figure S4 Funding: This study is part of the RenovRisk-Erosion project. This project is funded by Reunion Island's Regional Council, the European Union (FEDER), the French Government and BRGM.