# Proxy Data of Surface Water Floods in Rural Areas: Application to the Evaluation of the IRIP Intense Runoff Mapping Method Based on Satellite Remote Sensing and Rainfall Radar

Arnaud Cerbelaud, Pascal Breil, Gwendoline Blanchet, Laure Roupioz, Xavier Briottet

► To cite this version:

Arnaud Cerbelaud, Pascal Breil, Gwendoline Blanchet, Laure Roupioz, Xavier Briottet. Proxy Data of Surface Water Floods in Rural Areas: Application to the Evaluation of the IRIP Intense Runoff Mapping Method Based on Satellite Remote Sensing and Rainfall Radar. Water, 2022, 14 (3), pp.393. 10.3390/w14030393. [hal-03549056](https://hal.science/hal-03549056)

HAL Id: hal-03549056

https://hal.science/hal-03549056

Submitted on 31 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Proxy Data of Surface Water Floods in Rural Areas: Application to the Evaluation of the IRIP Intense Runoff Mapping Method Based on Satellite Remote Sensing and Rainfall Radar

Arnaud Cerbelaud<sup>1,2,3,*</sup> and Pascal Breil<sup>3</sup>

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Cerbelaud, A.; Breil, P.; Blanchet, G.; Roupioz, L.; Briottet, X. Proxy Data of Surface Water Floods in Rural Areas: Application to the Evaluation of the IRIP Intense Runoff Mapping Method Based on Satellite Remote Sensing and Rainfall Radar. *Water* **2022**, *14*, 393. https://doi.org/10.3390/w14030393

Academic Editor: Francesco Viola

Received: 21 December 2021

Accepted: 20 January 2022

Published: 27 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** Along with fluvial floods (FFs), surface water floods (SWFs) caused by extreme overland flow are one of the main flood hazards occurring after heavy rainfall. Using physics-based distributed hydrological models, surface runoff can be simulated from precipitation inputs to investigate regions prone to soil erosion, mudflows or landslides. Geomatics approaches have also been developed to map susceptibility towards intense surface runoff without explicit hydrological modeling or event-based rainfall forcing. However, in order for these methods to be applicable for prevention purposes, they need to be comprehensively evaluated using proxy data of runoff-related impacts following a given event. Here, the IRIP geomatics mapping model, or “Indicator of Intense Pluvial Runoff”, is faced with rainfall radar measurements and damage maps derived from satellite imagery and supervised classification algorithms. Six watersheds in the Aude and Alpes-Maritimes departments in the South of France are investigated over more than 2000 km<fnref n="2" /> of rural areas during two flash-flood events. The results of this study show that the greater the IRIP susceptibility scores, the more SWFs are detected by the remote sensing-based detection algorithm. Proportions of damaged plots become even larger when considering areas which experienced heavier precipitations. A negative relationship between the mean IRIP accumulation scores and the intensity of rainfall is found among damaged plots, confirming that SWFs preferably occur over potentially riskier areas where rainfall is lower. Land use and soil hydraulic conductivity are identified as the most relevant indicators for IRIP to define production areas responsible for downslope deteriorations. Multivariate logistic regression is also used to determine the relative weights of upstream and local topography, uphill production areas and rainfall intensity for explaining SWF occurrence. This work overall confirms the relevance of IRIP methodology while suggesting improvements to its core framework to implement better prevention strategies against SWF-related hazards.

**Keywords:** surface water flood; intense pluvial runoff; susceptibility mapping model; rainfall radar; Sentinel-2; logistic regression

## 1. Introduction

With more extreme precipitations expected in the 21st century due to climate change [1,2], increased attention has to be paid to the understanding and modeling of floods as of now. Comparably to river overflowing, inland flood events occurring outside the vicinity of active waterways have had devastating effects worldwide in the past decades. They have been observed all around Europe, both in urban and rural areas, and not only around the Mediterranean border where severe storms happen on a more frequent basis (see [3] for some examples). In France and

<footnote n="1">¹ ONERA, The French Aerospace Lab., Département Optique et Techniques Associées (DOTA), Université de Toulouse, F-31055 Toulouse, France; laure.roupioz@onera.fr (L.R.); xavier.briottet@onera.fr (X.B.)</footnote>

<footnote n="2">² Centre National d'Etudes Spatiales (CNES), Earth Observation (EO) Lab., 31400 Toulouse, France; gwendoline.blanchet@cnes.fr</footnote>

<footnote n="3">³ Institut National de Recherche pour l'Agriculture, l'Alimentation et l'Environnement (INRAE), Unité de Recherche RiverLy, 69625 Villeurbanne, France; pascal.breil@inrae.fr</footnote>

<footnote>* Correspondence: arnaud.cerbelaud@onera.fr</footnote>
Switzerland, surface water floods have been estimated accountable for respectively 50% and 45% of all flood damage claims in the past years [3,4]. Triggered by more or less intense rainfall, such floods can appear anywhere on a given territory depending on land use, soil type and topography. In this work, the phrase surface water floods (SWFs [3]) will be referred to when considering the various observed impacts of intense overland flow of rainwater on the land surface, as opposed to fluvial floods (FFs) which are confined to the proximity of streams. Another commonly used term is pluvial flood, which is defined as the result of an extreme precipitation event that creates a flood independent of an overflowing water body. Even though some differences might exist within these two terminologies (pluvial floods are often used in the urban context), both will be used as synonyms here. However, this paper will only focus on rural areas, in relation with the implemented methodology (see Section 2 for details).

Contrary to FFs which can be thoroughly studied in instrumented basins thanks to accurate time series of streamflow, data associated with SWFs are much harder to come by. Scarcity of observational data concerning intense overland flow [5–7] assumably stems from multiple combined factors: (i) its potential to appear anywhere in a given area; (ii) its ephemeral nature and (iii) the small spatial scale over which it can materialize [3,8,9]. Indirect measures of SWF impacts are thus used to estimate their occurrence in space and time as well as their severity. We refer to these types of observation as “proxy” data of intense runoff-related impacts. Many proxies exist. Some are obtained from field missions providing ground truth information. Others consist in damage claims contained in insurance records and disaster databases (both agricultural and urban). The latter have been fairly used in the literature to both analyze the spatio-temporal likelihood of SWF occurrence and the validation of pluvial flood prediction models (e.g., [4] in France, but for both FFs and SWFs, [3] in Switzerland and [10] in the Netherlands). Other proxy sources include aerial photos, videos, interviews or local reports. However, none of these sources can be considered completely exhaustive. Insurance records, depending on country legislation and structure, tend to be a very profitable and relevant tool, but suffer from spatial biases towards economically valuable zones [11]. Hence, SWF occurring in uninhabited and non-cultivated areas are not recorded. With the growing availability of satellite imagery at high spatial resolution and high revisit frequency, the detection of SWF footprints on the ground becomes more and more accessible, at least in rural areas [11]. This paves the way for large scale, exhaustive, automatic and repeatable identification of land surface damages in the following days of an extreme weather event, both for operational needs (damage assessment and restoration works) and evaluation of predictive hydrological models.

Among the great variety of hydrological models, some were specifically designed to study and simulate overland flow and intense runoff-related processes such as erosion and sediment yield [12]. Some are physics-based, such as KINEROS [13,14], conceptual, e.g., CEQEROSS [15], or based on a hybrid approach (LISEM [16] and WEPP [17]). The 2D hydraulic FloodAreaHPC model was for instance used for pluvial flood hazard public mapping by Tyrna et al. [18] in Germany. Several geomatics methods have also been developed to map the risks associated with intense runoff. They mostly rely on decision trees depending on topography, land use and soil type, and do not perform any hydrological simulation per se (e.g., STREAM [19]). Other more urban-specific approaches exist. They usually associate susceptibility maps to vulnerability indicators specific to urban networks (population density, buildings distribution, role of various infrastructures). Other flood models were tested in this context, for example on small scale events in Switzerland [20,21]. Bernet et al. [20] assessed the relevance of four models (among which the FLO-2D distributed physics-based model and the MFD geomatics method) to predict SWF occurrence in rural areas while Zischg et al. [21] investigated the performance of the BASEMENT 2D inundation model rather in urbanized areas. Various sources such as event documentations, field observations, aerial photos or insurance claim data were used for performance evaluation.

The IRIP method (French acronym for Indicator of Intense Pluvial Runoff) [12] constitutes a static approach to mapping susceptibility towards the three processes associated with overland flow: production, transfer and accumulation. Without explicit hydrological
modeling, IRIP simply consists in building score maps ranging from 0 (no susceptibility) to 5 (high susceptibility) through the attribution of binary indicators (0 or 1 value) representing various spatially distributed landscape factors such as topography, land use and soil type. Contrary to most distributed hydrological models, IRIP does not rely on an extended set of input parameters and does not require any calibration. Since this method does not simulate runoff dynamics per se, its main goal is rather to allow for the implementation of quick and large-scale preventive strategies against intense runoff hazards. In order to be fully operational, evaluation using observational data is crucial. If spatially comprehensive in situ measurements of overland flow are burdensome to come by, Dehotin et al. [22] managed to implement a field network to monitor rainfall events, soil water content and surface runoff occurrence during a whole year at nine locations in a small rural watershed. Along with multi-factor analysis of runoff occurrence, they confirmed that accumulation and production-prone areas in the IRIP susceptibility maps were associated with a higher surface runoff frequency (see Figure 18 of [22]). However, they did not investigate further the relationship between runoff occurrence and IRIP scores. In other works, risk regulatory zonings of surface runoff and soil erosion as well as post-event survey data of impacts on infrastructures and transport networks such as roads and railways were used to evaluate the IRIP model on several case studies in France [23–25]. IRIP was also applied in the same context to other regions (near Antibes in [26] and the Bayonne area in [27]). Recently, Braud et al. [6] suggested a robust and generic method to evaluate the IRIP susceptibility maps using runoff-related proxy data in the railway context. They developed a workflow based on contingency tables ($\chi^2$ test and false alarm ratios) between areas at risk and observed impacts taking into account hazard levels. However, none of these works could pretend using largely exhaustive and unbiased proxy data on extensive zones to thoroughly evaluate the IRIP method.

The goal of this paper is to use areas impacted by all kinds of extreme rainfall-related disturbances such as gullies, deposits, landslides or mudflows, as proxy data for SWF occurrence. Affected zones were identified using an optical and radar remote sensing-based probabilistic change detection method at 10 m spatial resolution developed in [28,29]. This new repeatable approach was designed specifically to allow for identification of ground deteriorations far from the main hydrographic network. Based on land cadastre plot delineation, land cover mapping [30] and supervised classification methods [31], it allows for the production of maps giving plot-wise probability of damage detection. For this, the machine learning algorithm relies on intra-plot change pixel statistics of vegetation-specific spectral indices. Thanks to these proxy data, this study seeks to comprehensively investigate the quality and limitations of IRIP maps over six watersheds during two extreme Mediterranean events for which the satellite-based method was successfully validated [29]. In order to study the significance of different factors such as land use, pedology or topography on SWF occurrence, the sensitivity to the different IRIP input indicators is also analyzed through contingency statistics and logistic regression, allowing us to suggest improvements to the current methodology. In other works for instance, Niehoff et al. [32] examined the role of land cover conditions on storm-runoff generation using spatially distributed process-based hydrological modeling and found that land-use impacts were specifically relevant during high intensity storms and in small-scale watersheds.

Since the IRIP method only indicates susceptibility to overland flow processes, comparison to actual damages can only be performed if spatially distributed rainfall measurements are available. Indeed, regardless of a site's susceptibility to runoff, the presence or absence of soil degradations following a given event should be conditional to the actual amount and intensity of rainfall received. Bernet et al. [33] and Spekkers et al. [10] worked on characterizing rainfall events leading to SWF on a country scale (Switzerland for the former and the Netherlands for the latter). They both found that high maximum rainfall intensities over short time windows were a major self-sufficient damage predictor. However, Spekkers et al. [10] emphasized the need to add runoff generating mechanisms as explanatory variables in his models to account for the remaining significant amount of unexplained variance. Therefore, in
order for this study to be as comprehensive as possible, analysis of the IRIP intense runoff model is performed considering not only the location of SWF damages but also the spatial distribution of rainfall intensities through radar measurements over the studied events.

To our knowledge, this paper is the first attempt to evaluate a surface runoff susceptibility mapping model using SWF damage proxies derived from satellite remote sensing. Its strength lies in the combination of large scale validated damage maps and rainfall radar measurements to estimate the relevance and suggest improvements to the IRIP method with an operational prospective objective. In addition, the damage proxy maps along with the methodology applied in this paper should be adequate to evaluate the performance and limitations of any other surface runoff mapping model. The paper is organized as follows: Section 2 describes the study areas as well as the multiple materials and methods that were set up in this work. IRIP evaluation results are outlined in Section 3, along with a sensitivity analysis of the model's input indicators for the production map. The relevance of IRIP's accumulation factors and precipitation intensity to predict SWF proxy locations is also investigated through logistic regression in that section. Results are discussed in Section 4 while Section 5 presents the concluding statements and opens up on future methodological developments for both remote sensing of SWF and modeling of pluvial runoff.

## 2. Materials and Methods

### 2.1. Study Areas

This paper focuses on two extreme flash-flooding Mediterranean events. One of them occurred in the Aude department on 15 October 2018, while the other took place in the mountainous region of Alpes-Maritimes on 2–3 October 2020 (“tempête Alex”). Although landscape and topography greatly differ between these two regions, both episodes induced heavy losses from FFs and SWFs due to high intensity precipitations (up to 250 mm locally in 6 h in Aude [34] and 90 mm in only 1 h at Coursegoules rain station in Alpes-Maritimes [35]). Three watersheds in Aude were selected to be modeled by IRIP over a region of interest of more than 1700 km²: Fresquel, Aude-Sou-Lauquet and Aude-centre. For the Alpes-Maritimes study area, the Tinée, Vésubie and Roya basins, representing 1400 km² of steep mountainous woodlands were examined. Selection of these 6 catchments corresponded to the regions of interest over which identification of SWF damages had already been conducted in [29] (see Section 2.2). In this paper, in relation to the work of Cerbelaud et al. [29], evaluation of the IRIP method relied on plot territorial subdivision. For this, the Aude and Alpes-Maritimes official department land cadastre datasets were retrieved. Characteristics of the study areas with regard to SWF detection are detailed in Table 1 (unaffected versus damaged plots, see Section 2.2).

<table><caption>**Table 1.** Characteristics of the study areas depending on satellite detection of SWF damages. Urban areas and plots within 50 m of an active stream were removed from the analysis (Section 2.5).</caption><thead><tr><th rowspan="2"></th><th colspan="2">Aude</th><th colspan="2">Alpes-Maritimes</th></tr><tr><th>Unaffected</th><th>Damaged</th><th>Unaffected</th><th>Damaged</th></tr></thead><tbody><tr><td>Number of plots</td><td>73,023</td><td>12,477</td><td>87,758</td><td>6591</td></tr><tr><td>Total area (km²)</td><td>626</td><td>191</td><td>628</td><td>242</td></tr><tr><td>Plot size (ha): 25–50–75th</td><td>0.14–0.32–0.68</td><td>0.36–0.75–1.62</td><td>0.03–0.10–0.28</td><td>0.09–0.25–0.88</td></tr><tr><td>Slope (°): mean (sd)</td><td>6.8 (6.8)</td><td>3.4 (5.2)</td><td>27.4 (12.0)</td><td>24.3 (15.7)</td></tr><tr><td>Mean—max—sd rainfall over:</td><td></td><td></td><td></td><td></td></tr><tr><td>15 min (mm)</td><td>12–23–4</td><td>13–23–4</td><td>17–39–7</td><td>18–37–7</td></tr><tr><td>60 min (mm)</td><td>30–62–10</td><td>31–62–10</td><td>38–107–18</td><td>41–107–20</td></tr><tr><td>24 h (mm)</td><td>128–233–42</td><td>129–233–41</td><td>209–529–80</td><td>234–529–87</td></tr></tbody></table>

In addition, a map describing the location and topography of all the watersheds in both study sites is provided in Figure 1.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Location and maps of the study areas in the South of France. Top left: Aude (Fresquel, Aude-Centre and Aude-Sou-Lauquet watersheds); Bottom right: Alpes-Maritimes (Tinée, Vésubie and Roya watersheds).</figcaption></figure>

## 2.2. Sentinel Plot-Based Change Detection (SPCD) Method: SWF Damage Proxy Maps

This study is directly subsequent to the works of Cerbelaud et al. [29]. Because actual intense rainwater runoff leading to SWFs is hardly ever observable, they used an original approach aiming at identifying the related footprints on the ground, especially at a distance from the main hydrographic network. The idea was to associate the plot delineation from the official French land cadastre, land use (LU) data and closely timed images from optical (Sentinel-2) and SAR (Sentinel-1) remote sensing to develop a repeatable plot-based change detection method. Identification of damaged areas at 10 m spatial resolution relied on using a machine learning (ML) algorithm, namely the Gaussian Process Classifier (GPC [31]). Explanatory variables consisted of plot-wise statistics depicting various spatial patterns measured from the temporal variations of vegetation-specific spectral indices. The model was run with a single training sample from one study site, while validation was performed on three distinct events in contrasted environments or at a different season. For this, numerous ground truths confirmed by photo-interpretation of very high resolution (VHR) optical imagery (Pléiades satellite, IGN orthophotos) were collected. Interested readers will refer to [11,28,29] for complete details about the method and the associated results.

In this work, the maps produced by this satellite-based method (hereafter referred to as SPCD method, for "Sentinel Plot-based Change Detection" method) will be considered as proxy data of SWFs for the evaluation of IRIP susceptibility maps (Figure 2). Although damage identification of ground truths was not completely exhaustive on the Aude and Alpes-Maritimes validation samples (producer accuracies from 70% to 85%), false discovery rates (or commission errors) inferior to 9% partly ensured that most areas detected by the
SPCD method could be considered with great confidence as proxy data for IRIP maps. Moreover, the very large areas covered (more than 1000 km² for both study sites) could reasonably guarantee consistent and significant evaluation statistics.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Damage maps produced by the SPCD method on the (a) Aude and (b) Alpes-Maritimes events. The color bar indicates the probability given by the GPC for the damaged class (0.5 low confidence, 1.0 very high confidence). Both SWFs and FFs are identified in these maps. The black rectangular box represents the imagery subsets chosen in [29].</figcaption></figure>

## 2.3. *The IRIP Susceptibility Mapping Method*

The present description of the IRIP model relies mostly on the original methodolog-
ical paper from Dehotin and Breil [12]. It is inspired by [6,24,25], but with up-to-date
improvements especially concerning the incorporation of rainfall climatological data in
the methodology. The IRIP method was designed to represent all the processes related to
intense surface runoff occurrence from (i) the generation to (ii) the transfer and (iii) the
accumulation of rainwater flow anywhere on a given territory. It is a geomatics method that
performs generic mapping at the basin scale without any explicit hydrological modeling.
It provides distributed information on the three processes associated with SWF hazards
by defining susceptibility levels from 0 (no risk) to 5 (very high risk). For this, 5 binary
indicators (=1 if favorable to runoff, =0 otherwise) are computed and summed up for
each of the three susceptibility maps, using four types of input data commonly found in
hydrological modeling. As the IRIP model focuses on SWF hazards occurring outside the
river network, it is complementary to FF risk mapping. First of all, in order to initialize
IRIP, a Digital Elevation Model (DEM) of both study sites (IGN BD ALTI® 25 m here) needs
to be processed to determine watershed boundaries and flow directions. The final spatial
resolution of the IRIP maps is constrained by that of the DEM.

### 2.3.1. Production Map

Overland flow of rainwater is more likely to arise over poorly infiltrating surfaces, shallow or saturated soils, and anywhere with little topography. These conditions can result in runoff buildup by infiltration excess and/or saturation excess following heavy rainfall [36]. Land use (LU) data allows discriminating impervious bare surfaces, prone to the production of surface runoff, from infiltrating and rough surfaces that can absorb rainwater and slow down overland flow. Here, each of the 23 classes of the OSO French land cover product [30] available at a 10 m spatial resolution are designated as prone to surface runoff (=1: urban and agricultural, bare grounds mostly) or unfavorable (=0: forests, woody heaths, grasslands) to define a first binary indicator (Figure 3a).
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Framework of the IRIP method for generation of intense runoff susceptibility maps: (a) production, (b) transfer and (c) accumulation. * "Permeability" and "Saturated water content" indicators are computed by comparing the values found respectively for the saturated hydraulic conductivity ($K_s$) and the saturated water content over the first 30 cm of the soil ($\theta_{sat}$) (using ESDAC hydrodynamic characteristics) to the 30-year 1-h and 24-h rainfall intensities (30-year recurrence interval). ** "Slacking crust" and "Erodibility" indicators are computed according to Cerdan et al. [37] using a soil texture triangle and a simple threshold (favorable to runoff if $\ge 3$). *** For topography indicators, slope and (Beven) topographic index are compared to thresholds $\delta_s$ and $\delta_t$ derived using a k-means clustering algorithm ($k = 2$).</figcaption></figure>

Three other indicators for the production map are derived using soil type (ST) information from the European Soil Data Centre (ESDAC) [38]: slacking rate, permeability and saturated water content. Within the ESDAC database, clay, sand and silt proportions (%) for soil texture are available at 500 m resolution for all of Europe. Cerdan et al. [37] methodology is then used to determine slacking (and erodibility for the transfer map, see next subsection) rates using soil texture triangles. Hydrodynamic characteristics of the top 30 cm soil layer such as saturated hydraulic conductivity (cm.day⁻¹) and saturated water content (cm³.cm⁻³) can be directly downloaded at 250 m resolution through ESDAC (thanks to hydraulic pedotransfer functions developed by Toth et al. [39]). In turn, these two attributes are faced respectively with the 1-h (for infiltration excess) and 24-h (saturation excess) mean rainfall intensities corresponding to a 30-year return period for each study area. This period was selected because most substantial disorders seem to start appearing on hillsides at such a recurrence interval. For this, a large 10-year time series (2006–2016) of quantitative precipitation estimations from Panthere radar-based gridded product (see Section 2.4 for more details) was collected over France at a 5-min time step and 1 km² spatial resolution.

Furthermore, the DEM is processed to consider a combination of the local slope and Beven topographic index [40,41] as a 5th indicator related to topography (Figure 3a). For these two variables, a k-means clustering algorithm is applied making use of both the iterative minimum distance [42] and the hill-climbing [43] methods to derive respective thresholds $\delta_s$ and $\delta_t$ and divide the grid values into two classes (=1 if $> \delta$, i.e., favorable to runoff, and =0 if $\le \delta$). This last indicator is then assigned a score of 1 if either the slope or the topographic index is favorable, and 0 if both are not. The Beven topographic index is
used because it reveals areas prone to buildup of surface water, where both the upslope drained area is high (resulting in large potential for overland flow) and the local slope is small (preventing a quick water discharge). It is defined as follows:

$$TI = \ln\left(\frac{a}{\tan b}\right), \qquad (1)$$

where a stands for the upslope drained area and tan b the local slope in radians. Eventually, urban areas identified in the OSO map are systematically given a favorable value of 1 for all four LU/ST indicators for the computation of the production map.

### 2.3.2. Transfer and Accumulation Maps

The first input indicator for the two other maps of susceptibility to transfer and accumulation reflects the need for significant uphill runoff production to induce SWFs. Therefore, for each grid cell, the modal value (i.e., the most represented value) of all IRIP production pixels is estimated among the associated uphill drained area. If this value is greater or equal to 3, the first indicator for both the transfer and the accumulation maps is deemed favorable (=1), and unfavorable otherwise (Figure 3b,c).

Runoff transfer occurs in areas where water can build up speed and flow downwards [44]. Consequences are eroded bedrocks and flattened/uprooted vegetation. A soil erodibility indicator is thus introduced by applying Cerdan et al. [37] method (as for slacking).

Runoff accumulation materializes where water can slow down or even get trapped and where flows can gather and intensify to induce SWF damages such as mud and sediment load deposits. Consequently, the Beven topographic index is adopted alone using the same threshold δt to determine a 0 or 1 contribution to the accumulation map.

Susceptibility maps to transfer and accumulation of surface runoff are further computed through 3 topography-related indicators (Figure 3b,c): drained area, local slope and break of slope. The drained area indicator reflects the need for important uphill contribution, with a 2.5 ha threshold derived from Lagadec [5]. The local slope and break of slope indicators have opposite conditions for transfer and accumulation. Indeed, steep slopes favor transfer while low slopes facilitate water accretion. The same threshold δs was used in both cases. Eventually, convex breaks of slope are deemed favorable to overland flow transfer, as opposed to concave breaks for accumulation.

The final accumulation maps on the six studied watersheds can be seen in Figure 4 (25 m resolution), along with the Sentinel imagery subset chosen for the SPCD proxy maps (Section 2.2 and Figure 2). For each study area, all plots featuring the same maximum IRIP accumulation score X are referenced in the remainder of this paper as “IRIPX” plots.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** IRIP accumulation maps produced on the (a) Aude (Fresquel, Aude-Centre and Aude-Sou-Lauquet watersheds) and (b) Alpes-Maritimes (Tinée, Vésubie and Roya watersheds) areas. Susceptibility scores to surface runoff accumulation range from 0 (no or low susceptibility) to 5 (very high susceptibility).</figcaption></figure>
## 2.4. Rainfall Radar Data

Quantitative precipitation estimations (QPEs) by Panthere radar-based gridded products from Météo-France were retrieved on 48 h windows around each of the two events (Figure 5, [45]). These rainfall estimates rely on the same basic physical principles as all active radar-based instruments. Looking up from the ground, the goal is to use relationships between signal attenuation mostly in the X- or C-band and the drop-size distribution (DSD) to estimate the latter. The Panthere QPE uses a dual-polarized C-band antenna to estimate rainfall at a 1 km² resolution every 5 min. Since cloud physics and weather radar polarimetry are still vast open research activities, radar signal inversion to estimate precipitations still suffers from various sources of uncertainties and constitutes a very challenging topic in hydrometeorology [46,47]. Nevertheless, radar data was shown to be more suitable than rain gauge stations to distinctly capture and characterize extreme precipitation at fine scales and over short durations [48].

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Maximum 60-min quantitative precipitation estimation (QPE) measured by Panthere radar-based gridded products from Météo-France on the (a) Aude and (b) Alpes-Maritimes events.</figcaption></figure>

For this work, maximum rainfall amounts were computed over consecutive periods of 5, 15, 30, 60, 180, 360, 720 and 1440 min on the overall 48 h data windows available, for each 1 km² grid cell of the study areas. It allowed capturing the different maximum rainfall intensities experienced by every 1 km² area over a given period of time, even though each grid cell may have experienced them at different times within the event.

## 2.5. Validation Workflow and Evaluation Metrics

Many methods have been used for the evaluation of flood mapping models using localized flood impacts [1, 18, 20, 24]. The vast majority of them are based on the use of contingency tables and the computation of detection rates and false-alarm ratio. For IRIP, all works to date have used buffers around localized impacts (25 or 50 m) to compensate for uncertainties and approximations in flow paths derived from DEM data (1 or 2 pixels of the BD ALTI® 25 m). In [6], for instance, if at least one IRIP transfer or accumulation pixel with a score ≥ 4 out of 5 was found within the buffer zone around an impact, it was considered as a “hit”. In this study, following the SPCD damage identification operating mode, detection rates and false alarms were defined and measured at the plot scale (see Table 1 for size characteristics). Because of the spatial resolution of Sentinel satellite imagery (10 m), most SWF damages identified by the SPCD method resulted from accumulation mechanisms. Indeed, erosion footprints are on average more elongated (only a few meters wide) and thus often require finer image resolution for identification. Therefore, only IRIP accumulation pixels were considered for evaluation in this work.

When facing IRIP with localized SWF impacts, information on rainfall spatial and temporal distribution is critical to better delineate the evaluation area, as no surface runoff
should be observed from little rainfall accumulation or low-intensity precipitation. Given that the events are often very localized in time, the use of distributed products with short time steps (e.g., 5 min) appears also essential [49].

In this work, we thus propose a quantitative comparison between the IRIP accumulation maps and the SPCD-derived SWF damages taking into account the maximum rainfall intensities over short time steps (Figure 6).

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Evaluation framework of this study. Both SWFs and FFs are identified in the SPCD maps.</figcaption></figure>

Within each cadastre plot of the regions of interest, we selected the maximum IRIP accumulation score ($25 \times 25 \text{ m}^2$ pixels) for evaluation. The maximum rainfall intensity measured over 60 consecutive minutes during a 48 h window around the event was determined at its native spatial scale of $1 \text{ km}^2$. For the SPCD probability maps, a simple 0.5 threshold was set to determine whether the plot belonged to the “damaged” or “unaffected” class. In order to be sure that most damages were due to surface runoff (SWF) and not to river overflowing (FF), a 50 m buffer zone was applied around all active streams using the French BD Carthage 2017 [50]. All plots within this buffer were then removed from the analysis (damaged or not). In addition, urban areas were also removed from the evaluation framework given that SPCD does not perform detection on this particular land cover type (mainly because of its spatial resolution). Other variables used to determine the relative influence of each IRIP input indicator to explain SFW damage location were also aggregated at the plot scale (see Section 3.3 for the associated results).

## 2.6. Multivariate Logistic Regression

Logistic regression is a simple and very efficient type of machine learning algorithm for categorical classification [51]. It consists in modeling the probability of a discrete outcome Y (binary here), called dependent variable, given a set of input data $X = (x_1, ..., x_n)$, called explanatory variables or predictors (Figure 7). For this, the conditional probability of Y being equal to 1 (as opposed to 0), knowing X, is assumed to be given by a sigmoid-type function called logistic function f: $t \rightarrow f(t)$ (see Figure 7). The regression coefficients $\beta_i$ are usually estimated using a maximum likelihood (ML) estimation as loss function. A seemingly analogous measure of the $R^2$ found in linear regression is McFadden's pseudo
R-squared [52], $\rho^2$, which gives a representation of the extent to which the maximized likelihood of the fitted model is close to 1 using the available predictors. It is important to be aware that, empirically, $\rho^2$ usually takes much lower values than those found for $R^2$ in linear regression, as stated in [53] (p. 306): "while the $R^2$ index is a more familiar concept to planner who are experienced in OLS, it is not as well behaved as the $\rho^2$ measure, for ML estimation. Those unfamiliar with $\rho^2$ should be forewarned that its values tend to be considerably lower than those of the $R^2$ index. [...] For example, values of 0.2 to 0.4 for $\rho^2$ represent excellent fit."

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Representation of the functioning of a logistic regression. For each explanatory variable, the odds ratio OR quantifies the change in the odds of the outcome variable Y being equal to 1 (here of a plot being damaged) for a $\Delta_x$ unit increase in the predictor x.</figcaption></figure>

Results are commonly discussed through the computation of odds ratios (ORs). Depending on the nature of predictors, they need to be computed by taking into account the range over which they roughly evolve. In the Results Section 3, OR represent, for each explanatory variable, the change in the odds of Y being equal to 1 (as opposed to 0) when the predictor varies between its 25th and 75th percentile values ($\Delta_x$ variation in Figure 7).

## 3. Results

### 3.1. Evaluation of IRIP Maps through Rainfall Estimation and Damage Identification

Spatial statistics and contingency graphs were produced following the evaluation framework on the two studied events and are summarized in Figure 8. Firstly, 52% of the SPCD-identified damaged plots during the Aude event featured an IRIP accumulation score of at least 4 out of 5 compared to 26% for the rest of the unaffected area. These ratios amounted to 41% versus 13% for Alpes-Maritimes (Figure 8, left panel). Considering the number of plots involved (slightly less than 100,000 for each event), IRIP score differences between damaged and unaffected areas were extremely significant ($\chi^2$ test, $p < 0.0001$). In addition, for the latter event, 75% of damaged plots were tagged with a 3+ score, as opposed to 38% for the remaining lands. As a consequence, these first simple results reveal that IRIP scores do provide a fairly good indication towards susceptibility to SWF occurrence. Nevertheless, IRIP clearly does not give a perfectly discriminating breakdown of affected areas. The Aude region in particular seemed to be overly tagged at medium or high risk of intense overland flow occurrence (IRIP scores $\ge$ 3 for more than 70% of the unaffected plots) in the current IRIP framework, especially regarding level 3. This calls for improvements to the methodology to better discriminate certain susceptibility levels, for instance with a more extended score scale (e.g., from 1 to 10).
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Left panel: Distribution of IRIP accumulation scores among unaffected and SPCD-identified damaged plots in both study events. Right panel: Ratio of SPCD-identified damaged plots among all plots featuring a given IRIP accumulation score. The X-axis represents the minimum threshold considered for rainfall intensity over 60 min (Panthere data). Color bars refer for each scatter point to the number of plots over which the ratio is computed. Upper and lower panels correspond respectively to Aude and Alpes-Maritimes areas. Color legend for IRIP scores applies to all graphs. Lecture note, top right panel: for instance, 25% of IRIP₄ plots (in red) were identified as damaged over the Aude study area. This ratio increases up to 36% when considering all plots that experienced at least a maximum rainfall intensity of 50 mm within 60 min.</figcaption></figure>

Secondly, once again with statistics performed over a great number of observations, a clear and significant segregation appeared in terms of damage proportion among the various IRIP accumulation scores (Figure 8, right panel). For the Aude (Alpes-Maritimes) event, 13% (9%), 25% (17%) and 56% (55%) of respectively all IRIP₃, IRIP₄ and IRIP₅ plots were detected as damaged by the SPCD method. Moreover, these proportions increased for both events when considering the areas that experienced higher rainfall intensities. Indeed, the proportion of damaged IRIP₅ (IRIP₄) plots was found to increase up 20% (10%) for the highest quantiles of precipitation (Figure 8, right panel). For both events, following a rather flat and stable part, the 35–40 mm.h⁻¹ threshold seemed to trigger an increase in the proportion of damaged plots among those indicated at risk by IRIP. Interestingly, this threshold corresponds to a saturated hydraulic conductivity of around 1.10⁻⁵ m.s⁻¹, exceeding most commonly found values except for sandy soils [54]. Such rainfall intensity also corresponds to between 5 and 10-year return periods for Aude and Alpes-Maritimes regions [45]. More interpretations are provided in Sections 4.1 and 4.2. Consequently, IRIP's score methodology can be confirmed as relevant to define susceptibility levels towards intense surface runoff occurrence.

Eventually, a negative relationship between the mean IRIP accumulation scores and rainfall intensities was found over damaged plots (Figure 9). This tendency was present for all time periods of rainfall measurement, but was the strongest for the lowest ones (the 5, 15 and 30-min precipitation intensities). Therefore, SWFs rather occur over more susceptible areas when rainfall is lower. This confirms that the higher the intensity of rainfall, especially during short time windows, the more areas should be prone to intense
overland flow whatever their susceptibility, and thus the less high IRIP scores should be relevant to predict SWF occurrence.

<figure><img src="image_13.png" /></figure>

**Figure 9.** Mean IRIP accumulation scores of SPCD-identified damaged plots (5th, 25th, 50th, 75th and 95th percentiles) as a function of the maximum rainfall amount measured in 5 min (Aude and Alpes-Maritimes, ~19,000 plots). Where rainfall was comparatively lower, SWFs were detected on average over more susceptible areas.

## 3.2. Sensitivity Analysis of IRIP through Land Use and Pedological Characteristics Used for the Production Map

Land use (LU) and pedological characteristics (ST) are crucial factors controlling the generation of overland flow following heavy rainfall. In this section, a sensitivity analysis of IRIP final evaluation scores is performed on the 4 LU/ST-related input binary indicators for the IRIP production map (Figure 3a): land cover, permeability, saturated water content and slacking crust. For this, each of these indicators was removed from the computation of the production map, leading to lower final scores for the accumulation map.

Land cover and permeability were found to be the highest contributors to overall IRIP accumulation scores (Table 2, middle columns). The detection rate, defined as the % of IRIP$_{≥4}$ plots within all SPCD-identified damaged plots, was particularly lower when these two indicators were removed from the production score board (Table 2, left columns). On the contrary, the saturated water content and slacking crust contributed little to the final evaluation scores. This is consistent with the flash-flood dimension of the studied events, where Hortonian runoff is expected to play a major role, along with land use characteristics [32]. The false alarm rate on IRIP$_{≥5}$ plots, i.e., the % of unaffected plots within all IRIP$_{≥5}$ plots, was around 45% for both events, and down to around 25% when considering the areas that experienced the most precipitations. Oddly, false alarm rates decreased when removing certain indicators, in some cases dramatically (e.g., -14% on Aude with no permeability in high rainfall plots) pointing out that these factors were responsible for wrongly indicating high susceptibility in areas that were left unaffected. The negative effects associated with the removal of the most contributing indicators, i.e., land use and permeability, were found to add themselves up and even be amplified (as in Aude) or overlap (Alpes-Maritimes) when removed together. No systematic conclusion could thus be drawn from it. These outcomes call for a more thorough sensitivity analysis of IRIP regarding production factors. Since IRIP currently only considers a binary production contribution depending on a score threshold of 3 out of 5, a better discrimination of the influence of this factor on the final accumulation map should be implemented (see Section 4.3 for suggestions).
<table><caption>**Table 2.** IRIP detection rate (number of IRIP$_{≥4}$ plots among SPCD-identified damaged plots), ratio of overall IRIP$_{≥4}$ plots and false alarm rate on IRIP$_{≥5}$ plots on both study events in the default framework. Differences measured on these ratios are given in the following rows, depending on which indicator was removed from the computation of the production map.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Aude</th><th colspan="3">Alpes-Maritimes</th></tr><tr><th>Detection Rate (%)</th><th>Ratio of IRIP<sub>≥4</sub> Plots Overall (%)</th><th>False Alarm on IRIP<sub>≥5</sub> Plots <sup>1</sup> (%)</th><th>Detection Rate (%)</th><th>Ratio of IRIP<sub>≥4</sub> Plots Overall (%)</th><th>False Alarm on IRIP<sub>≥5</sub> Plots <sup>1</sup> (%)</th></tr></thead><tbody><tr><td>Default framework: All indicators</td><td>52%</td><td>29%</td><td>44% (26%)</td><td>41%</td><td>15%</td><td>45% (26%)</td></tr><tr><td>No land use</td><td>-7%</td><td>-6%</td><td>-1% (-10%)</td><td>-12%</td><td>-9%</td><td>+2% (+2%)</td></tr><tr><td>No permeability</td><td>-7%</td><td>-6%</td><td>-2% (-14%)</td><td>-12%</td><td>-9%</td><td>-5% (-3%)</td></tr><tr><td>No (land use + permeability)</td><td>-30%</td><td>-19%</td><td>-0% (-2%)</td><td>-16%</td><td>-10%</td><td>-5% (-1%)</td></tr><tr><td>No saturated water content</td><td>-2%</td><td>-2%</td><td>+0% (+6%)</td><td>-1%</td><td>-1%</td><td>+0% (+3%)</td></tr><tr><td>No slacking crust</td><td>-0%</td><td>-0%</td><td>-1% (-2%)</td><td>-3%</td><td>-1%</td><td>-2% (+2%)</td></tr></tbody></table>

### 3.3. Relevance of Input Indicators for IRIP Accumulation Maps and Rainfall Estimates to Explain SPCD-Derived SWF Locations Using Multivariate Logistic Regression

In parallel to IRIP's formal evaluation, a multivariate logistic regression was run to try to explain the location of damaged plots as detected by the SPCD method. To this end, the input indicators for the IRIP accumulation map were used before being binarized. Prior to this, correlation statistics were performed on all plots within each study site to avoid collinearity issues within the predictors. The drained area was found correlated (>0.4, *p*-value < 0.001, not shown) with the Beven topographic index and was thus dropped from the set of explanatory variables. Since the model seeks to explain the location of deteriorations that were event-specific and over a large region, rainfall distribution needed to be added among the predictors so as to determine the genuine influence of each IRIP indicator towards SWF occurrence. The list of input variables for the logistic regression was thus eventually composed of:

i. Modal value for IRIP uphill production score, otherwise called mode (from 0 to 5, depending mostly on land cover and soil type, see Figure 3a);

ii. Slope in radians (roughly from 0 to 30° in Aude and to 70° in Alpes-Maritimes);

iii. Break of slope (convex, flat or concave);

iv. Beven topographic index (ranges roughly from 0 to around 25);

v. Maximum rainfall intensity measured over 60 consecutive minutes during a 48 h window around the event.

As the multivariate logistic regression was performed at a plot scale, plot-specific predictors needed to be computed. They consisted of the maximum pixel value contained inside each plot, except for the break of slope for which the presence of a concave break was simply searched for.

Results are displayed in Figure 10. Pseudo R squared values reached 0.12 for Aude and 0.16 for Alpes-Maritimes, which indicated fair to good fit of the model [53]. Substantial conclusions can be drawn from the odds ratios. An OR greater than 1 implies that the odds of a plot being damaged increases when the predictor varies between its 25th and 75th percentile values. Despite the differences in landscape and topography between the two sites, odds ratios were very similar for the Beven topographic index (around 2.20), the modal value for IRIP uphill production score (1.95) and the maximum rainfall intensity (1.25).

Interestingly, the local slope index behaved in opposite manners at the two study sites. For Aude, low local slopes were associated with a greater probability of SWF damage occurrence, as opposed to high slopes for Alpes-Maritimes. This is highly relevant as to the nature of the deteriorations detected in the two sites. For the Aude floods of 2018, most were due to water accretion in low areas within agricultural plots. On the contrary,

<footnote n="1">¹ Numbers in parentheses for the false alarm rate refer to the cases for which a minimum threshold of maximum rainfall intensity (90th percentile) was considered for contingency statistics.</footnote>
during the Alex storm in Alpes-Maritimes, a great number of SWFs were observed through
landslides and mudflows on steep mountain slopes. Either way, the odds ratios for the
maximum local slope predictor were rather low, especially in Alpes-Maritimes (1.05), which
pointed to a smaller weight among the overall influence of topography and rainfall factors
for intense overland flow occurrence. Eventually, the presence of concave breaks of slope
was found either non-significant (p-value > 0.1 on Aude) or with a fairly low odds ratio
(19% increase in odds on Alpes-Maritimes).

<figure><img src="image_14.png" /></figure>

**Figure 10.** Odds ratios representing the change in the odds of being damaged (as opposed to non-affected) when the predictor varies between its 25th and 75th percentile values. Box height indicates the 95% confidence interval for the OR. McFadden pseudo R squared ρ² are given in the figure subtitles. All coefficients are statistically significant at the 1‰ confidence interval, except for the concave break of slope indicator in the Aude regression. For instance, the odds of being damaged are multiplied by 2.3 (or increased by 130%) when the Beven topographic index increases by roughly 5 points on the Aude study event.

Examples of SWF damages identified by the SPCD method are displayed in Figure 11 (Aude) and Figure 12 (Alpes-Maritimes). In Figure 11, the plots received a maximum 46 mm of rainfall within 60 min during the event. Gullies can be observed downhill from high local slopes and within a water pathway materialized by flow direction arrows and high values of the Beven topographic index. Furthermore, IRIP uphill production mode is close to its maximum value of 5. The resulting logit score (probability of damage *p*) from the fitted model is 0.89, indicating the likely presence of an identified SWF, as confirmed by the SPCD damage map. In Figure 12, a mudflow detected by the SPCD method in Alpes-Maritimes was similarly found a high logit score of 0.78, in relation with a 41 mm.h⁻¹ rainfall intensity and incriminating topographical and runoff production factors.

Overall, the conclusions from the logistic regression pointed to a larger influence of the
Beven topographic index and the IRIP uphill production mode in explaining the location of
SWF damages. Consequently, greater weights and larger segmentation should be attributed
to these indicators in comparison with the slope and break of slope factors (see Section 4.3).
<figure><img src="image_15.png" /></figure>

**Figure 11.** Examples of lands affected by SWF with gullies in Aude. Maximum rainfall intensity in 60 min (Panthere data), associated characteristics used by IRIP and resulting logit score on an affected plot (black contour) can be seen in the top left panel. Top left figure: Pléiades image from 3 November 2018. Top right figure: Probability of damage detection by SPCD method. Bottom left figure: Flow direction, local slope (°) and Beven topographic index (DEM-derived). Bottom right figure: IRIP uphill production mode used in the generation of IRIP accumulation map.

<figure><img src="image_16.png" /></figure>

**Figure 12.** Examples of lands affected by SWF with mudslides in Alpes-Maritimes. Maximum rainfall intensity in 60 min (Panthere data), associated characteristics used by IRIP and resulting logit score on an affected plot (black contour) can be seen in the top left panel. Top left figure: IGN orthophotos from 5 October 2020. Top right figure: Probability of damage detection by SPCD method. Bottom left figure: Flow direction, local slope (°) and Beven topographic index (DEM-derived). Bottom right figure: IRIP uphill production mode used in the generation of IRIP accumulation map.
# 4. Discussion

## 4.1. IRIP SWF Evaluation Results

SWF inundation models have usually featured rather low performances when faced with localized impacts. Indeed, countless factors have to be taken into account simultaneously, making it very hard to implement comprehensive evaluation frameworks [10]. In Bernet et al. [20], four different models, from distributed physics-based to “only” DEM-based, were evaluated at the grid cell level. All models had a low evaluation performance with detection rates around 20–30% for half of the case studies, and around 40–50% for the best models on some other events. Only one event was well captured by all the models (see Figure 6 of [20]). The results were limited due to poor observational data based on an external map that indicated ponded water areas but omitted overland flow paths. The very small areas covered in each case study (~0.5 km²) did not allow either for a good rainfall representativeness with radar measurements being thus constant at the 1 km² scale.

Concerning IRIP, results may differ significantly between pixel-based and buffer-based evaluation methodologies. But there is always a trade-off between detection rates and false alarm rates. Lagadec et al. [24] evaluated IRIP maps on the Lézarde catchment (210 km², northern France) using buffers around proxies derived from witness interviews and field expertise. From 41 incidents listed on 21 railway sections from 1995 to 2012, promising probabilities of detection were achieved (73–90%) using IRIP ≥ 4 transfer and accumulation scores. High false alarm ratios were yet obtained (77–92%), which was explained by the presence of hydraulic mitigation structures as well as the lack of data concerning rainfall spatial variability of the associated events. In another work, detection rates between 94 and 100% and false alarm rates from 20 to 26% were reached by comparing IRIP maps to 17 railway sections selected for hydraulic regeneration works (out of 61 sections) on the Bréauté to Fécamp railway line (20 km long) [25]. These great values partly resulted from the limited number of observations combined to a very good field knowledge that allowed to “calibrate” IRIP, using specific input data and parameterizations adapted from the default framework. Braud et al. [6] eventually assessed the IRIP model on an 80 km railway line in Normandy (north of France) using a comprehensive database of runoff-related impacts (59 sections) covering the whole 20th century. To overcome high false alarm ratios, they combined susceptibility levels to vulnerability stakes and mitigation measures in a unique evaluation framework. By doing so, they improved the false alarm ratio from 62% to 28% while keeping a steady detection rate of 86%.

In this work, no adaptation was made to IRIP’s default framework and very generic countrywide data were used for both sites. Fair evaluation scores were obtained (around 50% detection rates with 25% false alarms on high precipitation areas), but on an infinitely more significant and exhaustive database (1000 km² + area with some 90,000 plots for each event) and taking into account rainfall intensity distribution. Furthermore, as SWF evaluation proxies obtained through the SPCD method are derived solely using radiometric information from multispectral images (in addition to not relying on damage claims), they are independent from IRIP’s input factors (which are mostly DEM-based), making this comparison highly meaningful.

Still, rather low results may be explained on IRIP’s part by several rationales: (i) the coarse resolution of IRIP input data (25 m for the DEM, 10 m for land use and 500 m for pedology); (ii) potentially inadequate land use and soil texture susceptibility knowledge; (iii) the presence of mitigation structures such as hedges, terraces and hydraulic networks that are not taken into account in IRIP and (iv) antecedent soil moisture conditions which are overlooked. The very large dataset of SWF occurrences derived from satellite imagery also suffers from uncertainties inherited from the Sentinel plot-based change detection algorithm: (i) land cadastre subdivision and miscellaneous plot sizes; (ii) LU homogeneity hypothesis within plots and (iii) instrument spatial resolution and location accuracy (see [29] discussion for more details). Nevertheless, because IRIP does not aim at exactly predicting SWF damages but rather at giving a broad representation of the areas that are more or less susceptible, all the results of the present study already appear very satisfying.
## 4.2. Space-Time Distribution of Rainfall and SWF Occurrence

In this work, the IRIP method was confirmed as relevant to indicate susceptibility towards surface runoff occurrence considering damage detection by satellite imagery but also rainfall distribution during two extreme events (Figures 8 and 9). The 35–40 mm.h⁻¹ threshold triggering an increase in the proportion of damaged plots among all high IRIP score plots was consistent with the findings of Dehotin et al. [22]. Indeed, from the results of a 1-year field experiment, they found that surface runoff was observed (either from saturation or infiltration excess) in around 40% of measured events when rainfall intensity was lower than 30 mm.h⁻¹, 80% between 30 and 60 mm.h⁻¹, and up to 100% for 60 mm.h⁻¹ and more. However, these observations did not mention whether significant deteriorations of the land surface had resulted from the occurrence of surface runoff. Marra et al. [49] studied the relationship between debris flow and distribution of precipitation in the eastern Italian Alps. They revealed that the event-cumulated rainfall field (from radar and rain gauge network) systematically exhibited a peak corresponding or close to the debris flow initiation points, with an average 33 mm triggering rainfall mean in one hour (see Figure 5 of [49]). Bernet et al. [33] analyzed SWF damage in Switzerland using geolocated insurance claims on a national scale. Time series of precipitation were linked to the location and date of deteriorations over a large and complex region in terms of topography, pedology and geology. With a gridded precipitation product based on both rain gauges and radar data, they showed that a thresholding method based on local quantiles of maximum intensity (90th) and the total sum of event precipitation (98th) allowed for an excellent identification of events leading to significant SWF occurrence. Naturally, a single absolute threshold could not be set for a countrywide prevention strategy.

Considering the return periods associated with the SWF-triggering rainfall intensities (5 to 10 years), the 30-year return period currently used to compute IRIP soil indicators (Ks and θsat) could be lowered. Such thresholds derived from radar-based data yet have to be taken with caution. Indeed, if the spatial distributions obtained from radar are generally more satisfactory than from rain gauge interpolation during heavy storms [48], many studies suggest that radars tend to underestimate absolute rainfall peaks compared with rain gauges up to more than 50% [55–59]. Moreover, it is known that peak intensity biases from radar products can also be important when considering low aggregation timescales (one hour or smaller) in heavy rain ([60] see their Figure 9). Eventually, finer spatial distribution of rainfall could be obtained using downscaling techniques (for example deep learning) or kriging methods with rain gauge data. This would allow better grasping the extent to which the total amount, the timing, or the maximum intensity of precipitation influence SWF occurrence.

## 4.3. Improvement Suggestions to IRIP Accumulation Map Computation

Logistic regression is a very common and practical tool that has already been used in related geosciences studies. It has been achieved in particular for landslide hazard assessment from aerial photographs and field surveys, in order to design classification models [61]. In that work for instance [61], landslide locations were found to be explained by, in descending order of the most influential factors: (i) land use (from Landsat images); (ii) slope; (iii) topographic aspect; (iv) lithology; (v) Normalized Difference Vegetation Index (NDVI, obtained by the authors from SPOT images) and (vi) topographic curvature. Such results can then be used to tune weights for optimal combinations of landscape determinants in models.

It is very important to acknowledge that precipitation intensity location will always remain a high driver of SWF occurrence [10,33] (OR around 1.25 in logistic regression, Figure 10). IRIP is not designed to take it into account (except in terms of climatological likelihood). Therefore, IRIP's goal is not to give a perfect representation of where SWF damages have materialized or are going to in the future. However, no physics-based distributed hydrological model can pretend doing so exactly either, once again because the distribution of rainfall in space and time can never be known at the fine required scale.
Instead, IRIP can serve as a broad benchmark for prevention strategies on a region to country scale. Modifications ought to be tested to match as closely as possible proxy locations available on different large scale events. Consequently, results from general evaluation and sensitivity analysis of IRIP and most of all from logistic regression in Section 3 inspired the following suggestions to the accumulation map. Considering the values of odds ratios, the Beven topographic index could account for around half of the overall score. More importantly, it should be partitioned into more groups to better segregate risk levels. The same goes for the modal value of IRIP uphill production score, which should account for only slightly less than the topographic index. This way, better analysis of the quality of LU/ST information used as input to IRIP could also be made. Therefore, a new suggestion for the IRIP accumulation map would be to lay it on a 0 to 10 score basis, with the Beven index ranging from 0 to 5, the IRIP production mode from 0 to 3 and the slope and break of slope indicators from 0 to 1, unchanged (Figure 13). High susceptibility levels could then be set around a 6 or 7 score for instance, depending on regions and rainfall events.

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Newly suggested computation for the IRIP accumulation map. More thresholds are defined for the uphill production mode as well as for the Beven topographic index. The drained area indicator is dropped. Final score is given between 0 and 10 and provides larger discrimination depending on the production and topographic index ranges.</figcaption></figure>

Eventually, more precise inputs could be used for IRIP. In particular, DEM with finer spatial resolutions could be implemented (like the RGE ALTI® 5 m or LiDAR products ≤ 1 m), or even be replaced by Digital Surface Models (DSM) computed through VHR satellite stereo acquisitions (Pléiades or else CO3D in the coming years).

## 5. Conclusions

SWFs caused by extreme overland flow account for almost half of all flood damage claims every year in France, and most likely in many other regions worldwide. Intense runoff predictive mapping models are yet still very limited, mainly because of the lack of observational data to evaluate them. In this context, a satellite-based classification algorithm called SPCD provides reasonably exhaustive proxy data of SWFs on a very large scale. In this work, the IRIP geomatics approach, a simple score-based susceptibility method, was faced with rainfall radar measurements and damage maps derived from the SPCD method during two flash-flood events. A total of 6 watersheds in the Aude and Alpes-Maritimes departments in the South of France were thus investigated using a highly significant proxy database composed of almost 200,000 plots and over more than 2000 km² of rural areas. Given IRIP's objective of identifying areas at risk of SWFs, this work shows that the susceptibility levels are very relevant on these two distinct study sites. Indeed, the proportions of damaged plots increase accordingly with higher IRIP scores, and even
more so when focusing on where precipitations are the heaviest. Results also corroborate that where rainfall is comparatively lower, SWFs occur on average over riskier areas as defined by IRIP. Using a multivariate logistic regression, precipitation intensity distribution is confirmed as an important driver of SWF occurrence. However, the main predictors of SWF damages location turn out to be the Beven topographic index and the IRIP uphill production mode, which calls for improvements in IRIP's core framework.

Suggested modifications will have to be implemented and tested on various test sites in future studies. Essentially, we are confident that damage proxy maps derived from satellite imagery, including using the SPCD method, can partly solve the great lack in extensive observational data to evaluate and improve any kind of SWF modeling workflow. This work is a first step in the comprehensive evaluation of the quality and limits of the IRIP model. Although IRIP maps will always be conditional to where intense precipitations hit, they may be used as an operational tool to broadly identify areas that require further refined hydraulic expertise and simulation.

Eventually, as most footprints related to runoff transfer processes such as erosion and gullies are less likely to be identified at the 10 m spatial resolution of Sentinel imagery, IRIP transfer maps were not included in this evaluation framework. This should be addressed in future works by developing a semantic segmentation approach. Using a U-net convolutional neural network architecture based on the fusion of Sentinel-2 bi-temporal and VHR post event data (Pléiades), more acute detection capacity could be expected.

**Author Contributions:** Conceptualization, A.C., P.B., G.B., L.R. and X.B.; methodology, A.C. and P.B.; software, A.C. and P.B.; validation, A.C., G.B., L.R., G.B. and X.B.; formal analysis, A.C. and G.B.; investigation, A.C.; resources, A.C., G.B. and L.R.; data curation, A.C. and P.B.; writing—original draft preparation, A.C.; writing—review and editing, A.C.; visualization, A.C.; supervision, X.B.; project administration, G.B., L.R. and P.B.; funding acquisition, L.R., X.B., G.B. and P.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by DGPR/SRNH, grant number 21367400.

**Data Availability Statement:** The data presented in this study is available on request from the corresponding author. The original data is not publicly available under ONERA's internal policy.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## References

1. Kharin, V.V.; Zwiers, F.W.; Zhang, X.; Hegerl, G.C. Changes in temperature and precipitation extremes in the IPCC ensemble of global coupled model simulations. *J. Clim.* **2007**, *20*, 1419–1444. [CrossRef]

2. Wang, X.; Jiang, D.; Lang, X. Future extreme climate changes linked to global warming intensity. *Sci. Bullet.* **2017**, *62*, 1673–1680. [CrossRef]

3. Bernet, D.B.; Prasuhn, V.; Weingartner, R. Surface water floods in Switzerland: What insurance claim records tell us about the damage in space and time. *Nat. Hazards Earth Syst. Sci.* **2017**, *17*, 1659–1682. [CrossRef]

4. Moncoulon, D.; Labat, D.; Ardon, J.; Leblois, E.; Onfroy, T.; Poulard, C.; Aji, S.; Rémy, A.; Quantin, A. Analysis of the French insurance market exposure to floods: A stochastic model combining river overflow and surface runoff. *Nat. Hazards Earth Syst. Sci.* **2014**, *14*, 2469–2485. [CrossRef]

5. Lagadec, L.-R. Développement et Evaluation de la Méthode IRIP de Cartographie Du Ruissellement Intense Pluvial, Application Au Contexte Ferroviaire. Ph.D. Thesis, Doctoral School Earth, Universe, Environment Communauté Université Grenoble Alpes, Saint-Martin-d'Hères, France, 2017; 336p. Available online: https://tel.archives-ouvertes.fr/tel-01778052 (accessed on 27 January 2021).

6. Braud, I.; Lagadec, L.-R.; Moulin, L.; Chazelle, B.; Breil, P. A method to use proxy data of runoff-related impacts for the evaluation of a model mapping intense storm runoff hazard: Application to the railway context. *Nat. Hazards Earth Syst. Sci.* **2020**, *20*, 947–966. [CrossRef]

7. Yu, D.; Coulthard, T.J. Evaluating the importance of catchment hydrological parameters for urban surface water flood modelling using a simple hydro-inundation model. *J. Hydrol.* **2015**, *524*, 385–400. [CrossRef]

8. Spekkers, M.H.; Kok, M.; Clemens, F.H.L.R.; ten Veldhuis, J.A.E. Decision-tree analysis of factors influencing rainfall-related building structure and content damage. *Nat. Hazards Earth Syst. Sci.* **2014**, *14*, 2531–2547. [CrossRef]
9. Steinbrich, A.; Leistert, H.; Weiler, M. Model-based quantification of runoff generation processes at high spatial and temporal resolution. *Environ. Earth Sci.* **2016**, *75*, 1423. [CrossRef]

10. Spekkers, M.H.; Kok, M.; Clemens, F.H.L.R.; ten Veldhuis, J.A.E. A statistical analysis of insurance damage claims related to rainfall extremes. *Hydrol. Earth Syst. Sci.* **2013**, *17*, 913–922. [CrossRef]

11. Cerbelaud, A.; Favro, A.; Roupioz, L.; Blanchet, G.; Briottet, X.; Delvit, J.-M.; Breil, P. Potential of high resolution satellite optical imagery to detect damages following extreme rainfall events. *Houille Blanche* **2020**, *6*, 66–74. [CrossRef]

12. Dehotin, J.; Breil, P. Technical Report of the IRIP Project: Surface Runoff Flood Hazard Mapping; Technical Report; IRSTEA Hydrology-Hydraulic Research Unit: Lyon, France, 2011.

13. Smith, R.E. A kinematic model for surface mine sediment yield. *Trans. ASAE* **1981**, *24*, 1508–1514. [CrossRef]

14. Woolhiser, D.A.; Smith, R.E.; Goodrich, D.C. *KINEROS: A Kinematic Runoff and Erosion Model: Documentation and User Manual*; U.S. Department of Agriculture, Agricultural Research Service, ARS-77: Washington, DC, USA, 1990; p. 130.

15. Duchemin, M.; Lachance, M.; Morin, G.; Lagacé, R. Approche géomatique pour simuler l'érosion hydrique et le transport des sédiments à l'échelle des petits bassins versants. *Water Qual. Res. J. Can.* **2001**, *36*, 435–473. [CrossRef]

16. de Roo, A.P.J.; Wesseling, C.G.; Ritsema, C.J. LISEM: A single-event physically based hydrological and soil erosion model for drainage bassin. I-Theory, input and output. *Hydrol. Process.* **1996**, *10*, 1107–1117. [CrossRef]

17. Forest, G.R.; Lane, L.J. User Requirements. USDA-Water Erosion Prediction Project. (NSERL Report #1); USDA-ARS National Soil Erosion Research Laboratory: West Lafayette, IN, USA, 1987.

18. Tyrna, B.; Assmann, A.; Fritsch, K.; Johann, G. Large-scale high-resolution pluvial flood hazard mapping using the raster-based hydrodynamic two-dimensional model FloodAreaHPC. *J. Flood Risk Manag.* **2016**, *11*, S1024–S1037. [CrossRef]

19. Cerdan, O.; Souchère, V.; Lecomte, V.; Couturier, A.; Le Bissonnais, Y. Incorporating soil surface crusting processes in an expert-based runoff model: Sealing and Transfer by Runoff and Erosion related to Agricultural Management. *Catena* **2002**, *46*, 189–205. [CrossRef]

20. Bernet, D.B.; Zischg, A.P.; Prasuhn, V.; Weingartner, R. Modeling the extent of surface water floods in rural areas: Lessons learned from the application of various uncalibrated models. *Environ. Model. Softw.* **2018**, *109*, 134–151. [CrossRef]

21. Zischg, A.P.; Mosimann, M.; Bernet, D.B.; Röthlisberger, V. Validation of 2D flood models with insurance claims. *J. Hydrol.* **2018**, *557*, 350–361. [CrossRef]

22. Dehotin, J.; Breil, P.; Braud, I.; de Lavenne, A.; Lagouy, M.; Sarrazin, B. Detecting surface runoff location in a small catchment using distributed and simple observation method. *J. Hydrol.* **2015**, *525*, 113–129. [CrossRef]

23. Lagadec, L.-R.; Breil, P.; Chazelle, B.; Braud, I.; Moulin, L. Use of post-event surveys of impacts on railways for the evaluation of the IRIP method for surface runoff mapping. In Proceedings of the E3S Web Conference 3rd European Conference on Flood Risk Management (FLOODrisk 2016), Lyon, France, 17–21 October 2016. [CrossRef]

24. Lagadec, L.-R.; Patrice, P.; Braud, I.; Chazelle, B.; Moulin, L.; Dehotin, J.; Hauchard, E.; Breil, P. Description and evaluation of a surface runoff susceptibility mapping method. *J. Hydrol.* **2016**, *541*, 495–509. [CrossRef]

25. Lagadec, L.-R.; Moulin, L.; Braud, I.; Chazelle, B.; Breil, P. A surface runoff mapping method for optimizing risk assessment on railways. *Saf. Sci.* **2018**, *110*, 253–267. [CrossRef]

26. Breil, P.; Lagadec, L.-R.; Gonzalez-Sosa, E. Inondation pluviale et développement périurbain. In Proceedings of the Novatech 2016, International Conference on Planning and Technologies for Sustainable Management of Water in the City, Lyon, France, 28 June–1 July 2016.

27. Dehotin, J.; Chazelle, B.; Laverne, G.; Hasnaoui, A.; Lambert, L.; Breil, P.; Braud, I. Mise en oeuvre de la méthode de cartographie du ruissellement IRIP pour l'analyse des risques lies aux écoulements sur l'infrastructure ferroviaire. *Houille Blanche* **2015**, *6*, 56–64. [CrossRef]

28. Cerbelaud, A.; Roupioz, L.; Blanchet, G.; Breil, P.; Briottet, X. Supervised classification methods for automatic damage detection caused by heavy rainfall using multitemporal high resolution optical imagery and auxiliary data. *Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.* **2021**, *XLIII-B3-2021*, 693–700. [CrossRef]

29. Cerbelaud, A.; Roupioz, L.; Blanchet, G.; Breil, P.; Briottet, X. A repeatable change detection approach to map extreme storm-related damages caused by intense surface runoff based on optical and SAR remote sensing: Evidence from three case studies in the South of France. *ISPRS J. Photogramm. Remote Sens.* **2021**, *182*, 153–175. [CrossRef]

30. Inglada, J. Contextual version of the OSO 2018 product. Centre d'Expertise Scientifique « CES Occupation des sols ». Available online: http://osr-cesbio.ups-tlse.fr/~|oso/ (accessed on 27 January 2021).

31. Rasmussen, C.E.; Williams, C.K.I. *Gaussian Processes for Machine Learning*; The MIT Press: Cambridge, MA, USA, 2006; ISBN 0-262-2518253-X.

32. Niehoff, D.; Fritsch, U.; Bronstert, A. Land-use impacts on storm-runoff generation: Scenarios of land-use change and simulation of hydrological response in a meso-scale catchment in SW-Germany. *J. Hydrol.* **2002**, *267*, 80–93. [CrossRef]

33. Bernet, D.B.; Trefalt, S.; Martius, O.; Weingartner, R.; Mosimann, M.; Röthlisberger, V.; Zischg, A.P. Characterizing precipitation events leading to surface water flood damage over large regions of complex terrain. *Environ. Res. Lett.* **2019**, *14*, 064010. [CrossRef]

34. Lebouc, L.; Payrastre, O.; Bourgin, F. Reconstitution Des Débits De Pointe Des Crues Du 15 Octobre 2018 Dans Le Bassin De l'Aude. (Convention DGPR-IFSTTAR 2018 No. 2201132931 Du 22 Mai 2018-Action 7 Appui Au SCHAPI); Rapport de Recherche IFSTTAR: Lyon, France, 2019; p. 14.

35. Carrega, P.; Michelot, N. Une catastrophe hors norme d'origine météorologique le 2 octobre 2020 dans les montagnes des Alpes-Maritimes. *Physio-Géo* **2021**, *16*, 1–70. [CrossRef]
36. Horton, R.E. The Role of Infiltration in the Hydrologic Cycle. In Proceedings of the Transactions of the American Geophysics Union, 14th Annual Meeting, Washington, DC, USA, 27–29 April 1933; pp. 446–460.

37. Cerdan, O.; Le Bissonnais, Y.; Souchère, V.; King, C.; Antoni, V.; Surdyk, N.; Dubus, I.; Arrouays, D.; Desprats, J.-F. Guide Méthodologique Pour un Zonage Départemental de L'érosion des Sols. Rapport n° 3: Synthèse et Recommandations Générales (No. BRGM-RP-55104-FR), BRGM-INRA. Available online: http://infoterre.brgm.fr/rapports/RP-55104-FR.pdf (accessed on 27 January 2021).

38. European Soil Database (Distribution Version v2.0). In *Maps of Indicators of Soil Hydraulic Properties for Europe*; European Commission Joint Research Centre: Ispra, Italy, 2016. Available online: https://esdac.jrc.ec.europa.eu/content/maps-indicators-soil-hydraulic-proPERTIES-EUROPE (accessed on 27 January 2021).

39. Toth, B.; Weynants, M.; Nemes, A.; Mako, A.; Bilas, G.; Toth, G. New generation of hydraulic pedotransfer functions for Europe. *Eur. J. Soil Sci.* **2014**, *66*, 226–238. [CrossRef]

40. Beven, K.; Kirkby, M.J. A physically-based variable contributing area model of basin hydrology. *Hydrol. Sci. Bull.* **1979**, *24*, 43–69. [CrossRef]

41. Sørensen, R.; Zinko, U.; Seibert, J. On the calculation of the topographic wetness index: Evaluation of different methods based on field observations. *Hydrol. Earth Syst. Sci.* **2006**, *10*, 101–112. [CrossRef]

42. Forgy, E. Cluster Analysis of multivariate data: Efficiency vs. interpretability of classifications. *Biometrics* **1965**, *21*, 768–780.

43. Rubin, J. Optimal classification into groups: An approach for solving the taxonomy problem. *J. Theor. Biol.* **1967**, *15*, 103–144. [CrossRef]

44. Horton, R.E. Erosional development of streams and their drainage basins; Hydrophysical approach to quantitative morphology. *Geol. Soc. Am. Bull.* **1945**, *56*, 275–370. [CrossRef]

45. Météo-France, Mosaïque RADAR de Cumul Lame d'Eau 5min-PANTHERE. Available online: https://radarsmf.aeris-data.fr/ (accessed on 27 January 2021).

46. Steiner, M.; Smith, J.A.; Uijlenhoet, R. A Microphysical Interpretation of Radar Reflectivity–Rain Rate Relationships. *J. Atmos. Sci.* **2004**, *61*, 1114–1131. [CrossRef]

47. Zhang, G.; Mahale, V.N.; Putnam, B.J.; Youcun, Q.; Qing, C.; Andrew, D.B.; Petar, B.; Dusan, S.Z.; Jidong, G.; Ming, X.; et al. Current Status and Future Challenges of Weather Radar Polarimetry: Bridging the Gap between Radar Meteorology/Hydrology/Engineering and Numerical Weather Prediction. *Adv. Amos. Sci.* **2019**, *36*, 571–588. [CrossRef]

48. Lengfeld, K.; Kirstetter, P.-E.; Fowler, H.J.; Yu, J.; Becker, A.; Flamig, Z.; Gourley, J. Use of radar data for characterizing extreme precipitation at fine scales and short durations. *Environ. Res. Lett.* **2020**, *15*, 085003. [CrossRef]

49. Marra, F.; Nikolopoulos, E.I.; Creutin, J.D.; Borga, M. Space-time organization of debris flows-triggering rainfall and its effect on the identification of the rainfall threshold relationship. *J. Hydrol.* **2016**, *541*, 246–255. [CrossRef]

50. Service D'administration Nationale Des Données Et Référentiels Sur L'eau (Sandre). Cours D'eau-Métropole 2017-BD Carthage. Available online: https://geo.data.gouv.fr/fr/datasets/ee5c709c9b7ff928ab2529b79ce6e879c4de6950 (accessed on 27 January 2021).

51. Berkson, J. Application of the Logistic Function to Bio-Assay. *J. Am. Stat. Assoc.* **1944**, *39*, 357–365. [CrossRef]

52. McFadden, D. Conditional logit analysis of qualitative choice behavior. In *Frontiers in Econometrics*; Zarembka, P., Ed.; Academic Press: New York, NY, USA, 1974; pp. 105–142.

53. McFadden, D. Quantitative Methods for Analyzing Travel Behaviour on Individuals: Some Recent Developments. In *Behavioural Travel Modelling*; Hensher, D., Stopher, P., Eds.; Groom Helm London: London, UK, 1979; p. 306.

54. U.S. Department of Agriculture; Natural Resources Conservation Service. National Soil Survey 430-VI. Available online: http://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/ref/?cid=nrcs142p2_054242 (accessed on 27 January 2021).

55. Smith, J.A.; Seo, D.J.; Baeck, M.L.; Hudlow, M.D. An intercomparison study of NEXRAD precipitation estimates. *Water Resour. Res.* **1996**, *32*, 2035–2045. [CrossRef]

56. Smith, J.A.; Baeck, M.L.; Villarini, G.; Welty, C.; Miller, A.J.; Krajewski, W.F. Analyses of a long-term, high-resolution radar rainfall data set for the Baltimore metropolitan region. *Water Resour. Res.* **2012**, *48*, W04504. [CrossRef]

57. Overeem, A.; Buishand, T.A.; Holleman, I. Extreme rainfall analysis and estimation of depth-duration-frequency curves using weather radar. *Water Resour. Res.* **2009**, *45*, W10424. [CrossRef]

58. Peleg, N.; Marra, F.; Fatichi, S.; Paschalis, A.; Molnar, P.; Burlando, P. Spatial variability of extreme rainfall at radar subpixel scale. *J. Hydrol.* **2018**, *556*, 922–933. [CrossRef]

59. Goudenhoofdt, E.; Delobbe, L.; Willems, P. Regional frequency analysis of extreme rainfall in Belgium based on radar estimates. *Hydrol. Earth Syst. Sci.* **2017**, *21*, 5385–5399. [CrossRef]

60. Schleiss, M.; Olsson, J.; Berg, P.; Niemi, T.; Kokkonen, T.; Thorndahl, S.; Nielsen, R.; Nielsen, J.E.; Bozhinova, D.; Pulkkinen, S. The accuracy of weather radar in heavy rain: A comparative study for Denmark, the Netherlands, Finland and Sweden. *Hydrol. Earth Syst. Sci.* **2020**, *24*, 3157–3188. [CrossRef]

61. Lee, S. Application of logistic regression model and its validation for landslide susceptibility mapping using GIS and remote sensing data. *Int. J. Remote Sens.* **2005**, *26*, 1477–1491. [CrossRef]