23 Grdr Migrations, Citoyenneté, Développement, 93100 Montreuil, France; gwenaelle.dejacquelot@grdr.org (G.D.); eugenia.gallese@grdr.org (E.G.)

24 OMVS Organisation Pour la Mise en Valeur du Bassin du Fleuve Sénégal, Dakar 12500, Senegal;
diakitemm2@gmail.com

25 FPFD, Fédération des Paysans du Fouta Djallon BP 52 PITA, Timbi Madina 7311, Guinea;
kourahoye@yahoo.fr

26 OMVG Organisation pour la Mise en Valeur du Bassin du Fleuve Gambie, Dakar 12500, Senegal;
lkonate@omvg.sn

27 Institut de Recherches Géologiques et Minières, Centre de Recherches Hydrologiques, Yaoundé BP 4110,
Cameroon; bnnomo@gmail.com

28 IRD, Cirad, INRAE, Montpellier SupAgro, UMR Eco&Sols, University of Montpellier, 34196 Montpellier,
France; didier.orange@ird.fr

29 SENASOL, Service National des Sols, Conakry 001, Guinea; mouctaryaya3@gmail.com

30 UMR IGE Univ Grenoble/Alpes, IRD, CNRS, Grenoble INP, Université Grenoble Alpes, 38058 Grenoble,
France; jean-pierre.vandervaere@univ-grenoble-alpes.fr

* Correspondence: luc.descroix@ird.fr; Tel.: +33-678-92-06-84

Received: 31 July 2020; Accepted: 10 September 2020; Published: 22 October 2020

**Abstract:** A large share of surface water resources in Sahelian countries originates from Guinea's Fouta Djallon highlands, earning the area the name of "the water tower of West Africa". This paper aims to investigate the recent dynamics of the Fouta Djallon's hydrological functioning. The evolution of the runoff and depletion coefficients are analyzed as well as their correlations with the rainfall and vegetation cover. The latter is described at three different space scales and with different methods. Twenty-five years after the end of the 1968–1993 major drought, annual discharges continue to slowly increase, nearly reaching a long-term average, as natural reservoirs which emptied to sustain streamflows during the drought have been replenishing since the 1990s, explaining the slow increase in discharges. However, another important trend has been detected since the beginning of the drought, i.e., the increase in the depletion coefficient of most of the Fouta Djallon upper basins, as a consequence of the reduction in the soil water-holding capacity. After confirming the pertinence and significance of this increase and subsequent decrease in the depletion coefficient, this paper identifies the factors possibly linked with the basins' storage capacity trends. The densely populated areas of the summit plateau are also shown to be the ones where vegetation cover is not threatened and where ecological intensification of rural activities is ancient.

**Keywords:** Fouta Djallon; water tower; depletion (or recession) coefficient; runoff coefficient; soil water holding capacity; basement; sandstone; Fula society

# 1. Introduction: Statements and State of the Art

The Fouta Djallon receives some of the greatest rainfall in West Africa. Although rainfall decreases inland from the coast, where Conakry receives more than 4000 mm yearly, to the mountains, the highlands, which separate the coastal basins from the great "sudano-sahelian" basins of Senegal and Niger Rivers, have annual precipitation ranging from 1800 to 2300 mm. On the northern slope, annual rainfall ranges from 1300 to 2000 mm. This is sufficient to ensure significant streamflows on both the coastal and continental hillslopes of the Fouta Djallon and the "Guinean Dorsale", its natural extension southeastward [1,2]. At the continental scale, the Fouta Djallon is one of the two great "water towers" of Northern sub-Saharan Africa (Figure 1), the other one being in Ethiopia [3,4].
<figure><img src="image_44.png" /></figure>

37/44

<figure><img src="image_45.png" /><figcaption>(b)</figcaption></figure>

Figure 1. *Cont.*
<figure><img src="image_46.png" /><figcaption>**Figure 1.** The two main natural water towers of northern sub-Saharan Africa and their 1951–2000 mean annual discharges; (a) comparison of the water towers (total discharges of French rivers is given for comparison purposes); (b) the Guinean water tower; (c) the Ethiopian water tower.</figcaption></figure>

The Fouta Djallon highlands are unanimously recognized as the ‘water tower’ of West Africa as they constitute the source of many transboundary rivers and their tributaries in the region. Table 1 gives the dependence rate of some countries on the FD rivers. In this paper, the acronym ‘FD’ refers to the Fouta Djallon and the broader Guinean Dorsale (Figure 2) that runs southeastwards from the Fouta Djallon mountains, in Guinea, to the Nimba mountains, on the border with Côte d’Ivoire. By contrast, the acronym ‘FDss’ refers to the Fouta Djallon stricto sensu, which only constitutes the northwestern side of the Guinean Dorsale.

<table><caption>**Table 1.** Dependence rate (%) from water coming of the Fouta Djallon highlands of the concerned West African countries [5].</caption><thead><tr><td>Country</td><td>% of Surface Water Resources Originating in the Fouta Djallon</td><td>Concerned River Basins</td></tr></thead><tbody><tr><td>Mauritania</td><td>96</td><td>Senegal</td></tr><tr><td>Gambia</td><td>63</td><td>Gambia</td></tr><tr><td>Senegal</td><td>60</td><td>Senegal, Gambia, Anyamba/Geba</td></tr><tr><td>Guinea-Bissau</td><td>48</td><td>Corubal/Koliba, Anyamba/Geba</td></tr><tr><td>Mali</td><td>55</td><td>Senegal, Niger</td></tr><tr><td>Niger</td><td>70</td><td>Niger</td></tr></tbody></table>
<figure><img src="image_47.png" /><figcaption>**Figure 2.** The main geological parts of the Fouta Djallon and Guinean Dorsale ensemble.</figcaption></figure>

The range gathers the sources of the three main West Sudano-Sahelian rivers: the Gambia River, the sources of the Bafing-Falémé complex which form the Senegal River after the confluence of the Bafing and the Bakoye in Mali, and those of the Niger River and its Guinean tributary on the left bank, the Tinkisso [2]. The range also feeds the two main rivers of the Atlantic side of Guinea, the Kouba/Koumba and the Tominé, which become the Koliba and then the Rio Corubal in Guinea Bissau, as well as the Konkouré complex, and the main rivers crossing Sierra Leone and parts of Liberia.

## 1.1. What Does Not Change

The geological context of the Fouta Djallon differs from a natural water tower. The entire range can be divided into a Precambrian and mostly granitic basement in the south east and a sedimentary ensemble in the north west. The latter is mostly composed of sandstone with doleritic sills and intercalations of argillites (Figure 2). The separation follows a virtual Conakry–Bamako line [6]. In the FDss, the granite is widely covered by hundreds of meters of sandstone, dolomites, and argillites, with frequent sills of dolerites, none of which constitutes a good natural reservoir. Mamedov et al. (2010) [7] estimated the sandstone-quartzite-argillite thickness as ranging from 1500 to 4000 m. The FDss is itself divided into three geological contexts: a small granitic area southward, the western sandstone area (soft sandstone and shale with some injections of dolerites), and the northwestern doleritic FD, which is also mostly composed of sandstone and shale, but which features a large number of outcropping sills and dolerite veins. Boulvert [6] noticed that these thick sandstone plateaus are separated by the incisions (of more than 700 m) of the so-called 'Atlantic southern rivers' (i.e., the Kokoulo, the Kakrima, and the Koumba), highlighted by spectacular sandstone cornices, gorges, and waterfalls. The Guinean Dorsale is mostly an Archean granitic basement with few sills of dolerites and few outcrops of gneiss and quartzite [2]. All the Upper Niger basin and a part of the Konkouré River basin constitute a granitic basement, locally covered by sandstone and the typical 'lateritic pan'. The common characteristics of the whole mountain area of Guinea is the great extension of the 'bowe', i.e., lateritic plateaus composed of ferruginous and overall aluminous crusts [8].

However, no clear difference can be noticed in their hydrological behavior, most probably due to their similarly high levels of imperviousness. Infiltration exists at the FD scale (faults, cracks,
and fractures, very few permeable rocks, local shale, and green rocks), but not at the sample scale. Most of the rocks are either completely impervious (granite, dolerites, quartzite) or mainly impervious (sandstone), except their alterites. In the entire FD, the rocks have a very low water storage capacity [2]. However, Boulvert [8] highlighted the extension of "paleo-crypto karsts" in the Bowe area of Guinea, which covers most of the Fouta Djallon. Therefore, the role of possible occulted streamflows remains to be analyzed.

The main characteristics of the FDss, which are common to the entire FD, can be summarized by the following [9]:

- the total absence of a generalized water table and of continuous groundwater reservoir [9]. There are only small and localized aquifers, well known by the local people [10].

- the very low mineral concentration of water (lower than 50 mg L⁻¹). The average of dissolved load in total sampling observed by Orange (1990a) [9] is only 35 mg L⁻¹, due to the sandstone and/or granite predominance and the importance of stable lateritic crusts [9]. This relatively old finding remains overall coherent and in line with a series of small sampling tests made in the FDss by the authors in March and May 2019. In 28 samples taken at an elevation of less than 800 m, the mineral concentration of water was on average 55 mg·L⁻¹ (standard deviation = 74, with values ranging from 6.5 to 414), and in 14 samples taken at an elevation higher than 800 m, the mineral concentration was on average 24.5 mg·L⁻¹ (standard deviation = 14.8, values from 6.5 to 53).

- the stationary geochemical regime: the alteration profile deepening (4.6 mm·yr⁻⁴) is slightly higher than soil erosion (4.2 mm·yr⁻⁴) [9,10].

The FD's overall impermeability strongly impacts surface water behavior. All of the geological context display the same functioning: very weak discharges during the dry season (from December to June), because they are not supported by groundwater. In contrast, at the occurrence of the first monsoon rainfall events (in June or July), discharges increase suddenly to peak in September, with the rainfall peak [2].

The impervious geology of most of the FD mountains [9,11] leads to high runoff and recession coefficients in the main rivers, due to a very low soil water-holding capacity on the one hand and due to the above-mentioned absence of a generalized water tables on the other.

The monsoon lasts five to six months a year in the core part of the highlands, and annual precipitation is very low during the dry season. The low yet permanent streamflows of the main rivers are then fed by the emptying of the eight thousand small wetlands (commonly settled on alterites), of the thin aquifers, and of the free water of the lowlands alongside rivers and streams of the whole Fouta Djallon and Guinean Dorsale.

## 1.2. What Could Change

D'Aubreville (1947) [12] and then Maignien (1966) [13] were impressed by the great extension of "ferrallitic pans in Guinea". The former thought that the formation of these "bowe" (sing. "bowal"), or "bowalisation" (as the process), is still ongoing, yet most of them have proven to be fossils [8]. These are not evidence of any land degradation processes and results.

During the 1940s and 50s, geographers [14–16] such as Richard-Molard (1943), Pouquet (1956, 1960), agronomists [17] like Sudres (1947), or foresters [12] (e.g., D'Aubreville, 1947) were studying the impact of deforestation and erosion before the 70s [8].

The most densely populated area is the summit of the FDss (around Labé and Mali town) and its historical core. This location is very rainy, but it allows very poor water storage. People encounter severe problems in drinking water supply and even greater challenges to obtain water to irrigate a second or third yearly crop and produce enough food for the entire population.

As observed in the Sahel, the decrease in discharge during the prolonged rainfall deficit (1968–1993) of the Sahel's so-called "Great Drought" is approximately twice that of the rainfall [18]. During the
drought, the authors highlighted that the severe rainfall deficit provoked for the first time in the Senegal
River an exponential increase of the recession (or depletion) coefficient, which was stable for 70 years
(a = 0.0186 d⁻¹) and which strongly increased to reach 0.038 in 1984 [19]. This long "Great Drought"
(Figure 3) led to a lasting degradation of the hydrological system. In Guinean areas, the decrease in
discharges lasted 10 to 15 years. The increase in recession coefficients also lasted for years beyond the
end of the drought. In 1997, Bricquet et al. [20] showed that despite the observed rainfall recovery,
discharges did not begin to increase again at that time. This was due to the necessary refilling of the
various natural storages such as surface and sub-surface reservoirs, and streamflows recovered over
more than a decade [21]. Figure 4 [1] provides the spatial distribution of rainfall.

<figure><img src="image_48.png" /><figcaption>**Figure 3.** Rainfall evolution in the Fouta Djallon mountains since 1923; (a) Senegal and Niger upper basins. (b) Konkouré, Corubal, and Gambia rivers basins. Dots represent the annual total rainfall (ATR), in mm.</figcaption></figure>
<figure><img src="image_49.png" /><figcaption>**Figure 4.** Rainfall spatial distribution (1950–2000) [1] and location of rain gauges (white triangles) used to unskew CHIRPS data.</figcaption></figure>

However, Gomis observed in 2000 [22] that the rainfall recovery led to a decrease in the recession coefficient of the Gambia River. In the Bafing basin (Upper Senegal basin), this slow and progressive discharge recovery was also characterized by a decrease in the recession coefficient, from 0.025 d⁻¹ in 1972–1994 to 0.022 d⁻¹ during the 1995–2005 period at the Bafing Makana stream gauge station [23].

In some basins, dams have been built or are being planned, and they will strongly modify the annual regime and allow hydroelectric power generation, such as in the Konkouré River basin [24].

This paper aims to study the recent evolution of the FD as a natural water tower through the analysis of runoff and depletion coefficients, along other factors such as rainfall and land cover dynamics.

## 2. Material and Methods

This study of the permanency of the Guinean water tower is based on the analysis of the hydrological behavior evolution of 16 basins and sub-basins across these highlands. Basins and their downstream gauging stations are indicated in Figure 5.
<figure><img src="image_50.png" /><figcaption>Figure 5. Location of the 16 basins and sub-basins considered in this study.</figcaption></figure>

Excel file in supplementary materials gives the list of the 16 considered basins and the availability of their hydrological data, where vertical black lines separate the three climatic periods: 1947–1969, 1970–1993 and 1994–2019.

Depletion coefficient (or recession coefficient). The recession coefficient is used after the Maillet law (1905) and as cited by Tallaksen (1995) [25] as follows:

$$Q_t = Q_0 e^{-\alpha t} \qquad (1)$$

$Q_t$ = discharge at t time; $Q_0$ = initial discharge at the beginning of the receding waters; $\alpha$ = recession index (= DP). The latter is the inverse of the elapsed time in d⁻¹ and is expressed as

$$\alpha = -\ln(Q_t/Q_0)/t \qquad (2)$$

The depletion coefficient (DC) used here are the average value of the DC for three durations: one month, two months, and four months (three on the smaller basins). The determined value of the DC is the average value of the DC for each duration around the maximal value of the DC for this duration. The considered period is the one when, on average for the period, the DC is the highest; thus, it is the same for all the years of the period, in the ensuing “Results” part. In a “Supplementary Material” part, we considered, for each station, and for each year, and for three durations (1, 2, and 4 months), the period when the DC reaches its highest value.

### 3. Discharges and Runoff Coefficients

Daily discharge values are provided by the basin agencies (ABN/NBA Niger Basin Authority, for the Niger basin; OMVS—Organisation pour la Mise en Valeur du bassin du fleuve Sénégal- for the Senegal basin; OMVG—Organisation pour la Mise en Valeur du bassin du fleuve Gambie- for the Gambia and Corubal basins) and the National Hydraulic Service of Guinea for the Konkouré basin.
Rainfall data of 23 stations are provided by SIEREM (Système d'Informations Environnementales sur les Ressources en Eau et leur Modélisation) [26,27], AMMA Catch (AMMA = African Monsoon Multidisciplinary Analysis) (http://www.amma-catch.org/), ANACIM (Agence Nationale de l'Aviation Civile et de la Météorologie) (Senegalese data), and from the CHIRPS (Climate Hazards Group InfraRed Precipitation with Station data) data system [28] for recent years (since 1980, 1994, or 2000 according to the basins). The CHIRPS products used here have a spatial resolution of 0.25° and a daily temporal resolution. The covered area is from 50° N to 50° S. They combine satellite data at 0.25° to soil rainfall data [28]. In order to improve the data accuracy, the CHIRPS data were unskewed again using observed rainfall data for a few stations (see Figure 4): Labé, Koundara, Mamou, Dalaba, Faranah, Tougué, Siguiri, Mali, Gaoual, Conakry, Kindia, Dabola, Kissidougou, Kankan, and Kerouané in Guinea; Kedougou, Tambacounda, and Kolda in Senegal; Kayes, Kenieba, Bafing Makana, Bamako, and Falea in Mali.

# 4. Land Use/Land Cover Changes (LULCC) are Determined at Three Scales

## 4.1. At the Regional Scale

At the regional scale, the tree cover variations between 2000 and 2019 over the 16 basins were extracted from the Global Forest Change (GFC2019) datasets produced by the University of Maryland, based on Landsat 4, 5, 7, and 8 imagery. These provide 30 m spatial resolution data of tree gains and tree losses compared to land use cover of the year 2000. Cloud-free observations for each pixel, based on quality assessment layers, are used to create composite imagery in bands 3, 4, 5, and 7 during the growing season. Image interpretation relies on a decision tree supervised learning approach and training data, including very high resolution Quickbird imagery as described in Hansen et al. (2013) [29]. Google Earth Engine© (Google, Silicon Valley, CA, USA) geoprocessing capabilities [30] were used here to reduce the time to extract and combine pixel data for the 16 basins.

## 4.2. At the Basin Scale

At the basin scale, a special analysis was conducted in the upper Bafing basin. At the Upper Bafing basin scale, changes in vegetation were searched using a time series analysis (TSA) of the normalized difference vegetation index (NDVI) based on a MODIS (Moderate-Resolution Imaging Spectroradiometer; EOS Earth Observing System, NASA, Washington, D.C., USA) (MOD13 dataset (250 m resolution; 16 days temporal resolution; from February 2000 to February 2020). Mann-Kendall's tau test is an efficient trend detection method [31] and has been used for this data type for other regions of western Africa [32–34]. Mann-Kendall's tau test can be submitted to null hypothesis (*P* value) testing as the data come from a population with independent realizations and are identically distributed. A significance at 0.001 has been retained for this paper meaning that a sample as extreme as that observed would occur 1 time in 1000 if the null hypothesis is true.

## 4.3. At the Local Scale

At the local scale, some small wetlands of river sources are compared in the central area of the FD. Here, we compare the current status (via Google Earth© images; Google, Silicon Valley, CA, USA) with aerial imagery acquired by the IGN (Institut National de l'Information Géographique et Forestière, Saint Mandé, France) 1953 mission over West Africa.

# 5. Results

## 5.1. Changes of Depletion Coefficients

Figure 6a–e and Table 2A–D gather the results about changes in the DC (depletion coefficient). In these figures, the DC is calculated for a period of four months (and three months for the two Sokotoro and Pont de Linsan smaller basins) for each year on the same date, i.e., the average date of the
maximum depletion of the discharge (in the Supplementary Materials, Figures S1a,b–S16a,b illustrate the evolution of the DC at the highest yearly depletion rate of each year).

<figure><img src="image_51.png" /><figcaption>Figure 6. *Cont.*</figcaption></figure>
<figure><img src="image_52.png" /><figcaption>**Figure 6.** Dynamic of depletion coefficients for the period 1947–2018, by basin; (a) Konkouré River basin; (b) Koliba/Tominé/Corubal River basin; (c) Gambia River basin; (d) Senegal River basin; (e) Niger River basin.</figcaption></figure>
<table><caption>Table 2. Depletion coefficient by basins, by periods, and for different durations.</caption><thead><tr><th rowspan="2">A</th><th rowspan="2">Decreasing DC</th><th rowspan="2"></th><th colspan="3">4 Months Depletion Coef</th><th colspan="3">2 Months Depletion Coef</th><th colspan="3">1 Month Depletion Coef</th><th colspan="3">Yearly</th></tr><tr><th>Mean 1947-70</th><th>Mean 71-93</th><th>Mean 94-2019</th><th>Mean 1947-70</th><th>Mean 71-93</th><th>Mean 94-2019</th><th>Mean 1947-70</th><th>Mean 71-93</th><th>Mean 94-2019</th><th>Mean 1947-70</th><th>Mean 71-93</th><th>Mean 94-2019</th><th>Deforestation %</th></tr></thead><tbody><tr><td rowspan="4">Senegal</td><td rowspan="2">Gourbassi</td><td>Faleme coef</td><td>0.0339</td><td>0.0460</td><td>0.0412</td><td>0.0433</td><td>0.0533</td><td>0.0465</td><td>0.0550</td><td>0.0622</td><td>0.0554</td><td>0.0554</td><td>0.0554</td><td>0.0554</td></tr><tr><td>day</td><td>28-November</td><td>25-November</td><td>1-December</td><td>27-October</td><td>3-November</td><td>3-Nov</td><td>18-October</td><td>25-October</td><td>18-October</td><td>18-October</td><td>18-October</td><td>18-October</td></tr><tr><td rowspan="2">Senegal</td><td>Daka</td><td>0.0238</td><td>0.0243</td><td>0.0237</td><td>0.0298</td><td>0.0321</td><td>0.0317</td><td>0.0352</td><td>0.0395</td><td>0.0392</td><td>0.0392</td><td>0.0392</td><td>0.0392</td></tr><tr><td>Saidou</td><td>day</td><td>29-November</td><td>8-November</td><td>3-December</td><td>7-November</td><td>5-November</td><td>17-November</td><td>5-November</td><td>18-November</td><td>11-November</td><td>11-November</td><td>11-November</td></tr><tr><td rowspan="4">Niger</td><td rowspan="2">Dahola</td><td>Tinkissaoef</td><td>0.0264</td><td>0.0228</td><td>0.0252</td><td>0.0381</td><td>0.0366</td><td>0.0351</td><td>0.0485</td><td>0.0456</td><td>0.0424</td><td>0.0424</td><td>0.0424</td><td>0.0424</td></tr><tr><td>day</td><td>3-December</td><td>20-November</td><td>9-December</td><td>4-November</td><td>25-October</td><td>15-November</td><td>27-October</td><td>17-October</td><td>2-November</td><td>2-November</td><td>2-November</td><td>2-November</td></tr><tr><td rowspan="2">B</td><td rowspan="2">Stabilized DC</td><td>4 Months Depletion Coef</td><td>Mean 1947-70</td><td>Mean 71-93</td><td>Mean 94-2019</td><td>Mean 1947-70</td><td>Mean 71-93</td><td>Mean 94-2019</td><td>Mean 1947-70</td><td>Mean 71-93</td><td>Mean 94-2019</td><td>Mean 1947-70</td><td>Mean 71-93</td><td>Mean 94-2019</td><td>Deforestation %</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="4">Gambia</td><td rowspan="2">Missirah</td><td>Kouloumou</td><td>0.0407</td><td>0.0394</td><td>0.0394</td><td>0.0434</td><td>0.0443</td><td>0.0443</td><td>0.0525</td><td>0.0625</td><td>0.0625</td><td>0.0625</td><td>0.0625</td><td>0.0625</td><td>0.0625</td></tr><tr><td>day</td><td>21-November</td><td>3-December</td><td>4-December</td><td>6-November</td><td>15-November</td><td>15-November</td><td>24-October</td><td>8-November</td><td>8-November</td><td>8-November</td><td>8-November</td><td>8-November</td><td>8-November</td></tr><tr><td rowspan="2">Gambia</td><td>Mako</td><td>0.0294</td><td>0.0308</td><td>4-November</td><td>0.0381</td><td>0.0401</td><td>0.0401</td><td>0.0485</td><td>0.0500</td><td>0.0500</td><td>0.0500</td><td>0.0500</td><td>0.0500</td><td>0.0500</td></tr><tr><td>Gambiaconf</td><td>day</td><td>29-November</td><td>4-November</td><td>27-Oct</td><td>5-November</td><td>5-November</td><td>22-October</td><td>22-October</td><td>29-October</td><td>29-October</td><td>29-October</td><td>29-October</td><td>29-October</td></tr><tr><td rowspan="4">Corubal</td><td rowspan="2">Gaoual K</td><td>Koumbaoef</td><td>0.0201</td><td>0.0197</td><td>9-December</td><td>0.0284</td><td>0.0267</td><td>9-November</td><td>0.0388</td><td>0.0388</td><td>0.0371</td><td>0.0371</td><td>0.0371</td><td>0.0371</td><td>0.0371</td></tr><tr><td>Kolibaklay</td><td>30-November</td><td>9-December</td><td>30-October</td><td>30-October</td><td>9-November</td><td>31-October</td><td>31-October</td><td>14-November</td><td>14-November</td><td>14-November</td><td>14-November</td><td>14-November</td><td>14-November</td></tr><tr><td rowspan="2">Senegal</td><td>Sokotor</td><td>0.0269</td><td>0.0269</td><td>16-November</td><td>0.0314</td><td>16-November</td><td>0.0314</td><td>0.0403</td><td>30-Oct</td><td>30-Oct</td><td>30-Oct</td><td>30-Oct</td><td>30-Oct</td><td>30-Oct</td></tr><tr><td>Bafing coef</td><td>27-November</td><td>27-November</td><td>0.0271</td><td>0.0302</td><td>0.0321</td><td>0.0312</td><td>0.0312</td><td>0.0383</td><td>10-December</td><td>10-December</td><td>10-December</td><td>10-December</td><td>10-December</td></tr><tr><td rowspan="4">Niger</td><td rowspan="2">Kouroussa</td><td>Niger coef</td><td>0.0249</td><td>0.0269</td><td>18-December</td><td>28-December</td><td>25-December</td><td>14-December</td><td>1-December</td><td>26-November</td><td>10-December</td><td>10-December</td><td>10-December</td><td>10-December</td><td>10-December</td></tr><tr><td>day</td><td>27-December</td><td>28-December</td><td>28-December</td><td>25-December</td><td>14-December</td><td>1-December</td><td>26-November</td><td>10-December</td><td>10-December</td><td>10-December</td><td>10-December</td><td>10-December</td><td>10-December</td></tr><tr><td rowspan="2">C</td><td rowspan="2">Increasing DC &lt; 20%</td><td>4 Months Depletion Coef</td><td>Mean 1947-70</td><td>Mean 71-93</td><td>Mean 94-2019</td><td>Mean 1947-70</td><td>Mean 71-93</td><td>Mean 94-2019</td><td>Mean 1947-70</td><td>Mean 71-93</td><td>Mean 94-2019</td><td>Mean 1947-70</td><td>Mean 71-93</td><td>Mean 94-2019</td><td>Deforestation %</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="4">Ddef</td><td rowspan="2">Konkouré</td><td>Nianso</td><td>0.0297</td><td>0.0305</td><td>0.0305</td><td>0.0281</td><td>0.0334</td><td>0.0334</td><td>0.0368</td><td>0.0418</td><td>0.0418</td><td>0.0418</td><td>0.0418</td><td>0.0418</td><td>0.0418</td></tr><tr><td>day</td><td>30-November</td><td>12-December</td><td>0.0295</td><td>14-November</td><td>24-November</td><td>14-November</td><td>1-November</td><td>19-November</td><td>19-November</td><td>19-November</td><td>19-November</td><td>19-November</td><td>19-November</td></tr><tr><td rowspan="2">Konkouré</td><td>Kaba</td><td>0.0278</td><td>0.0295</td><td>2-December</td><td>0.0320</td><td>0.0344</td><td>0.0344</td><td>0.0399</td><td>0.0400</td><td>0.0400</td><td>0.0400</td><td>0.0400</td><td>0.0400</td><td>0.0400</td></tr><tr><td>day</td><td>9-December</td><td>2-December</td><td>29-November</td><td>29-November</td><td>17-November</td><td>17-November</td><td>15-November</td><td>15-November</td><td>15-November</td><td>15-November</td><td>15-November</td><td>15-November</td><td>15-November</td></tr><tr><td rowspan="2">Corubal</td><td rowspan="2">Gaoual T</td><td>Tominéconf</td><td>0.0201</td><td>0.0241</td><td>0.0241</td><td>0.0279</td><td>0.0318</td><td>0.0318</td><td>0.0393</td><td>0.0428</td><td>0.0428</td><td>0.0428</td><td>0.0428</td><td>0.0428</td><td>0.0428</td></tr><tr><td></td><td></td><td></td><td></td><td>30-November</td><td>4-December</td><td>4-December</td><td>30-October</td><td>6-November</td><td>6-November</td><td>6-November</td><td>6-November</td><td>6-November</td><td>6-November</td></tr></tbody></table>
<table><caption>Table 2. Cont.</caption><thead><tr><th rowspan="2">D</th><th rowspan="2">Increasing DC &gt; 20%</th><th rowspan="2">4 Months Depletion Coef mean 1947-70</th><th rowspan="2">2 Months Depletion Coef mean 71-93</th><th rowspan="2">1 Month Depletion Coef mean 94-2019</th><th colspan="6">Yearly</th></tr><tr><th>mean 1947-70</th><th>mean 71-93</th><th>mean 94-2019</th><th>mean 1947-70</th><th>mean 71-93</th><th>mean 94-2019</th><th>mean 1947-70</th><th>mean 71-93</th><th>mean 94-2019</th><th>Deforestation %</th></tr></thead><tbody><tr><td>Niger</td><td>Tinkissò</td><td>Tinkissò coef</td><td>0.0269</td><td>0.0333</td><td>0.0365</td><td>0.0342</td><td>0.0374</td><td>0.0390</td><td>0.0459</td><td>0.0465</td><td>0.0505</td><td>0.235</td><td></td></tr><tr><td rowspan="2">Niger</td><td rowspan="2">Baro</td><td>Niandamoef</td><td>18-December</td><td>21-December</td><td>1-Jan</td><td>24-November</td><td>19-November</td><td>2-December</td><td>14-November</td><td>2-November</td><td>19-November</td><td>0.626</td><td></td></tr><tr><td>14-December</td><td>0.0230</td><td>0.0279</td><td>0.0314</td><td>0.0279</td><td>0.025</td><td>0.0383</td><td>0.0356</td><td>0.0393</td><td>0.0443</td><td>0.626</td><td></td></tr><tr><td rowspan="2">Niger</td><td rowspan="2">Kankan</td><td>Milo</td><td>8-December</td><td>21-December</td><td>22-December</td><td>10-December</td><td>5-December</td><td>27-December</td><td>30-November</td><td>25-November</td><td>24-December</td><td>0.634</td><td></td></tr><tr><td>8-December</td><td>0.0227</td><td>0.0256</td><td>0.0321</td><td>0.0285</td><td>0.0310</td><td>0.0419</td><td>0.0347</td><td>0.0397</td><td>0.0517</td><td>0.634</td><td></td></tr><tr><td>Konkouré</td><td>Pt Rie</td><td>Konkouré coef</td><td>0.0132</td><td>0.0212</td><td>15-November</td><td>16-November</td><td>11-November</td><td>14-November</td><td>2-November</td><td>9-November</td><td>11-November</td><td>0.881</td><td></td></tr><tr><td rowspan="3">Konkouré</td><td rowspan="3">Pt de Linsan</td><td>Upper</td><td>7-January</td><td>26-December</td><td>13-December</td><td>0.0153</td><td>0.0252</td><td>0.0268</td><td>0.0191</td><td>0.0281</td><td>0.0359</td><td>0.731</td><td></td></tr><tr><td>Upper coef</td><td>0.0234</td><td>0.0303</td><td>0.0334</td><td>0.0294</td><td>3-December</td><td>29-November</td><td>26-November</td><td>19-November</td><td>11-November</td><td>0.0460</td><td></td></tr><tr><td>day</td><td>7-November</td><td>14-November</td><td>22-November</td><td>22-October</td><td>6-November</td><td>7-November</td><td>23-October</td><td>31-October</td><td>11-November</td><td>3.731</td><td></td></tr></tbody></table>

Basins are classified into four classes, according to the temporal evolution of the depletion coefficient (DC): decreasing, stabilized, increasing until 20% from the first to third period, and increasing above 20%. Deforestation is here one of the tree cover variations between 2000 and 2019 over the 16 basins, using data of Maryland University (the GFC2019), according to the method proposed by [29] Hansen et al. (2013).
Data covering the 1947–2019 timeline was divided into three periods, based on previous segmentation by several authors ([35–38]):

- the first one is the 'hyper humid' 1947–1969 period, when rainfall was significantly higher than the long-term average across the whole of West Africa.

- the second period is the well-known “Great Drought” of West Africa (1970–1993). The rainfall decreased significantly in 1968, but we arbitrarily made the second period begin in 1970 because a significant number of stream gauge stations have data starting from 1970. The drought is considered to last until 1993.

- the last period spans the years 1994–2019. The year 1994 is when annual rainfall again reached again its long-term mean value. During this period, rainfall is approximately of the long-term trend.

Rainfall predictions in West Africa indicate a regional increase in rainfall during the ongoing 21st century due to the activation of monsoons as a result of the increase in the difference between the continental temperature and SST (sea surface T°) [39]. However, locally, in the western part of West Africa, rainfall is expected to decrease significantly after 2030 or 2035 (Senegambia, Mauritania, and western Mali), while further south, an increase in ATR (annual total rainfall) across most of the Guinean territory is expected [40,41].

The dynamics of the depletion coefficients (see Figure 6a–e and Table 2A–D; see also the Supplementary Materials) shows that:

- the evolution through the three periods is not easy to observe for each basin due to the lack of data (some of the data begins in 1970). The only basins where an analysis over the entire period of 1950–2019 is possible are the largest ones (Senegal and Niger Rivers).

From period 1 to period 2:

- an increase in the DC between period 1 and period 2 is clearly observed in the Upper Niger (at Kouroussa, where it begins in 1965, and at Tinkisso, where an increase in the DC appeared after 1970) and in the Senegal basin. This is clear in the Faleme basin (and probably in the upper Bafing at Sokotoro, though data are lacking). This increase is also observed in the Corubal basin (but particularly over the first drought period, 1970–1975). Despite the lack of data, such a trend also appears in the Gambia River basins (Mako and Missira stream gauges, the second one showing much higher variability). However, it is very difficult to observe this increasing trend in the Konkouré basin, and clearly, the Bafing basin at Daka Saidou does not show any such trends across all available observations.

- only the Niandan at Baro and, to a lesser extent, the upper Tinkisso at Dabola (Figure 5e, but with very limited data) showed a decrease in the DC.

From period 2 to period 3:

- except the Bafing at Daka Saidou, the three main ‘Sudano-Sahelian’ rivers (the Gambia at Mako, the Niger at Kouroussa, as well as the Faleme at Gourbassi) show a decrease in the DC between the dry 1970–1993 period and the following 1994–2019 period. The DC of Koliba Koumba River also decreased slowly after the drought.

- for 8 out of the 16 basins, a decrease (3) or a stabilization (5) in the DC is observed after 1994 and the rainfall recovery. For the other eight stations, an increase in the DC is ongoing. The increase over the whole period is lower than 20% in three basins (Table 2C) and higher than 20% in five others (Table 2D).

- some basins (Niandan, Milo, Gambia, Tominé) show a marked increase in the DC at the end of the last period: in all of these basins, depletion was higher during the 2005–2015 period than during the Great Drought period.
- Extreme values can be observed in the last years (e.g., in the Milo basin), and in 1983–1984 in the Faleme. As described by Olivry (1986), this happens when the higher value of depletion coefficient occurs when the river dries up completely. Consequently, the DC becomes in some cases the drying up coefficient.

General observations

- in some basins, the depletion coefficient can be different depending on the method selected or on the context.

- except the increase in the DC during the drought and decrease after the drought in the large basins, there are no general trends to detect, reflecting great variability across the basins.

- overall, contrary to previous observations by Descroix [5] and Descroix et al. [42], which were based on the behavior of the Upper Niger basin tributaries, there is not a general trend of increasing DC, which could demonstrate a general reduction in the soil water-holding capacity. This trend only appears in the Konkouré basin and in the Milo and Niandan basins.

In Table 2, the last column provides the evolution of forest cover by basin based on the Global Forest Change (GFC2019) methodology.

The day when the maximal depletion occurs seems to follow a general behavior, occurring earlier when the depletion coefficient increases (mostly during the drought), later during the first period (yet with numerous gaps in the data), and overall during the last period. This follows an intuitive pattern, with a depletion stage of the river starting later when discharges are abundant and earlier during the dry years. The exception is the case of the Konkouré basin, where the upper basin (Pont de Linsan station) is the one wherein the discharge has decreased the most. The same evolution is shown at the downstream station "Pont de la route de Télimélé", but this station has been influenced since 1999 by the filling and management of the Garafiri dam. This is why its DC decreased strongly after 1999.

5.2. *Evolution of Runoff Coefficients*

Table 3 gathers the basins’ hydrological characteristics and their changes over the three determined periods: discharge (Q), specific discharge (SD), runoff coefficient (RC), as well as the maximal discharge of the year and its day of occurrence.
<table><caption>Table 3. Hydrological characteristics of the 16 classified basins.</caption><thead><tr><th rowspan="2">B</th><th rowspan="2">Class B<br>Station</th><th rowspan="2">Basin</th><th rowspan="2">Area km<sup>2</sup></th><th colspan="2">Discharge m<sup>3</sup>s<sup>-1</sup></th><th colspan="2">Specific Discharge l<sup>s</sup>-1.km<sup>-2</sup></th><th colspan="2">Runoff Coefficient</th><th colspan="2">Max Daily Discharge/+ Day</th></tr><tr><th>47-70</th><th>71-93</th><th>94-2019</th><th>47-70</th><th>71-93</th><th>94-2019</th><th>47-70</th><th>71-93</th><th>94-2019</th></tr></thead><tbody><tr><td rowspan="2">Senegal</td><td>Daka</td><td>Bafing</td><td>15700</td><td>306.6</td><td>182.4</td><td>225.9</td><td>19.5</td><td>11.6</td><td>14.4</td><td>0.390</td><td>0.264</td><td>0.310</td><td>1750/16<br>September</td><td>1074/29<br>August</td><td>71-93</td><td>94-2019</td></tr><tr><td>Saidou</td><td>Up Bafing</td><td>1750</td><td>29.6</td><td>87.1</td><td>14.0</td><td>16.9</td><td>8.6</td><td>13.7</td><td>0.296</td><td>0.197</td><td>0.308</td><td>337/31<br>August</td><td>1365<br>September</td><td></td><td></td></tr><tr><td rowspan="2">Senegal</td><td>Sokotoro</td><td>Tinkissi</td><td>6370</td><td>89.0</td><td>55.0</td><td>14.5</td><td>11.8</td><td>11.2</td><td>11.5</td><td>0.244</td><td>0.226</td><td>0.242</td><td>67.1/4<br>September</td><td>244/19<br>September</td><td>332/20<br>September</td><td></td></tr><tr><td>Niger</td><td>Up Tinkissi</td><td>1260</td><td>14.9</td><td>14.1</td><td>18.1</td><td>15.5</td><td>10.0</td><td>10.9</td><td>0.286</td><td>0.201</td><td>0.212</td><td>1188/29<br>September</td><td>65.1/5<br>September</td><td>64/2/21<br>September</td><td></td></tr><tr><td rowspan="2">Niger</td><td>Dabola</td><td>Tinkissi</td><td>16560</td><td>256.0</td><td>165.0</td><td>181</td><td>21.5</td><td>16.1</td><td>16.1</td><td>0.348</td><td>0.281</td><td>0.284</td><td>1192/21<br>September</td><td>728/15<br>September</td><td>840/25<br>September</td><td></td></tr><tr><td>Niger</td><td>Niger</td><td>12770</td><td>274.2</td><td>205.2</td><td>205.4</td><td>21.5</td><td>16.1</td><td>16.1</td><td>0.348</td><td>0.281</td><td>0.284</td><td>1192/21<br>September</td><td>865/13<br>September</td><td>922/7<br>September</td><td></td></tr><tr><td rowspan="2">Niger</td><td>Baro</td><td>Niandan</td><td>9620</td><td>206.4</td><td>149.4</td><td>144.9</td><td>21.0</td><td>15.5</td><td>15.1</td><td>0.349</td><td>0.284</td><td>0.262</td><td>771/18<br>September</td><td>665/10<br>September</td><td>712/19<br>September</td><td></td></tr><tr><td>Niger</td><td>Karikan</td><td>Milo</td><td>9620</td><td>206.4</td><td>149.4</td><td>144.9</td><td>21.0</td><td>15.5</td><td>15.1</td><td>0.349</td><td>0.284</td><td>0.262</td><td>771/18<br>September</td><td>665/10<br>September</td><td>712/19<br>September</td><td></td></tr><tr><th rowspan="2">A</th><th rowspan="2">Class A<br>Station</th><th rowspan="2">Basin</th><th rowspan="2">Area km<sup>2</sup></th><th colspan="2">Discharge m<sup>3</sup>s<sup>-1</sup></th><th colspan="2">Specific Discharge l<sup>s</sup>-1.km<sup>-2</sup></th><th colspan="2">Runoff Coefficient</th><th colspan="2">Max Daily Discharge/+ Day</th></tr><tr><th>47-70</th><th>71-93</th><th>94-2019</th><th>47-70</th><th>71-93</th><th>94-2019</th><th>47-70</th><th>71-93</th><th>94-2019</th><th>47-70</th><th>71-93</th><th>94-2019</th></tr><tr><td rowspan="2">Senegal</td><td>Gourhassi</td><td>Faleme</td><td>17100</td><td>166.5</td><td>64.8</td><td>100.4</td><td>9.7</td><td>3.8</td><td>5.9</td><td>0.225</td><td>0.102</td><td>0.149</td><td>1322/1<br>January</td><td>746/6<br>September</td><td>977/10<br>September</td><td></td></tr><tr><td>Gambia</td><td>Mako</td><td>10540</td><td>79.3</td><td>116.6</td><td>79.3</td><td>116.6</td><td>7.6</td><td>11.2</td><td>0.184</td><td>0.247</td><td></td><td></td><td>699/6<br>September</td><td>922/10<br>September</td><td></td></tr><tr><td rowspan="2">Gambia</td><td>Missirah</td><td>Koulountou</td><td>6200</td><td>27.3</td><td>32.9</td><td>32.9</td><td>4.4</td><td>5.3</td><td>0.118</td><td>0.131</td><td></td><td></td><td>170/17<br>September</td><td>170/17<br>September</td><td>209/11<br>September</td><td></td></tr><tr><th rowspan="2">C</th><th rowspan="2">Class C<br>Station</th><th rowspan="2">Basin</th><th rowspan="2">Area km<sup>2</sup></th><th colspan="2">Discharge m<sup>3</sup>s<sup>-1</sup></th><th colspan="2">Specific Discharge l<sup>s</sup>-1.km<sup>-2</sup></th><th colspan="2">Runoff Coefficient</th><th colspan="2">Max Daily Discharge/+ Day</th></tr><tr><th>47-70</th><th>71-93</th><th>94-2019</th><th>47-70</th><th>71-93</th><th>94-2019</th><th>47-70</th><th>71-93</th><th>94-2019</th><th>47-70</th><th>71-93</th><th>94-2019</th></tr><tr><td rowspan="2">Konkouré</td><td>Pt de</td><td>Konkouré</td><td>402</td><td>17.5</td><td>13.1</td><td>12.2</td><td>32.5</td><td>30.5</td><td>0.655</td><td>0.528</td><td>0.476</td><td>113/20<br>August</td><td>88/26<br>August</td><td>93/3/6<br>September</td><td></td></tr><tr><td>Linsan</td><td>Kokoulo</td><td>2260</td><td>60.8</td><td>60.8</td><td>62.1</td><td>26.9</td><td>27.5</td><td></td><td>0.486</td><td>0.439</td><td></td><td>533/28<br>August</td><td>441/29<br>August</td><td>441/29<br>August</td><td></td></tr><tr><td rowspan="2">Konkouré</td><td>Niamsa</td><td>Kakrima</td><td>3190</td><td>72.4</td><td>77.6</td><td>72.4</td><td>22.7</td><td>24.3</td><td></td><td>0.407</td><td>0.392</td><td></td><td>468/27<br>August</td><td>496/26<br>August</td><td>496/26<br>August</td><td></td></tr><tr><td>Kabà</td><td>Pt</td><td>10210</td><td>255.0</td><td>262.0</td><td>262.0</td><td>41.9</td><td>25.0</td><td>25.7</td><td>0.604</td><td>0.410</td><td>0.401</td><td>1971/19<br>August</td><td>1409/23<br>August</td><td>1422/5<br>September</td><td></td></tr><tr><td rowspan="2">Konkouré</td><td>Télémé</td><td>Konkouré</td><td>427.4</td><td>255.0</td><td>262.0</td><td>262.0</td><td>41.9</td><td>25.0</td><td>25.7</td><td>0.604</td><td>0.410</td><td>0.401</td><td>1971/19<br>August</td><td>565/30<br>August</td><td>642/13<br>September</td><td></td></tr><tr><td>Corubal</td><td>Gaoual T</td><td>3300</td><td>98.1</td><td>124.2</td><td>29.7</td><td>37.6</td><td>37.6</td><td>0.529</td><td>0.471</td><td>0.529</td><td>0.529</td><td>0.533</td><td>879/6<br>September</td><td>972/8<br>September</td><td></td></tr><tr><td rowspan="2">Corubal</td><td>Gaoual K</td><td>Koumba</td><td>6200</td><td>174.1</td><td>201.5</td><td>28.5</td><td>33.0</td><td>33.0</td><td>0.500</td><td>0.533</td><td>0.533</td><td>0.533</td><td>0.533</td><td>879/6<br>September</td><td>972/8<br>September</td><td></td></tr></tbody></table>

Basins are classified according to the specific discharges (SD) and runoff coefficient (RC). Class A: SD below 12, RC below 0.2. Class B: SD ranging from 12 to 22, RC ranging from 0.2 to 0.4. Class C: SD above 22, RC above 0.4.
The general dynamic follows the rainfall trends with a decrease in streamflow between the first period (1947–1969) and the second (1970–1993), followed by an increase in discharge between the second period and the third (1994–2019). Yet two important exceptions have to be highlighted:

- the Konkouré basin, where runoff did not rise again after the dry period, but in this basin, data is not available after 2002, preventing analysis of the long term and obscuring a possible recovery of streamflows. However, since rainfall increased significantly between period 2 and 3, a decrease in the runoff coefficient is noticed.

- the Upper Niger basin where runoff increased only in the Tinkisso basin (except for its upper part) over the last period.

This could be explained by a longer time being necessary to refill the natural water reservoir.

The maximum daily discharge observed each year shows approximately the same trends (Table 3, last column): the day when this maximum is observed is in most cases later when the discharge increases and earlier when streamflows are low.

### 5.3. Any Explanation from the LULCC?

The analysis of the relationship between, on the one hand, runoff and depletion coefficients and, on the other hand, land use/land cover changes (LULCC) during the period of 2000–2019 is made with the data provided by the 30 m Global Forest Change datasets [29] (Figures 7 and 8). Considering no local field validation has been done and the inconsistencies in some results, the forest evolution shown here must be considered as relative and not as absolute.

<figure><img src="image_53.png" /><figcaption>**Figure 7.** Vegetation (forest density higher than 15%) in 2000 and 2019 and vegetation loss and gain (GFC (Global Forest Change) product [29]), for the 16 basins. The empty brown rectangle is the area covered by the map in Figure 8.</figcaption></figure>
<figure><img src="image_54.png" /><figcaption>**Figure 8.** Vegetation (density higher than 15%) in 2000 and 2019 and vegetation losses and gains (2001-2019) (GFC product, [29]) for the western basins of the FD, where vegetation losses are the highest. The area represented is the one of the empty brown square in Figure 7.</figcaption></figure>

One of the most important trends is the higher deforestation rate measured in the Upper Niger basin (Kouroussa stream gauge) and the Konkouré basin (mostly in the basins of Pont de Linsan and Pont de la Route de Télimélé stations, secondarily Kaba and Nianso stations). These basins are the ones where:

- RC paradoxically decreases or remains low during the last period, which is characterized by rainfall recovery, possibly because natural reservoirs continue to be filled. Can this hypothesis however still be valid after nearly 25 years of rainfall recovery?

- DC increased strongly (downstream Konkouré basin), moderately (Kaba and Nianso stations), or inversely, and remained unchanged (Upper Niger at Kouroussa station).

- in Table 3, the basins where DC increased more than 20% are clearly not the ones where the deforestation rate was the highest (except the small Pt de Linsan–Upper Konkouré River basin).

# 6. Analysis of LULCC

## 6.1. LULCC at the Regional Scale

GFC2019 rasters reveal that the whole south area of the FD has a high value of vegetation cover; only the Upper Bafing basin is extensively deforested, corresponding to the most densely populated area: this is the extension of the 18th and 19th century Peul/Fulani confederation of small theocratic
states which dominated West Africa. Deforestation in recent decades on the western side of the FD is not necessarily due to the extension of this high density of the upper plateaus area. It is more likely due to wood-charcoal extraction/forest exploitation. A very low number of blue pixels appear in Figures 7 and 8, reflecting the low vegetation gain detected by GFC2019, which is partly due to methodological reasons of vegetation density and vegetation height. The Bebele basin (Tene River at Bebele stream gauge station), which is presented in Figures 7 and 8, is not considered in this analysis, due to its very limited hydrological data set.

The northern areas on Figure 7 correspond to the driest Sudanian zone, covered with bushes, savannah, and low forest cover. Forest vegetation density falls below 15% here, and this limited vegetation is excluded by the vegetation cover threshold employed here.

## 6.2. LULCC at the Upper Bafing Basin Scale

Figure 9 is a map of the Mann-Kendall tau after applying, for each pixel, the P value test at 0.001 and classifying in three categories: significant positive trend, stable or not significant trend, and significant negative trend.

<figure><img src="image_55.png" /><figcaption>**Figure 9.** Evolution of vegetation cover (NDVI—normalized difference vegetation index) in the Upper Bafing Basin from February 2000 to February 2020. The area represented is the one of the empty blue rectangle in Figure 7.</figcaption></figure>

Eighty-six percent of the zone covered by the map (Figure 9) has shown stable values of NDVI from 2000 to 2020 or non-significant trends according to the P value test at 0.001, and 11% of the zone has shown a significant positive trend (greening; biomass increase) for the same period. Only 3% has shown a significant negative trend (browning; biomass decrease).

Unlike the whole Fouta Djalon region, the Upper Bafing area is not concerned by important changes in its vegetation over the last 20 years, and the main change is a positive one for vegetation cover, even if it is hard to say without field validation if it corresponds to higher vegetation cover, higher biomass (either for spontaneous vegetation or for crops), or simple betterment of vegetative activities (i.e., phenology with longer green season).

The basins are not very different regarding this analysis. Stability ranges from 84.1% in the Bafing basin to 90.5% in the Niger basin (Figure 10). Negative trends range from 0.4% in the Gambia basin to 4.4% in the Tinkisso basin. Positive trends range from 6.3% to 14.6% in the Niger and Bafing basins, respectively.
<figure><img src="image_56.png" /><figcaption>**Figure 10.** Part of the vegetation trends on the four main basins in the Bafing River basin.</figcaption></figure>

## 6.3. LULCC at the Sources Scale

The comparison between the 1953 aerial picture (IGN—see Section 4.3 aerial mission in Guinea) and the Google Earth satellite scene (2019) does not demonstrate any apparent deforestations (Figure 11). The quality of both products being very different, it is difficult to do a classical numerical comparison, but the current situation (2019, on the right) seems a priori to be at least as forested as the one in 1953 (Figure 11). This area is named here the “three sources area”: effectively, the two creeks (brooks) flowing through the north are the main sources of the Kouba (Koumba) River and the main upper reach of the Koliba River, which becomes the Rio Corubal in Guinea Bissau. The one flowing eastward is the Dima, the upper Gambia River. ‘Ore Dima’ means “the sources of the Dima River” in Pulaar language (Ore is the head). The other two rivers flowing southwards and the one flowing through the west are the main sources of the Kakrima River, the main tributary of the Konkouré River. This area is mostly composed by a doleritic bedrock.

<figure><img src="image_57.png" /><figcaption>**Figure 11.** The so-called “three sources” area (15 km northward from Labé, the main city of the Fouta Djallon) in 1953 (IGN, **left**) and in 2019 (Google Earth©, **right**). No apparent deforestation. The extension of the picture is 5 km from the west to the east by 4.5 km from the north to the south. The area represented is the one of the red square in Figure 7.</figcaption></figure>

Changes in vegetation cover evolution over time appear different at these three spatial scales but are difficult to compare considering the different analytical methodologies which are here applied
across different areas, periods, and resolutions. However, clearly, some areas are being deforested, though results do not conclude that there is an overall deforestation, contrary to what is commonly described in the literature. Our three studies notably disagree in identifying negative trends in vegetation cover.

### 6.4. Is There Any Synthesis Provided by Data Analysis?

A PCA (Principal Component Analysis) is realized in order to provide clearer information about the variables and observations (basins) sampled in the FD. Table 4 details the variables considered in this analysis. The observations (sampling) are the 16 basins studied.

<table><caption>**Table 4.** List of the variables considered in the PCA analysis.</caption><thead><tr><th>Q</th><th>Mean Yearly Discharge in m³/s</th></tr></thead><tbody><tr><td>QS</td><td>specific discharge (l·s<sup>-1</sup>·km<sup>-2</sup>)</td></tr><tr><td>QDX</td><td>maximum daily discharge</td></tr><tr><td>DC4</td><td>depletion coef in 4 months</td></tr><tr><td>ED41</td><td>DC4 evolution period 1 to per 2</td></tr><tr><td>ED42</td><td>DC4 evolution period 2 to per 3</td></tr><tr><td>DC1</td><td>depletion coef in 1 month</td></tr><tr><td>ED11</td><td>DC1 evolution period 1 to per 2</td></tr><tr><td>ED12</td><td>DC1 evolution period 2 to per 3</td></tr><tr><td>RAD</td><td>rainfall depth</td></tr><tr><td>RUD</td><td>runoff depth</td></tr><tr><td>RC</td><td>runoff coefficient</td></tr><tr><td>ER1</td><td>RC evolution period 1 to per 2</td></tr><tr><td>ER2</td><td>RC evolution period 2 to per 3</td></tr><tr><td>A</td><td>area</td></tr><tr><td>DEFO</td><td>% deforestation 2000–2020</td></tr></tbody></table>

No geological data are accounted for in this study, because it appears that there are no significant differences in the hydrological functioning of the different outcrops in the FD. In both Figures 12 and 13:

- axis 1 (component 1) explains 46.4% of the total variance; it is mostly provided by QS, RAD, RUD, and RC (The definition of all abbreviations is given in Table 4).

- axis 2 (component 2) explains 23.8% of the variance, provided by Q, QDX, and A.

- axis 3 (component 3) explains 11.5% of the variance, using DC4, ED41, ED11, and ER1.

- axis 4 (component 4) explains 7.6% of the variance, provided by ED42 and ED12.

<figure><img src="image_58.png" /><figcaption>**Figure 12.** Variables spaces of the PCA (Principal Component Analysis) analysis; (a) components 1 and 2; (b) components 3 and 4.</figcaption></figure>
<figure><img src="image_59.png" /></figure>

<figure><img src="image_60.png" /><figcaption>**Figure 13.** Observations spaces of the PCA analysis; (a): components 1 and 2; (b) components 3 and 4.</figcaption></figure>

Figure 12a shows some redundancy of, on the one hand, hydrological variables (RAD, RUD, QS, RC) and, on the other, the depletion coefficient, separated from their evolution in two different clusters. They are, by contrast, gathered according to axis 3 and axis 4 (Figure 12b). The LULCC is considered here through the regional mapping. Although its results are different (and locally contradictory) with the ones provided by the other LULCC analyses, we considered these to be suitable as relative values rather than as absolute values. The LULCC is also then determined with the same data and
methods for the 16 basins. Deforestation as shown in Figure 7 is correlated with the increase in runoff between period 1 and 2, but negatively correlated with the increase in runoff between period 2 and 3. Relative deforestation is clearly correlated with the increase in the DC (depletion coefficient) in four months during the second period, the one when rainfall recovery should have allowed a significant regreening. However, it is not so clearly correlated with the DC in one month (Figure 12a); yet following components 3 and 4 (Figure 12b), they seem more clearly correlated.

Figure 13 shows firstly a great dispersion of samples (basins) in all the statistical spaces, as evidence of their great geographical and hydrological diversity. Figure 13a seems logically dominated by runoff, and then, the basins of the western slope of the FDss and the ones coming off the Guinean Dorsale are clearly opposed to the ones of the northern part, where streamflow is weaker. No general trend appears in Figure 13b.

# 7. Discussion

The following points must be highlighted:

## 7.1. The DC

A strong interannual increase in the DC (see also [5,42]) is observed for the two main tributaries of the Upper Niger River from 1995. However, this is not observed in the Upper Niger (Kouroussa) nor in the Tinkisso River (Figure 6e). Otherwise, the DC, according to the three durations tested here, can be similar (e.g., in the Niandan River) or very different (e.g., in the Milo River).

A great deal was devoted to the DC during the 1990s and the 2000s by Mahé et al. [43,44] (1997 and 2001) and Olivry [18,19,21] (1987, 1993, 2002), as well as Bricquet et al. [20] (1997), Bamba et al. [45] (1996), and Sangaré et al. [46] (Figure 14). The latter showed that the DC increased from the beginning of the 1970s in all the basins, even if more slowly for the Sankarani River (one of the other main tributaries of the Upper Niger River, not investigated here). The exception was then the Milo River, the one whose DC has further increased in recent years. The DC of the Milo River only started its increase at the beginning of the 1980s [39]. A possible role of the basins’ water-holding capacity remains to be investigated because the geological context is probably more complex than the schematic regionalization in Figure 2. The Upper Niger River basin is mostly a granitic basement one, however, Brunet Moret et al. (1986) [47] showed that there are sometimes sandstone outcrops in the lower basins of the Milo River. This could have served as a better water reservoir than granite and thus reduced the impact of the first drought spell in the 1970s. Another explanation can be provided by the rainfall’s variation coefficient, which is lower in the Upper Milo basin than in most places in West Africa, as was noticed by Mahé et al. [44], most probably due to the specific relief orientation in the area. The Upper Milo basin runs through elevated areas both in the South and in the East, which seem to play a role in the reduction of the interannual variability, as is also observed in Central Ghana [44]. As can be seen (Figure 14) the DC of the Sankarani river also showed a reduced variability compared to that of other Upper Niger River tributaries, during the Great Drought. In the Senegal River basin, Bodian et al. (2020) [48] observed that the runoff recovery after the Great Drought lasted for five or ten years and is more significant for yearly and maximum discharge and less for low waters, which remain very weak. The same authors concluded that the high value of the DC for the Faleme River at Gourbassi (0.62) at the beginning of the 1980s (Figure 6d) can be attributed to a quasi-total drying up of the streamflow. A similar observation can be made for the Milo River during the 2010s (Figure 6e), but the reason remains to be explained since neither rainfall evolution nor the LULCC are satisfactory explaining factors for this behavior.
Niger River basin, tributaries, 1950-1990

<figure><img src="image_61.png" /><figcaption>**Figure 14.** Evolution of DC in some Sudano-Guinean tributaries of the Niger River before and during the drought [41].</figcaption></figure>

## 7.2. LULCC (Land Use/Land Cover Change)

As it has been shown above, vegetation cover in West Africa does not decline systematically though the population keeps on increasing at a rapid pace [5]. The Fouta Djallon is one of the best examples of this trend. A negative and Afro-pessimistic belief often links African deforestation with demographic growth. This stereotype anticipated catastrophic scenarios based on extensive peasant practices, which, against the background of strong demographic growth, should threaten the integrity of the Fouta Djallon ecosystems and therefore endanger the whole 'West African water tower'. In fact, local research has opened the debate on the reality and validity of such claims [49]. Brandt et al. (2017) [50] showed localized reductions in woodlands across West Africa but that farmland in Sahelian areas also promotes woody cover, contradicting simplistic ideas of a high negative correlation between population density and woody cover.

This strengthens the Boserupian point of view [51] (Boserup, 1965), increasingly accepted in the African scientific literature, as in Kenya [52], Niger [53], Northern Côte d'Ivoire [54], Northeastern Nigeria [55] or the whole Sudano-Sahelian strip [5,42]. This perspective is opposed to the catastrophist, Neo-Malthusian, and Afro-pessimistic general consensus of the second half of the 20th century.

Finally, it is worth noting that the 'tapades' traditional intensive gardening and tree cropping practices of the Fouta Djallon, which have played the role of man-made sponges though they have mostly been managed by women, are largely being abandoned. These practices should be rehabilitated in order to promote a sustainable agricultural intensification as the most useful way to infiltrate and store water in a catchment [5].

## 7.3. Some Hydrological Consequences

Except during the drought, the increase in the DC has only been observed in a few basins in recent years (the Konkouré and southern tributaries of the Upper Niger). The partial recovery in the RC (runoff coefficient), as well as the observed increase—at the regional scale—of the occurrence of extreme rainfall events [36,42,56,57], has consequences for flooding frequency and intensity.
This increases streamflow irregularities, accentuating low waters as well as floods. This could have impacts on hydraulic and civil engineering equipment and infrastructures. As illustrated in Figures 15 and 16, the bridges are threatened by some of these extreme flooding events.

<figure><img src="image_62.png" /><figcaption>**Figure 15.** Bridge of the Koundara-Youkounkoun road, northern Guinea, on the Koulountou River. Built in the 1950s, it was flooded by the streamflow in August 2018, depositing whole trees (bottom) and wood debris (on the bridge deck).</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>**Figure 16.** Pieces of branches blocked under the bridge deck of the new bridge of Komba on the Komba River (Upper Koliba/Corubal basin) after the flood of 2018. On the small picture: the new (2012) and, in the back, the ancient bridge (1950s).</figcaption></figure>
However, as this paper attempted to demonstrate, it is not necessarily the sign of a decrease in soil water-holding capacity attributable to LULCC. Depletion coefficients of the main rivers providing freshwater resources to the Sahel are decreasing due to the rainfall recovery. West African agro-ecosystems seem to be developing resilience, which could invert the observed trend and constitute rural resiliency models and sustainable development practices.

## 8. Conclusions

The Fouta Djallon seems to have suffered from a reduction in water-holding capacity during the drought (1970s and 1980s). However, after this dry period, in some basins and mostly in those flowing towards Sahelian areas, the expected decrease in the DC due to rainfall recovery is observed (the Niger at Kouroussa, the Falémé at Gourbassi, probably the Gambia and Corubal basin rivers). The Bafing at Daka Saidou did not show any increases in DC during the drought.

In the Upper Niger basin (the Milo and Niandan rivers), the Konkouré, and partially the Koliba/Tominé basins, a significant increase in the DC was noticed during the end of the recent 1994-2019 period. However, in the Upper Niger basin, the DC increased more in the basins where relative deforestation was lower (Milo and Niandan) and decreased slightly in the basin where relative deforestation was stronger (Upper Niger at Kouroussa). Following the PCA analysis, only the increase in the DC in four months between period 2 and period 3 is clearly correlated with relative deforestation. Besides the well-known and well-studied increase in rainfall intensity (increase in extreme events occurrence and increase in rainfall hourly mean intensity), the LULCC could be a discriminant factor of any increase in runoff and depletion coefficients, particularly in the Konkouré basin where deforestation remains an important dynamic. The PCA analysis only discriminates the southwestern basins (the Konkouré and Koliba/Tominé basins).

Contrary to popular belief, there is evidence of successful rural intensification practices in West Africa, which are generally characterized by an increasing vegetation cover and, locally, increasing tree cover. This kind of landscape should logically improve the soil and basins' water-holding capacity and eventually lead to a reduction in the DC (depletion coefficient). In the core of the FD (i.e., the FDss), the 'tapades' intensive system is disappearing, yet no deforestation is observed due to several combined patterns.

Since the sole increase in extreme rainfall events is unlikely to be responsible for the increase in the DC, LULCC seems to be a factor in the reduction in the basin water-holding capacity. However, it is statistically apparent only in the Konkouré and in the Koliba/Tominé basins, where annual rainfall is higher and where the runoff coefficients are significantly higher than in all the other basins, making their susceptibility to changes in runoff conditions more influential than in other basins.

The Fouta Djallon, as the West African water tower, does not appear as threatened. However, there is a stark difference between, on the one hand, the southwestern basins (the Konkouré and Upper Corubal), where LULCC seems to impact the basin water-holding capacity, and, on the other hand, basins such as the Gambia, Senegal, and Niger basins, which have been impacted more by the Great Drought but which are at least partially recovering their natural water storage capacity. This is an optimistic observation for the future of the West African Sahel, whose basins depend the most on the Fouta Djallon for its supply of fresh water resources.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/11/2968/s1. Excel 1: 16 considered basins and the availability of their hydrological data. Figure S1a: Evolution of depletion coefficient, Konkouré River basin. Figure S1b: Evolution of rainfall and runoff coefficient, Konkouré River basin. Figure S2a: Evolution of depletion coefficient, Upper Konkouré River basin. Figure S2b: Evolution of rainfall and runoff coefficient, Konkouré River basin. Figure S3a: Evolution of depletion coefficient, Kokoulo River basin. Figure S3b: Evolution of rainfall and runoff coefficient, Kokoulo River basin. Figure S4a: Evolution of depletion coefficient, Kakrima River basin. Figure S4b: Evolution of rainfall and runoff coefficient, Kakrima River basin. Figure S5a: Evolution of depletion coefficient, Tominé River basin. Figure S5b: Evolution of rainfall and runoff coefficient, Tominé River basin. Figure S6a: Evolution of depletion coefficient, Koliba/Koumba River basin. Figure S6b: Evolution of rainfall and runoff coefficient, Koliba/Koumba River basin. Figure S7a: Evolution of depletion coefficient, Koulountou River basin. Figure S7b: Evolution of rainfall and runoff coefficient, Koulountou River basin.
River basin. Figure S8a: Evolution of depletion coefficient, Gambia River basin. Figure S8b: Evolution of rainfall and runoff coefficient, Gambia River basin. Figure S9a: Evolution of depletion coefficient, Faleme River basin. Figure S9b: Evolution of rainfall and runoff coefficient, Faleme River basin. Figure S10a: Evolution of depletion coefficient, Upper Bafing River basin. Figure S10b: Evolution of rainfall and runoff coefficient, Upper Bafing River basin. Figure S11a: Evolution of depletion coefficient, Bafing River basin. Figure S11b: Evolution of rainfall and runoff coefficient, Bafing River basin. Figure S12a: Evolution of depletion coefficient, Upper Tinkisso River basin. Figure S12b: Evolution of rainfall and runoff coefficient, Upper Tinkisso River basin. Figure S13a: Evolution of depletion coefficient, Tinkisso River basin. Figure S13b: Evolution of rainfall and runoff coefficient, Tinkisso River basin. Figure S14a: Evolution of depletion coefficient, Upper Niger River basin. Figure S14b: Evolution of rainfall and runoff coefficient, Upper Niger River basin. Figure S15a: Evolution of depletion coefficient, Niandan River basin. Figure S15b: Evolution of rainfall and runoff coefficient, Niandan River basin. Figure S16a: Evolution of depletion coefficient, Milo River basin. Figure S16b: Evolution of rainfall and runoff coefficient, Milo River basin.

**Author Contributions:** Conceptualization, L.D., B.F., S.P.M. and A.B.D.; data curation, B.F., S.P.M., A.B.D., S.S. (Safietou Soumaré), F.B.S.D., A.D., K.D. (Kadiatou Diallo), J.A. (Julien Andrieu), I.A., J.-C.B., A.B. (Aliou Barry), A.B. (Ansoumana Bodian), Y.B., N.B., J.-L.C., H.D., G.D., M.D., K.D. (Kourahoye Diallo), A.F., E.G., L.F., L.K., B.N.N., Y.S., S.S. (Saly Sambou) and B.A.T.; formal analysis, B.F., S.P.M., A.B.D., J.A. (Julien Andrieu), A.O., J.A. (Jean Albergel), and J.-C.B.; funding acquisition, J.A. (Julien Andrieu), I.A., J.A. (Jean Albergel), J.-C.B., A.B. (Aliou Barry), A.B. (Ansoumana Bodian), and Y.B.; investigation, S.P.M., A.B.D., S.S. (Safietou Soumaré), B.F., F.B.S.D., A.D., K.D. (Kadiatou Diallo), and J.-P.V.; methodology, A.B.D., S.S. (Safietou Soumaré), J.A., G.M., A.O., D.O., and B.F.; project administration, L.D., J.A. (Julien Andrieu), G.M., J.A. (Jean Albergel), A.O., and D.O.; resources, F.B.S.D., A.D., and K.D. (Kadiatou Diallo); software, B.F., F.B.S.D., and K.D. (Kourahoye Diallo); supervision, J.-P.V., L.A.L., J.-C.O. and L.D.; validation, G.M., J.A. (Julien Andrieu), A.O., D.O., L.A.L., and J.-P.V.; visualization, F.B.S.D., A.D., and K.D. (Kadiatou Diallo); writing—original draft, L.D., B.F., S.P.M., and A.B.D.; writing—review and editing, L.D., B.F., S.P.M., and A.B.D.; G.M., J.A. (Julien Andrieu), A.O., J.-C.B., D.O., L.A.L., and J.-P.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the European Union, via AICS (Agenzia Italiana de Cooperazione per il Sviluppo, grant number ENV/2017/383-744/AP n) 01/2019/WEFE SENEGAL for the OMVS, Senegal River Basin Agency.

**Acknowledgments:** The authors warmly acknowledge the Rivers Basins Agency (NBA, OMVS and OMVG) as well as the Direction of Hydraulics of Guinea for providing the daily discharge data, and the SIEREM and AMMA Catch databases as well as the ANACIM (Senegalese Met Office) for the daily rainfall data. We kindly acknowledge the reviewers who allowed us to propose a deeply improved version of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. L'Hôte, J.-P.; Mahé, G. *Afrique Centrale et de l'Ouest: Précipitations Moyennes Annuelles* (Période 1951-1989); IRD (Orstom): Bondy, France, 1996; 1 map.

2. Aurouet, A.; Ferry, L.; Cougny, G.; ECOWAS. *Atlas de l'eau du Massif du Fouta Djalon-Le Château d'eau de l'Afrique de l'Ouest/Fouta Djallon Highland-Water Atlas*; CEDEAO (ECOWAS), IWA, Anteagroup, The World Bank: Washington, DC, USA, 2017; p. 114. Available online: http://documents.banquemondiale.org/curated/fr/844511520436091243/pdf/123986-WP-P150210-FRENCH-PUBLIC.pdf (accessed on 11 April 2020).

3. Gebrehiwot, S.G.; Gardenas, A.I.; Bewket, W.; Seibert, J.; Ilstedt, U.; Bishop, K. The long-term hydrology of East Africa's water tower: Statistical change detection in the watersheds of the Abbay Basin. *Reg. Environ. Chang.* **2013**. [CrossRef]

4. Blanc, P. De l'Égypte à l'Éthiopie, quand la puissance se déplace en Afrique Nilotique. In *Confluences Méditerranée*; L'Harmattan: Paris, France, 2014; ISSN 1148-2664; ISBN 9782343045627. Available online: https://www.cairn.info/revue-confluences-mediterranee-2014-3-page-123.htm (accessed on 21 June 2020).

5. Descroix, L. *Processus et Enjeux d'eau en Afrique de l'Ouest Sahélo-Soudanienne*; Editions Des Archives Contemporaines: Paris, France, 2018; 320p, ISBN 9782813003140. [CrossRef]

6. Boulvert, Y. Approche Synthétique des Aplanissements Cuirassés de Centrafrique et de Guinée (Conakry). 2005. Available online: http://horizon/documentation.ird.fr/exl-doc/pleins_textes/divers19-10/010038259.pdf (accessed on 18 April 2020).

7. Mamedov, V.I.; Boufeev, Y.V.; Nikitine, Y.A. *Geology of the Republic of Guinea*; Geoprospects and State University Lomonosov of Moscow: Moscow, Russia, 2010; 326p.
8. Boulvert, Y. Repères Historiques sur la Caractérisation du Cuirassement et de la Géomorphologie Guinéenne. Conférence à l'Université de Conakry, Conakry (GIN), 1998/04. Available online: https://numerisud.ird.fr/documents-et-films/publications/Reperes-historiques-sur-la-caractérisation-du-cuirassement-et-de-la-geomorphologie-guinéenne (accessed on 18 April 2020).

9. Orange, D. Hydroclimatologie du Fouta Djallon et Dynamique Actuelle d'un Vieux Paysage Latéritique. Ph.D. Thesis, Université Louis Pasteur de Strasbourg, Strasbourg, France, 1990; 232p. Available online: https://www.persee.fr/doc/sgeol_0302-2684_1992_mon_93_ (accessed on 30 October 2018).

10. Orange, D.; Gac, J.-Y. Reconnaissance géochimique des eaux du Fouta Djalon (Guinée), flux de matières dissoutes et en suspension en Haute-Gambie. *Géodynamique* **1990**, *5*, 35–49.

11. Michel, P. *Les Bassins des Fleuves Sénégal et Gambie, Etude Géomorphologique*; ORSTOM: Bondy, France, 1973; 810p.

12. D'Aubreville, A. Erosion et "bovalisation" en Afrique noire française. *L'Agron. Trop.* **1947**, *2*, 339–357.

13. Maignien, R. Compte-rendu de recherches sur les latérites. In *Colloquium on Laterites, Madagascar* 1964; Reports on Natural Ressources, n°4; UNESCO: Paris, France, 1966; 155p.

14. Richard-Molard, J. Les traits d'ensemble du Fouta-Djalon. *Rev. Géograph. Alp.* **1943**, XXXI, 199-213.
[CrossRef]

15. Pouquet, J. Le plateau de Labé (Guinée française, A.O.F). Remarques sur le caractère dramatique des phénomènes d'érosion des sols et sur les remèdes proposés. *Bull. l'IFAN (Inst. Fr. d'Afr. Noire)* **1956**, *18*, 1-16.

16. Pouquet, J. Quelques types d'évolution du relief en Guinée française: Processus hydrographiques et phénomènes de cuirassement sur les hautes surfaces du Fouta-Dialon (AO.F.). In *Proceedings of the Eighteenth International Congress of Geography, Rio de Janeiro, Brazil, 9-18 August 1956*; pp. 350-361.

17. Sudres, A. La dégradation des sols au Foutah Djalon. *L'Agron. Trop.* **1947**, I, 227-246.

18. Olivry, J.-C. *Synthèse des Connaissances Hydrologiques et Potentiel en Ressources en Eau du Fleuve Niger*; Provisional Report; World Bank, Niger Basin Authority: Niamey, Niger, 2002; 160p.

19. Olivry, J.-C. Les conséquences durables de la sécheresse actuelle sur l'écoulement du fleuve Sénégal et l'hypersalinisation de la Basse-Casamance. In *The Influence of Climate Change and Climatic Variability on the Hydrologic Regime and Water Resources; Proceedings of the Vancouver Symposium*; IAHS Publ. No. 168: Wallingford, UK, 1987; 12p.

20. Bricquet, J.-P.; Bamba, F.; Mahé, G.; Touré, M.; Olivry, J.-C. Evolution récente des ressources en eau de l'Afrique atlantique. *Rev. Sci. l'Eau* **1997**, *3*, 321-337. [CrossRef]

21. Olivry, J.-C.; Bricquet, J.-P.; Mahé, G. Vers un appauvrissement durable des ressources en eau de l'Afrique humide? In *Hydrology of Warm Humid Regions. Proceedings of the 4th Assembly IAHS, Yokohama, Japan, 13-15 July 1993*; Gladwell, J.S., Ed.; Publication IAHS 216: Wallingford, UK, 1993; pp. 67-78.

22. Gomis, D.E.R. *Synthèse Hydrologique Du Fleuve GAMBIE en Amont de Gouloumbou*; Master Memory; UCAD: Dakar, Senegal, 2000; 166p.

23. Bodian, A.; Dacosta, H.; Dezetter, A. Analyse des débits de crues et d'étiages dans le bassin versant du fleuve sénégal en amont du barrage de Manantali. *Clim. Dev.* **2013**, *15*, 46–56.

24. Samoura, K. Contributions Méthodologiques à L'évaluation Environnementale Stratégique De L'exploitation du Potentiel Hydroélectrique des Bassins Côtiers en Milieu Tropical: Cas du Konkouré, en Guinée. Ph.D. Thesis, UQAM, Montreal, QC, Canada, 2011; 304p.

25. Tallaksen, L.M. A review of baseflow recession analysis. *J. Hydrol.* **1995**, *165*, 349-370. [CrossRef]

26. Boyer, J.F.; Dieulin, C.; Rouché, N.; Crès, A.; Servat, E.; Paturel, J.E.; Mahé, G. SIEREM: An environmental information system for water resources. In *Water Resource Variability: Hydrological Impacts. Proceedings of the 5th FRIEND World Conference, Havana, Cuba, November 2016*; IAHS Publ.: La Havana, Cuba, 2006; Volume 308, pp. 19-25.

27. Dieulin, C.; Mahé, G.; Paturel, J.-E.; Ejjiyar, S.; Tramblay, Y.; Rouché, N.; El Mansouri, B. A new 60-year 1940-1999 monthly gridded rainfall data set for Africa. *Water* **2019**, *11*, 387. [CrossRef]

28. Funk, C.; Peterson, P.; Landsfeld, M.; Pedreros, D.; Verdin, J.; Shukla, S.; Husak, G.; Rowland, J.; Harrison, L.; Hoell, A.; et al. The climate hazards infrared precipitation with stations—A new environmental record for monitoring extremes. *Sci. Data* **2015**, 2, 150066. [CrossRef] [PubMed]

29. Hansen, M.C.; Potapov, P.V.; Moore, R.; Hancher, M.; Turubanova, S.A.; Tyukavina, A.; Thau, D.; Stehman, S.V.; Goetz, S.J.; Loveland, T.R.; et al. High-Resolution Global Maps of 21st-Century Forest Cover Change. *Science* **2013**, *342*, 850–853. [CrossRef] [PubMed]
30. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. *Remote Sens. Environ.* **2017**, *202*, 18–27. [CrossRef]

31. Gilbert, R.O. *Statistical Methods for Environmental Pollution Monitoring*; Van Nostrand Reinhold: New York, NY, USA, 1987; ISBN 0-442-23050-8.

32. San Emeterio, J.-L.; Alexandre, F.; Andrieu, J.; Génin, A.; Mering, C. Changements socio-environnementaux et dynamiques des paysages ruraux le long du gradient bioclimatique nord-sud dans le sud-ouest du Niger (régions de Tillabery et de Dosso. *VertigO Rev. Electron. Sci. L'Environ.* **2013**, *13*. Available online: http://journals.openedition.org/vertigo/14456 (accessed on 12 October 2015). [CrossRef]

33. Andrieu, J. Phenological analysis of the savanna-forest transition from 1981 to 2006 form Côte d'Ivoire to Benin with NDVI NOAA time series. *Eur. J. Remote Sens.* **2017**. [CrossRef]

34. Andrieu, J. Vegetation change analysis in Côte d'Ivoire during conflicts using a phenological metric and Kendall correlation of two NDVI time series (Text in French). *Tropicultura* **2018**, *36*, 258–270. Available online: http://www.tropicultura.org/text/v36n2/258.pdf (accessed on 9 November 2019).

35. Nicholson, S.E. On the question of the "recovery" of the rains in the West African Sahel. *J. Arid Environ.* **2005**, *63*, 615–641. [CrossRef]

36. Ali, A.; Lebel, T. The Sahelian standardized rainfall index revisited. *Int. J. Climatol.* **2009**, *29*, 1705–1714. [CrossRef]

37. Nicholson, S.E. The West African Sahel: A review of recent studies on the rainfall regime and its interannual variability. *Int. Sch. Res. Not.* **2013**, *2013*, 453251. [CrossRef]

38. Descroix, L.; Diongue Niang, A.; Panthou, G.; Bodian, A.; Sané, T.; Dacosta, H.; Malam Abdou, M.; Vandervaere, J.-P.; Quantin, G. Evolution Récente de la Mousson en Afrique de l'Ouest à Travers Deux Fenêtres (Sénégambie et Bassin du Niger Moyen). *Climatologie* **2015**, *12*, 25–43.

39. Ardoin-Bardin, S.; Dezetter, A.; Servat, E.; Paturel, J.-E.; Mahé, G.; Niel, H.; Dieulin, C. Using general circulation model outputs to assess impacts of climate change on runoff for large hydrological catchments in West Africa. *Hydrol. Sci. J.* **2009**, *54*, 77–89. [CrossRef]

40. Biasutti, M. Forced Sahel rainfall trends in the CMIP5 archive. *J. Geophys. Res. Atmos.* **2013**, *118*, 1613–1623. [CrossRef]

41. IPCC. *Climate Change Synthesis Report*; IPCC, WMO: Geneva, Switzerland, 2014; 167p.

42. Descroix, L.; Guichard, F.; Grippa, M.; Lambert, L.A.; Panthou, G.; Gal, L.; Dardel, C.; Quantin, G.; Kergoat, L.; Bouaïta, Y.; et al. Evolution of surface hydrology in the Sahelo-Sudanian stripe: An updated synthesis. *Water* **2018**, *10*, 748. [CrossRef]

43. Mahé, G.; Bamba, F.; Diabaté, M.; Diarra, A.; Diarra, M. The reduction of the water resources on upper basins of the Niger river: Hydrological balances and analysis of the depletion curves (1951–1989). Poster proceedings, Sustainability of water resources under increasing uncertainty. In Proceedings of the 5th IAHS Assembly, Rabat, Maroc, 23 April–3 May 1997; pp. 9–12.

44. Mahé, G.; L'Hôte, Y.; Olivry, J.C.; Wotling, G. Trends and discontinuities in regional rainfall of west and central Africa–1951–1989. *Hydrol. Sci. J.* **2001**, *46*, 211–226. [CrossRef]

45. Bamba, F.; Mahé, G.; Bricquet, J.P.; Olivry, J.C. Changements climatiques et variabilité des ressources en eau des bassins du Haut Niger et de la Cuvette Lacustre. In *Réseaux Hydrométriques, Réseaux Télématiques, Réseaux Scientifiques: Nouveaux Visages de l'Hydrologie Régionale en Afrique*; Fritsch, J.M., Paturel, J.E., Servat, E., Eds.; XIIèmes Journées Hydrologiques de l'ORSTOM: Montpellier, France, 1996; 26p.

46. Sangaré, S.; Mahé, G.; Paturel, J.-E.; Bangoura, Y. Bilan hydrologique du fleuve Niger en Guinée de 1950 à 2000. *Sud Sciences et Technologies* **2002**, *9*, 21–33. Available online: https://horizon/documentation.ird.fr/exl-doc/pleins_textes/divers16-04/010034319.pdf (accessed on 18 June 2017).

47. Brunet Moret, Y.; Chaperon, P.; Lamagat, J.-P.; Molinier, M. *Monographie Hydrologique du Fleuve Niger*; Orstom: Paris, France, 1986; 800p.

48. Bodian, A.; Diop, L.; Panthou, G.; Dacosta, H.; Deme, A.; Dezetter, A.; Ndiaye, P.M.; Diouf, I.; Vischel, T. Recent Trend in Hydroclimatic Conditions in the Senegal River Basin. *Water* **2020**, *12*, 436. [CrossRef]

49. André, V.; Pestaña, G. Les visages du Fouta-Djalon. *Les Cah. d'Outre-Mer* **2002**, *217*. Available online: http://com.revues.org/index1038.html (accessed on 11 September 2014). [CrossRef]

50. Brandt, M.; Rasmussen, K.; Peñuelas, J.; Tian, F.; Schurgers, G.; Verger, A.; Mertz, O.; Palmer, J.R.B.; Fensholt, R. Human population growth offsets climate-driven increase in woody vegetation in sub-Saharan Africa. *Nat. Ecol. Evol.* **2017**, *1*, 0081. [CrossRef]
51. Boserup, E. *The Conditions of Agricultural Growth: The Economics of Agrarian Change under Population Pressure*; (Republished 1993: Earthscan Publications: London, UK); Allen and Unwin: London, UK, 1965.

52. Tiffen, M.; Mortimore, M.; Gichuki, F. *More People, Less Erosion: Environmental Recovery in Kenya*; John Wiley & Sons: London, UK, 1994; 311p.

53. Luxereau, A.; Roussel, B. *Changements Ecologiques et Sociaux au Niger*; Etudes Africaines; L'Harmattan Ed: Paris, France, 1997; 239p.

54. Demont, M.; Jouve, P. *Evolution d'Agro-Systèmes Villageois Dans la Region de Korhogo (ord Côte d'Ivoire)*: Boserup vs. Malthus, Opposition ou Complémentarité? *Dynamiques Agraires et Construction Sociale Du Territoire*; Séminaire CNEARC-UTM: Montpellier, France, 2000; pp. 93–108.

55. Mortimore, M.J.; Adams, W.M. Farmer adaptation, change and &crisis' in the Sahel. *Glob. Environ. Chang.* **2001**, *11*, 49–57.

56. Panthou, G.; Vischel, T.; Lebel, T. Recent trends in the regime of extreme rainfall in the Central Sahel. *Int. J. Climatol.* **2014**, *34*, 3998–4006. [CrossRef]

57. Panthou, G.; Lebel, T.; Vischel, T.; Quantin, G.; Sané, Y.; Ba, A.; Ndiaye, O.; Diongue-Niang, A.; Diop Kane, M. Rainfall intensification in tropical semi-arid regions: The Sahelian case. *Environ. Res. Lett.* **2018**, *13*, 064013. [CrossRef]

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

© 2020 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 (http://creativecommons.org/licenses/by/4.0/).
<figure><img src="image_65.png" /></figure>
<figure><img src="image_66.png" /></figure>

Article

# A Comparative Study of Statistical Methods for Daily Streamflow Estimation at Ungauged Basins in Turkey

Mustafa Utku Yilmaz ¹,* and Bihrat Onoz ²

¹ Department of Civil Engineering, Kirklareli University, 39100 Kirklareli, Turkey

² Department of Civil Engineering, Istanbul Technical University, 34469 Istanbul, Turkey; onoz@itu.edu.tr

* Correspondence: utkuyilmaz@klu.edu.tr; Tel.: +90-288-214-0514

Received: 20 November 2019; Accepted: 7 February 2020; Published: 9 February 2020

**Abstract:** In this study, a comparative evaluation of the statistical methods for daily streamflow estimation at ungauged basins is presented. The single donor station drainage area ratio (DAR) method, the multiple-donor stations drainage area ratio (MDAR) method, the inverse similarity weighted (ISW) method, and its variations with three different power parameters (1, 2, and 3) are applied to the two main subbasins of the Euphrates Basin in Turkey to estimate daily streamflow data. Each station in each basin is considered in turn as the target station where there are no streamflow data. The donor stations are selected based on the physical similarities between the donor and target stations. Then, streamflow data from the most physically similar donor station(s) is transferred to the target station using the statistical methods. In addition, the effect of data preprocessing on the estimation performance of the statistical methods is investigated. The preprocessing discussed in this study is streamflow data smoothing using the two-sided moving average (MA). Three statistical methods using the smoothed data by the MA, named as DAR-MA, MDAR-MA, and ISW-MA, are proposed. The estimation performance of the statistical methods is compared by using daily streamflow data with preprocessing and without preprocessing. The Nash–Sutcliffe efficiency (NSE), the ratio of the root mean square error (RMSE) to the standard deviation of the observed data (RSR), the percent bias (PBIAS), and the coefficient of determination (R²) are used to evaluate the performance of the statistical methods. The results show that MDAR and ISW give improved performances compared to DAR to estimate daily streamflow for 7 out of 8 target stations in the Middle Euphrates Basin and for 4 out of 7 target stations in the Upper Euphrates Basin. Higher NSE values for both MDAR and ISW are mostly obtained with the three most physically similar donor stations in the Middle Euphrates Basin and with the two most physically similar donor stations in the Upper Euphrates Basin. The best statistical method for each target station exhibits slightly greater NSE when the smoothed data by the MA is used for all target stations in the Middle Euphrates Basin and for 6 out of 7 target stations in the Upper Euphrates Basin.

**Keywords:** data preprocessing; donor selection; drainage area ratio; Euphrates basin; moving average; physical similarity; streamflow estimation; ungauged basins

## 1. Introduction

In recent years, several factors, such as climate change, global warming, drought, population growth, and industrialization, have led to a rapid increase in demand for water. Hence, issues related to the planning and implementation of the water budget become important. Measurements and estimates of streamflow play an important role in the stage of the planning and implementation of the water budget. Since drainage basins in many parts of the world are ungauged or poorly gauged, the International Association of Hydrological Sciences (IAHS) launched a scientific decade from 2003 to 2012 on Predictions in Ungauged Basins (PUB) [1]. It was an effort to improve streamflow
estimations for ungauged basins. Streamflow estimation at ungauged and poorly gauged basins is an important issue in growing economies countries such as Turkey because there are a limited number of stations in the streamflow gauging network of Turkey and streamflow estimates are often needed at ungauged basins where water resources projects are planned. Some stations in the river basins of Turkey contain large amounts of missing data during the observation period [2,3]. This lack of adequate data creates significant problems in the water resources projects for Turkey. For these reasons, accurate measurement and analysis of streamflow data and reliable streamflow estimates are needed.

Many methods have been used to improve the reliability and the accuracy in estimations for the development of streamflow estimation methods, and the research in this area still continues. In order to estimate streamflow, several researchers have suggested the artificial intelligence methods such as artificial neural networks [4–6], fuzzy logic [7–9], genetic programming [10–12], and machine learning [13–15]. In addition, artificial intelligence methods have been coupled with the data preprocessing methods to improve streamflow estimation accuracy and reliability in recent studies in the literature [16–19]. For example, Wu and Chau [20] used data preprocessing methods such as moving average (MA) and singular spectrum analysis (SSA) in order to improve the performance of artificial neural networks (ANN). Results showed that the MA was more effective than the SSA when they were coupled with the ANN. Moreover, ANN methods coupled with the MA performed the best among all methods. On the other hand, statistical methods such as interpolation by inverse distance weighted (IDW) [21,22] and kriging [23,24], regression analysis [25,26], flow duration curves [27,28], and information transfer methods [29,30] are widely used in the estimation of streamflow. This study focuses on statistical methods for improving estimation accuracy and reliability for ungauged basins.

Regionalization is a statistical process, which aims to estimate streamflow at ungauged basins. Various regional methods have been used for regional estimation of streamflow for the different time scales (i.e., daily, monthly, or annually) at ungauged basins in the literature [31–35]. Streamflow estimation at ungauged basins where streamflow data are not available requires the transfer of hydrologic information available at a donor station to the target station where only morphological and meteorological characteristics are available [36]. Drainage area ratio (DAR) method is one of the oldest information transfer methods for obtaining streamflow values at the target station from the donor station. This method is straightforward to apply and is in widespread use by hydrologists because it requires no additional information other than the streamflow values at the target station and the drainage areas of the donor and target stations. The DAR method has gained acceptance in Turkey as well, and it is widely used to estimate streamflow for ungauged basins in Turkey [2,3]. In the traditional application of this method, area-normalized streamflow values are transferred from only a single donor station to the target station. In addition to the drainage area, there are some other factors that have a significant influence on the unique streamflow characteristics of a station. Because the DAR method is used with only a single donor station, systematic errors can be encountered in the estimation of a target station [32]. When more streamflow gauging stations are used to estimate streamflow for the target station, this method is referred to as the multiple-donor stations drainage area ratio (MDAR) method [2,32]. The MDAR method assumes that the streamflow estimates at the target station can be computed as the weighted average of the estimates (produced by the DAR method) of the multiple donor stations selected.

The inverse distance weighted (IDW) method is one of the most widely used interpolation methods based on the geographical distance between the donor and target stations [21,22]. This method can be considered as a variant of the DAR method. The IDW method estimates the streamflow value for the target station by taking the geographical distance between the donor station and the target station as the weight. The closer the geographical distance between the donor station and the target station is, the larger the influence on the target station will be. That is, when the distance decreases, the weight coefficient increases. The IDW method, also called an inverse distance to power, is a weighted average interpolator, and the main factor affecting the accuracy of the IDW method is the value of the power parameter. As the power parameter increases, more influence is given to the donor stations
close to the target station. In the literature, the value of the power parameter is commonly chosen as 2, which is known as the inverse distance squared weighted [36]. Alternatively, the inverse similarity weighted (ISW) method [37,38], which is similar to the IDW method, can be applied on the basis of multiple donor stations. Unlike the IDW method, the ISW method uses physical similarity instead of the geographical distance between the target and the donor station. The ISW method with three different power parameters (1, 2, and 3) was used for daily streamflow estimation in this study, and area normalized streamflow values are directly transferred to a target station from multiple donor stations.

The streamflow characteristics at the ungauged basin are directly affected by the donor stations. Therefore, the selection of hydrologically similar donor stations is important for estimating streamflow values at the ungauged basin. In practical applications, the donor station is usually selected as the geographically nearest station to the ungauged basin [36,39,40]. However, the geographical distance may not always be correct for the selection of the donor stations [41,42]. In this study, the physical similarities between the donor and the target station were taken into account when selecting the appropriate donor station for the target station. Physical similarity defines which stations are most similar in terms of some physical characteristics such as drainage area, elevation, precipitation, temperature, latitude, and longitude. According to the physical similarity, donor stations were defined for each target station. This procedure is described in detail in the section "Selection of Donor Stations".

In this study, continuous daily streamflow data were used between 1986–2009 for selected streamflow gauging stations in two subbasins of the Euphrates basin. In order to estimate daily streamflow at the ungauged basin, the single-donor station drainage area ratio (DAR) method, the multiple-donor stations drainage area ratio (MDAR) method, and the inverse similarity weighted (ISW) methods were applied. Three different power parameters (1, 2, and 3) of the ISW method were compared to determine their accuracy and suitability for estimating daily streamflow values. In addition, the daily streamflow data were smoothed with symmetric two-sided moving average (MA) filtering in order to reduce noise. The observed (original) data (without data preprocessing) or the smoothed data (with data preprocessing by the MA) were used as inputs of the statistical methods for estimating daily streamflow values at the target station. In the former case, the estimated daily streamflow values at the target station were compared to the observed (original) daily streamflow values at the target station, while in the latter case, the estimated daily streamflow values at the target station were compared to the observed-MA (smoothed) daily streamflow values at the target station. These two approaches were presented to estimate the daily streamflow values with and without MA. It is believed that the results will help decision makers choose the best one for their objectives. In summary, the major objectives of this study can be listed as follows: 1) to test applicability of the statistical methods to two subbasins of the Euphrates basin in Turkey, 2) to evaluate the success of physical similarity approaches in selecting donor stations in this basin, and 3) to investigate the effect of the statistical methods coupled with the data preprocessing method of moving average (MA) on the accuracy of streamflow estimation.

## 2. Study Area and Data

### 2.1. *Study Area*

Turkey is divided into 25 hydrological river basins, where the Euphrates-Tigris (indicated with basin number 21) is regarded as one single basin (Figure 1). Euphrates-Tigris Basin is located in the eastern part of Turkey with a drainage area of 185,000 km², which is the largest basin of Turkey. It has also nearly 28.5% of the water potential of Turkey. As the biggest water source of the Euphrates-Tigris Basin, the Euphrates River is the longest and one of the most historically significant rivers of the Middle East. The total length of Euphrates is nearly 2800 km, and 40% of its length is in Turkey, 25% is in Syria, and 35% is in Iraq. The Euphrates River consists of two major tributaries, the Karasu River and the Murat River, which both originate in the Eastern Anatolia mountains of Turkey. These two rivers merge near the Keban Dam, which is one of the largest dams of Turkey. The Euphrates River
Basin is subdivided into the Upper Euphrates, the Middle Euphrates, and the Lower Euphrates basins, which have some distinctive physical features. The water regime of the Euphrates River Basin depends heavily on winter rainfalls and spring snowmelt.

<figure><img src="image_67.png" /><figcaption>Figure 1. Hydrological river basins in Turkey.</figcaption></figure>

In order to determine the water potential of the Euphrates-Tigris Basin, a large number of streamflow gauging stations were established on the Euphrates River and its tributaries. However, there was a large amount of missing data in the daily streamflow measurements of some streamflow gauging stations. These missing data lead to significant problems in hydrological modeling studies. Statistical estimation methods require the use of daily streamflow time series obtained from a large number of streamflow gauging stations within the study area. Also, the observation period should be the same for all these streamflow gauging stations. In the Euphrates Basin as a case study, a total of 15 streamflow gauging stations, which have 24 years (1986–2009) of common daily streamflow data, was selected. Eight of these stations are located in the Middle Euphrates Basin, and the other seven stations are located in the Upper Euphrates Basin (Figure 2). Moreover, they are not located downstream of a dam.

<figure><img src="image_68.png" /><figcaption>Figure 2. Location of the two selected basins with streamflow gauging stations.</figcaption></figure>
## 2.2. Hydrological and Meteorological Data

Eight streamflow gauging stations from the Middle Euphrates Basin and seven streamflow gauging stations from the Upper Euphrates Basin were selected for this case study. The stream networks of these two basins and the locations of the selected streamflow gauging stations in each basin are shown in Figure 3. Continuous daily streamflow data of the streamflow gauging stations operated by the General Directorate of State Hydraulic Works (DSI) were used. Each streamflow gauging station contains a 24-year period spanning from 1986 to 2009, and there is no missing data within the streamflow time series. The main characteristics of these streamflow gauging stations are listed in Tables 1 and 2. As shown by Table 1, drainage areas of the stations in the Middle Euphrates Basin vary between 65.3 and 25,515.6 km² whereas their elevations range between 852 and 1810 m above sea level. As shown by Table 2, drainage areas of the stations in the Upper Euphrates Basin vary between 233.2 and 15,562 km² whereas their elevations range between 840 and 1830 m above sea level.

<figure><img src="image_69.png" /><figcaption>Figure 3. Streamflow gauging stations in the study area.</figcaption></figure>

<table><caption>Table 1. Characteristics of streamflow gauging stations in the Middle Euphrates Basin.</caption><thead><tr><td>Station Number</td><td>Drainage Area (km²)</td><td>Elevation (m)</td><td>Long-Term Mean (m³/s)</td><td>Record Period (Years)</td></tr></thead><tbody><tr><td>D21A167</td><td>250</td><td>1650</td><td>3.55</td><td>1986–2009</td></tr><tr><td>D21A169</td><td>276.1</td><td>1600</td><td>3.35</td><td>1986–2009</td></tr><tr><td>D21A213</td><td>65.3</td><td>1810</td><td>0.74</td><td>1986–2009</td></tr><tr><td>E21A002</td><td>25,515.6</td><td>852</td><td>239.82</td><td>1986–2009</td></tr><tr><td>E21A022</td><td>5882.4</td><td>1552</td><td>48.20</td><td>1986–2009</td></tr><tr><td>E21A058</td><td>1577.6</td><td>1310</td><td>18.91</td><td>1986–2009</td></tr><tr><td>E21A064</td><td>2232</td><td>990</td><td>32.97</td><td>1986–2009</td></tr><tr><td>E21A077</td><td>2995.3</td><td>1452</td><td>29.94</td><td>1986–2009</td></tr></tbody></table>
<table><caption>**Table 2.** Characteristics of streamflow gauging stations in the Upper Euphrates Basin.</caption><thead><tr><td>Station Number</td><td>Drainage Area (km²)</td><td>Elevation (m)</td><td>Long-term Mean (m³/s)</td><td>Record Period (Years)</td></tr></thead><tbody><tr><td>D21A001</td><td>233.2</td><td>1830</td><td>2.75</td><td>1986–2009</td></tr><tr><td>D21A193</td><td>518.1</td><td>1000</td><td>6.31</td><td>1986–2009</td></tr><tr><td>E21A033</td><td>3284.8</td><td>875</td><td>89.38</td><td>1986–2009</td></tr><tr><td>E21A051</td><td>8185.6</td><td>1355</td><td>60.23</td><td>1986–2009</td></tr><tr><td>E21A054</td><td>2886</td><td>1675</td><td>19.68</td><td>1986–2009</td></tr><tr><td>E21A056</td><td>15,562</td><td>865</td><td>153.57</td><td>1986–2009</td></tr><tr><td>E21A066</td><td>5430</td><td>840</td><td>78.26</td><td>1986–2009</td></tr></tbody></table>

Basin characteristics such as geographical, topographical, and climate variables were considered for determining the physical similarity between the donor and target stations. Annual mean total precipitation and annual mean temperature were selected as climatic variables. Concurrent precipitation and temperature data of the meteorological stations operated by the Turkish State Meteorological Service (DMI) were used. The annual mean total precipitation and annual mean temperature values for each streamflow gauging station were calculated by the Thiessen polygon method (Figure 4). Thus, annual mean total precipitation and annual mean temperature values of the drainage area represented by each streamflow gauging station were obtained using the precipitation and temperature data of the meteorological stations. Drainage area, elevation, basin slope, and channel length were selected as topographical variables. Basin slope and channel length for the drainage basin of each streamflow gauging station were extracted using geographic information system (GIS) software. The latitude and longitude were selected as geographical variables because geographically nearby streamflow gauging stations could have similarities in hydrological behavior. They were converted to decimal degrees and then used to calculate the similarity coefficient. Since the latitude and longitude define the geographical location of the streamflow gauging stations, these selected basin characteristics combine the physical similarity approach with the geographical proximity approach [43]. Descriptive statistics of the selected physical characteristics are presented in Table 3 for the Middle Euphrates Basin and in Table 4 for the Upper Euphrates Basin.

<figure><img src="image_70.png" /><figcaption>**Figure 4.** Meteorological stations and Thiessen polygons of the study area.</figcaption></figure>
<table><caption>**Table 3.** Statistics of physical characteristics used in the Middle Euphrates Basin.</caption><thead><tr><th>Physical Characteristics</th><th>Maximum</th><th>Minimum</th><th>Mean</th></tr></thead><tbody><tr><td>Drainage Area (km²)</td><td>25,515.6</td><td>65.3</td><td>4849.3</td></tr><tr><td>Elevation (m)</td><td>1810</td><td>852</td><td>1402</td></tr><tr><td>Annual Mean Total Precipitation (mm)</td><td>939.50</td><td>431.20</td><td>679.68</td></tr><tr><td>Annual Mean Temperature (°F)</td><td>53.60</td><td>42.26</td><td>47.29</td></tr><tr><td>Basin Slope (%)</td><td>2.69</td><td>0.19</td><td>1.21</td></tr><tr><td>Channel Length (km)</td><td>565.11</td><td>14.75</td><td>142.27</td></tr><tr><td>Latitude (°)</td><td>39.54</td><td>38.69</td><td>39.19</td></tr><tr><td>Longitude (°)</td><td>42.78</td><td>39.93</td><td>41.58</td></tr></tbody></table>

<table><caption>**Table 4.** Statistics of physical characteristics used in the Upper Euphrates Basin.</caption><thead><tr><th>Physical Characteristics</th><th>Maximum</th><th>Minimum</th><th>Mean</th></tr></thead><tbody><tr><td>Drainage Area (km²)</td><td>15,562</td><td>233.2</td><td>5157.1</td></tr><tr><td>Elevation (m)</td><td>1830</td><td>840</td><td>1205.7</td></tr><tr><td>Annual Mean Total Precipitation (mm)</td><td>840.17</td><td>374.90</td><td>524.15</td></tr><tr><td>Annual Mean Temperature (°F)</td><td>56.48</td><td>42.26</td><td>48.52</td></tr><tr><td>Basin Slope (%)</td><td>2.82</td><td>0.16</td><td>0.96</td></tr><tr><td>Channel Length (km)</td><td>381.60</td><td>25.10</td><td>161.19</td></tr><tr><td>Latitude (°)</td><td>40.11</td><td>38.86</td><td>39.44</td></tr><tr><td>Longitude (°)</td><td>41.39</td><td>38.41</td><td>39.79</td></tr></tbody></table>

### 3. Methods

Estimation of daily streamflow time series at the target station consists of the following steps: (1) the selection of hydrologically similar donor stations to the target station and (2) the transfer of the daily streamflow time series from the donor station to target station by using statistical streamflow transfer methods. Proposed flowcharts for streamflow estimation at the target station are illustrated in Figure 5.

<figure><img src="image_71.png" /><figcaption>**Figure 5.** Proposed flowcharts for streamflow estimation at the target station.</figcaption></figure>
## 3.1. Statistical Information Transfer Methods

For daily streamflow estimation at the target stations, three statistical streamflow transfer methods are considered in this study. These methods include the single-donor station drainage area ratio (DAR) method, the multiple-donor stations drainage area ratio (MDAR) method, and the inverse similarity weighted (ISW) method. Moreover, the variations of the ISW method, in which constant power parameters are modified, are utilized as well. The DAR and the MDAR methods were based on the drainage area of the stations, whereas the ISW method was based on the physical similarity between the donor and target stations. Each statistical method is briefly described below.

### 3.1.1. Drainage Area Ratio (DAR) Method

The DAR method [39,40] assumes that the streamflow per unit drainage area for the target station equals that at the streamflow gauging station used as a donor station for a given day, as described in Equation (1).

$$Q_{target} = \frac{A_{target}}{A_{donor}} Q_{donor} \quad (1)$$

where $Q_{target}$ is the daily streamflow for the target station, $Q_{donor}$ is the streamflow at a donor station, and $A_{target}$ and $A_{donor}$ are the drainage areas for the target station and the donor station, respectively.

### 3.1.2. Multiple-Donor Stations Drainage Area Ratio (MDAR) Method

The MDAR method [2,32] generates the streamflow estimation at the target station as the weighted average of the DAR method estimations from the donor stations. The streamflow at the target station from the n donor stations can be calculated for a given day using Equation (2).

$$Q_{target} = \frac{\sum_{i=1}^{n} w_i \hat{Q}_{donor_i}}{\sum_{i=1}^{n} w_i} \quad (2)$$

where $w_i$ is the weight of the donor station i on the target station, $\hat{Q}_{donor_i}$ is the daily streamflow estimations from each donor station, and n is the total number of the donor stations. The values of the weights in Equation (2) which show the similarity between the target station and donor station can be calculated as Equation (3).

$$w_i = \frac{\frac{1}{d_i}}{\sum_{i=1}^{n} \frac{1}{d_i}} \quad (3)$$

where $d_i$ is the similarity distance between the target station and donor station i. The drainage area is frequently considered as the most important variable in many hydrological regionalization studies [39,40,44]. Moreover, the drainage area is also the only scaling factor used in the DAR method for streamflow estimation at the target stations. Therefore, it was used as the similarity distance in this study. It can be calculated using Equation (4).

$$d_i = |A_{target} - A_{donor_i}| \quad (4)$$

where $A_{donor_i}$ is the drainage area of the donor station i and $A_{target}$ is the drainage area of the target station.

### 3.1.3. Inverse Similarity Weighted (ISW) Method

The ISW method [37,38] estimates streamflow values at the target station as the weighted average of the streamflow values at n donor stations. The weights are inversely proportional to the power of
physical similarity from the target station. The mathematical expression of the ISW method is given by Equation (5).

$$q_{\text{target}} = \sum_{i=1}^{n} w_i q_{\text{donor}_i} \quad \text{and} \quad Q_{\text{target}} = q_{\text{target}} A_{\text{target}} \qquad (5)$$

where $q_{\text{target}}$ is the area normalized streamflow (m³/s/km²) at the target station and $q_{\text{donor}_i}$ is the area normalized streamflow (m³/s/km²) at the donor station $i$. The weights $w_i$ based on physical similarity can be calculated for all donor stations using Equation (6). The sum of the weights assigned to each donor station is equal to 1.

$$w_i = \frac{\frac{1}{s_i^p}}{\sum_{i=1}^{n} \frac{1}{s_i^p}} \quad \text{and} \quad \sum_{i=1}^{n} w_i = 1 \qquad (6)$$

where $s_i$ is the similarity coefficient between the target station and donor station $i$ and where the exponent $p$ is called a power parameter ($p > 0$). In this study, the estimation performance of the ISW method was evaluated using different power parameters from 1 to 3. For power parameters of 1, 2, and 3, the ISW method was referred to as ISW1, ISW2, and ISW3, respectively.

The similarity coefficient, $s$, is used to define the physical similarity between the target station and the donor station, which is calculated using Equation (7) [43,45]. Drainage area, elevation, annual mean total precipitation, annual mean temperature, basin slope, channel length, latitude, and longitude were considered as the basin characteristics in order to measure the physical similarity between the donor station and the target station. The station with the lowest similarity coefficient was selected as the donor station. The similarity coefficient was used both to select the donor stations and to transfer streamflow from several donor stations as the weight.

$$s = \sum_{i=1}^{k} \frac{|X_{i}^{\text{donor}} - X_{i}^{\text{target}}|}{\max(X_{i}) - \min(X_{i})} \qquad (7)$$

where $i$ indicates one of a total of $k$ selected basin characteristics; $X_{i}^{\text{donor}}$ and $X_{i}^{\text{target}}$ are the values of basin characteristic $i$ for the donor station and the target station, respectively; and $\max(X_i)$ and $\min(X_i)$ are the maximum and the minimum values of basin characteristics over the set of stations considered, respectively.

## 3.2. Selection of Donor Stations

For transferring the streamflow to the target (ungauged) station, the streamflow values of the donor stations are used. Therefore, the selection of the donor stations is an important step in estimating streamflow at the target station. In this study, the physical similarity approach was considered to identify donor stations. In the physical similarity approach, the station that minimizes the similarity coefficient defined in Equation (7) was used as the donor station. That is, the best donor station was given to the station having the smallest $s$ value. Although all stations are gauged, initially, each station in each basin was considered in turn as a target station and daily streamflow time series were estimated for all stations assumed as a target station within each basin. Subsequently, their actual streamflow time series were used in order to evaluate the performance of the streamflow estimation. When using one donor station, the most physically similar station was identified for each target station. On the other hand, when using more than one donor station, the two or three most physically similar stations were identified for each target station (Figure 5). Therefore, the DAR method was applied to each station using the most physically similar station as the donor station, while the MDAR and ISW methods were applied to each station for two different cases: 1) using the two most physically similar stations as the donor stations and 2) using the three most physically similar stations as the donor stations. The results for these two different cases were compared with the original observations to evaluate which one provides better estimation performance.
### 3.3. Data Preprocessing

The streamflow data may contain possible errors, and these errors are collectively referred to as noise. As the noise in the data increases, reliable results will be difficult to achieve. In this study, data preprocessing was conducted to remove noise and to improve the reliability of daily streamflow estimates. The data preprocessing discussed here was daily streamflow data smoothing using the moving average (MA). Each daily streamflow time series of all stations was smoothed by a centered (or symmetric two-sided) moving average of length $m = 2k + 1$, i.e., MA($m$), and then, the smoothed streamflow time series were used into the statistical methods. Hereafter, the statistical methods, DAR, MDAR, and ISW are referred to as DAR-MA, MDAR-MA, and ISW-MA, respectively. A centered moving average smooths data by replacing each observed daily streamflow value with the average of the current day, previous day, and subsequent days and is defined as Equation (8). For example, a centered moving average of length $m = 3$ (hence $k = 1$), i.e., MA(3) with equal weights, replaces the observed daily streamflow value $x_t$ at time $t$ with the averages of $x_{t-1}$, $x_t$, and $x_{t+1}$.

$$x_t^* = \frac{1}{2k+1} \sum_{i=-k}^{k} x_{t+i} \quad (8)$$

where $x_t^*$ is the smoothed streamflow value at time $t$ and $m = 2k + 1$ is the number of observed values that are averaged.

In order to smooth daily streamflow data, MA was applied with lengths of 3, 5, 7, 9, and 11 days in this study, and then, it was seen that the larger the length $m = 2k + 1$, the more the streamflow peaks (maximum values) and streamflow valleys (minimum values) were smoothed out. The peaks and valleys of streamflow are not well represented by the relatively high length of MA(5), MA(7), MA(9), and MA(11). For a better representation of streamflow peaks and valleys, MA(5), MA(7), MA(9), and MA(11) were excluded from the rest of the study.

### 3.4. Evaluation Criteria

A jackknife (leave one out) procedure was used for evaluating the performance of each method. In this procedure, each station in each basin was considered in turn as a target station and was removed from the database. This procedure was repeated for all stations considered for this study.

The performance of each statistical method was evaluated in terms of the Nash–Sutcliffe efficiency (NSE) [46], the ratio of the root mean square error (RMSE) to the standard deviation of the observed data (RSR) [47], the percent bias (PBIAS), and the coefficient of determination ($R^2$) between the estimated and observed streamflow values. NSE, RSR, PBIAS, and $R^2$ were calculated as follows:

$$NSE = 1 - \frac{\sum_{i=1}^{n} (X_i^{obs} - X_i^{est})^2}{\sum_{i=1}^{n} (X_i^{obs} - \bar{X}^{obs})^2} \quad (9)$$

$$RSR = \frac{RMSE}{\sqrt{\frac{\sum_{i=1}^{n} (X_i^{obs} - \bar{X}^{obs})^2}{n}}} = \frac{\sqrt{\sum_{i=1}^{n} (X_i^{obs} - X_i^{est})^2}}{\sqrt{\sum_{i=1}^{n} (X_i^{obs} - \bar{X}^{obs})^2}} = \sqrt{1 - NSE} \quad (10)$$

$$PBIAS = \frac{\sum_{i=1}^{n} (X_i^{obs} - X_i^{est})}{\sum_{i=1}^{n} (X_i^{obs})} \times 100 \quad (11)$$
$$R^2 = \left( \frac{\sum_{i=1}^{n} (X_{i}^{\text{obs}} - \bar{X}^{\text{obs}})(X_{i}^{\text{est}} - \bar{X}^{\text{est}})}{\sqrt{\sum_{i=1}^{n} (X_{i}^{\text{obs}} - \bar{X}^{\text{obs}})^2} \sqrt{\sum_{i=1}^{n} (X_{i}^{\text{est}} - \bar{X}^{\text{est}})^2}} \right)^2 \quad (12)$$

where $X_{i}^{\text{obs}}$ is the ith observed daily streamflow value; $X_{i}^{\text{est}}$ is the ith estimated daily streamflow value; $\bar{X}^{\text{obs}}$ and $\bar{X}^{\text{est}}$ are the mean of observed and estimated daily streamflow values, respectively; and n is the total number of observed daily streamflow values.

The NSE values range between -∞ and +1, where a value of 1 indicates a perfect agreement between estimated and observed streamflow values. The values closer to 1 indicate an increasingly better agreement, whereas the values far from 1 indicate an increasingly poor agreement. The RSR standardizes RMSE using the standard deviation of the observed data. The RSR varies from the optimal value of 0 to a large positive value. The lower the RSR, the better the performance of the method [47]. PBIAS measures the average tendency of the estimated values to be larger or smaller than corresponding observed values. The optimal value of PBIAS is 0, and the closer it is to 0, the more accurate the estimated values are to the observed values. Negative PBIAS values indicate overestimation, while positive PBIAS values indicate underestimation [47]. The coefficient of determination (R²) is the square of the correlation coefficient according to Pearson. The R² values range from 0 to 1, with higher values indicating better agreement between estimated and observed values. Generally, R² values greater than 0.5 are considered acceptable [47].

Moriasi et al. [47] suggested performance ratings of recommended statistics such as NSE, RSR, and PBIAS for monthly streamflow. According to Moriasi et al. [47], the performance of the method is considered satisfactory when the NSE is greater than 0.5, the RSR is less than 0.7, and the PBIAS ranges are less than ±25%. However, NSE values lower than 0.5 for daily streamflow can still be considered satisfactory [48]. Therefore, some of the constraints for the recommended statistics can be relaxed for daily streamflow. In this study, the adjusted performance ratings of the NSE and PBIAS statistics for the daily time scale proposed by Kalin et al. [49] were used to evaluate the performance of the statistical methods (Table 5).

<table><caption>**Table 5.** Performance ratings of the Nash–Sutcliffe efficiency (NSE) and percent bias (PBIAS) statistics for daily streamflow.</caption><thead><tr><td>Performance Rating</td><td>NSE</td><td>abs(PBIAS) %</td></tr></thead><tbody><tr><td>Very good</td><td>NSE ≥ 0.7</td><td>abs(PBIAS) ≤ 25</td></tr><tr><td>Good</td><td>0.5 ≤ NSE &lt; 0.7</td><td>25 &lt; abs(PBIAS) ≤ 50</td></tr><tr><td>Satisfactory</td><td>0.3 ≤ NSE &lt; 0.5</td><td>50 &lt; abs(PBIAS) ≤ 70</td></tr><tr><td>Unsatisfactory</td><td>NSE &lt; 0.3</td><td>abs(PBIAS) &gt; 70</td></tr></tbody></table>

# 4. Results and Discussion

## 4.1. Middle Euphrates Basin

The statistical methods described above were applied on each of eight streamflow gauging stations in the Middle Euphrates Basin, which were considered in turn as the target station. The donor station or stations based on physical similarity was selected to transfer daily streamflow data to the target station. Table 6 shows the sequence of the donor stations for each target station in the Middle Euphrates Basin which was determined according to the similarity coefficient.

In order to estimate daily streamflow at the target stations, the DAR method was applied by using the most physically similar station to each target station. In order to test the applicability of the donor station selection criteria for the study area, the NSE values were determined for the donor stations identified by the physical similarity and compared with the NSE values obtained from the donor stations traditionally selected as the geographically nearest stations. On the other hand, the
MDAR and ISW methods were applied by using the two and the three most physically similar stations to each target station. In order to test the effect of different power parameter selection in the use of the ISW method on the accuracy of daily streamflow estimation, the ISW method was applied with power parameters of 1, 2, and 3. In addition, the comparisons between the statistical methods with and without the MA were carried out to indicate the effectiveness of the MA-based preprocessing on the accuracy of daily streamflow estimation. Daily streamflow values estimated using both observed and smoothed data from the donor stations were compared with observed data at the target station. According to the selection criteria of the donor stations, the NSE values obtained for the target stations in the Middle Euphrates Basin are given in Table 7 for DAR and MDAR and in Table 8 for ISW with three different power parameters (1, 2, and 3). Performance evaluations of the best statistical method without and with the MA for each target station in the Middle Euphrates Basin are presented in Table 9 and in Table 10, respectively.

<table><caption>**Table 6.** Physically similar donor stations for each target station in the Middle Euphrates Basin.</caption><thead><tr><th rowspan="2">Target Station</th><th colspan="7">Donor Station</th></tr><tr><th>1st</th><th>2nd</th><th>3rd</th><th>4th</th><th>5th</th><th>6th</th><th>7th</th></tr></thead><tbody><tr><td>D21A167</td><td>D21A213</td><td>D21A169</td><td>E21A022</td><td>E21A058</td><td>E21A077</td><td>E21A064</td><td>E21A002</td></tr><tr><td>D21A169</td><td>E21A058</td><td>E21A077</td><td>D21A167</td><td>D21A213</td><td>E21A022</td><td>E21A064</td><td>E21A002</td></tr><tr><td>D21A213</td><td>D21A167</td><td>E21A022</td><td>D21A169</td><td>E21A058</td><td>E21A077</td><td>E21A064</td><td>E21A002</td></tr><tr><td>E21A002</td><td>E21A064</td><td>E21A058</td><td>E21A077</td><td>E21A022</td><td>D21A169</td><td>D21A167</td><td>D21A213</td></tr><tr><td>E21A022</td><td>D21A213</td><td>E21A077</td><td>D21A167</td><td>E21A058</td><td>D21A169</td><td>E21A064</td><td>E21A002</td></tr><tr><td>E21A058</td><td>E21A077</td><td>D21A169</td><td>E21A064</td><td>E21A022</td><td>D21A167</td><td>D21A213</td><td>E21A002</td></tr><tr><td>E21A064</td><td>E21A058</td><td>E21A077</td><td>D21A169</td><td>E21A002</td><td>E21A022</td><td>D21A167</td><td>D21A213</td></tr><tr><td>E21A077</td><td>E21A058</td><td>D21A169</td><td>E21A022</td><td>E21A064</td><td>D21A213</td><td>D21A167</td><td>E21A002</td></tr></tbody></table>

<table><caption>**Table 7.** NSE values for the drainage area ratio (DAR) and multiple-donor stations drainage area ratio (MDAR) methods in the Middle Euphrates Basin.</caption><thead><tr><th rowspan="2">Target Station</th><th>The Geographically Nearest Donor Station</th><th>The Most Physically Similar Donor Station</th><th>Two Most Physically Similar Donor Stations</th><th>Three Most Physically Similar Donor Stations</th></tr><tr><th>DAR</th><th>DAR</th><th>MDAR</th><th>MDAR</th></tr></thead><tbody><tr><td>D21A167</td><td>0.313</td><td>-0.187</td><td>0.312</td><td>0.316</td></tr><tr><td>D21A169</td><td>0.581</td><td>0.850</td><td>0.852</td><td>0.595</td></tr><tr><td>D21A213</td><td>0.718</td><td>0.354</td><td>0.369</td><td>0.608</td></tr><tr><td>E21A002</td><td>-0.171</td><td>-0.171</td><td>0.433</td><td>0.729</td></tr><tr><td>E21A022</td><td>-0.205</td><td>-0.205</td><td>0.622</td><td>0.718</td></tr><tr><td>E21A058</td><td>0.826</td><td>0.814</td><td>0.894</td><td>0.839</td></tr><tr><td>E21A064</td><td>0.649</td><td>0.724</td><td>0.694</td><td>0.706</td></tr><tr><td>E21A077</td><td>0.300</td><td>0.647</td><td>0.651</td><td>0.781</td></tr></tbody></table>

<table><caption>**Table 8.** NSE values for the inverse similarity weighted (ISW) with powers of 1, 2, and 3 in the Middle Euphrates Basin.</caption><thead><tr><th rowspan="2">Target Station</th><th colspan="3">Two Most Physically Similar Donor Stations</th><th colspan="3">Three Most Physically Similar Donor Stations</th></tr><tr><th>ISW1</th><th>ISW2</th><th>ISW3</th><th>ISW1</th><th>ISW2</th><th>ISW3</th></tr></thead><tbody><tr><td>D21A167</td><td>0.136</td><td>0.060</td><td>-0.009</td><td>0.317</td><td>0.184</td><td>0.065</td></tr><tr><td>D21A169</td><td>0.844</td><td>0.849</td><td>0.852</td><td>0.816</td><td>0.825</td><td>0.833</td></tr><tr><td>D21A213</td><td>0.501</td><td>0.463</td><td>0.429</td><td>0.608</td><td>0.556</td><td>0.502</td></tr><tr><td>E21A002</td><td>0.418</td><td>0.396</td><td>0.374</td><td>0.696</td><td>0.664</td><td>0.630</td></tr><tr><td>E21A022</td><td>0.432</td><td>0.375</td><td>0.316</td><td>0.636</td><td>0.591</td><td>0.535</td></tr><tr><td>E21A058</td><td>0.893</td><td>0.892</td><td>0.892</td><td>0.904</td><td>0.907</td><td>0.906</td></tr><tr><td>E21A064</td><td>0.696</td><td>0.701</td><td>0.706</td><td>0.712</td><td>0.713</td><td>0.714</td></tr><tr><td>E21A077</td><td>0.641</td><td>0.651</td><td>0.656</td><td>0.774</td><td>0.743</td><td>0.716</td></tr></tbody></table>
<table><caption>**Table 9.** Performance evaluation of the best statistical method without the moving average (MA) for each target station in the Middle Euphrates Basin.</caption><thead><tr><th>Target Station</th><th>Method</th><th>NSE</th><th>RSR</th><th>PBIAS</th></tr></thead><tbody><tr><td>D21A167</td><td>ISW1</td><td>0.317 <sup>3</sup></td><td>0.827</td><td>23.265 <sup>1</sup></td></tr><tr><td>D21A169</td><td>ISW3</td><td>0.852 <sup>1</sup></td><td>0.385</td><td>6.459 <sup>1</sup></td></tr><tr><td>D21A213</td><td>ISW1</td><td>0.608 <sup>2</sup></td><td>0.626</td><td>-6.122 <sup>1</sup></td></tr><tr><td>E21A002</td><td>MDAR</td><td>0.729 <sup>1</sup></td><td>0.520</td><td>-30.114 <sup>2</sup></td></tr><tr><td>E21A022</td><td>MDAR</td><td>0.718 <sup>1</sup></td><td>0.531</td><td>-39.089 <sup>2</sup></td></tr><tr><td>E21A058</td><td>ISW2</td><td>0.907 <sup>1</sup></td><td>0.304</td><td>3.123 <sup>1</sup></td></tr><tr><td>E21A064</td><td>DAR</td><td>0.724 <sup>1</sup></td><td>0.525</td><td>18.865 <sup>1</sup></td></tr><tr><td>E21A077</td><td>MDAR</td><td>0.781 <sup>1</sup></td><td>0.468</td><td>-11.015 <sup>1</sup></td></tr></tbody></table>

<table><caption>**Table 10.** Performance evaluation of the best statistical method with the MA for each target station in the Middle Euphrates Basin.</caption><thead><tr><th>Target Station</th><th>Method</th><th>NSE</th><th>RSR</th><th>PBIAS</th></tr></thead><tbody><tr><td>D21A167</td><td>ISW1-MA</td><td>0.326 <sup>3</sup></td><td>0.821</td><td>23.265 <sup>1</sup></td></tr><tr><td>D21A169</td><td>ISW3-MA</td><td>0.876 <sup>1</sup></td><td>0.352</td><td>6.459 <sup>1</sup></td></tr><tr><td>D21A213</td><td>ISW1-MA</td><td>0.618 <sup>2</sup></td><td>0.618</td><td>-6.122 <sup>1</sup></td></tr><tr><td>E21A002</td><td>MDAR-MA</td><td>0.767 <sup>1</sup></td><td>0.483</td><td>-30.114 <sup>2</sup></td></tr><tr><td>E21A022</td><td>MDAR-MA</td><td>0.725 <sup>1</sup></td><td>0.524</td><td>-39.089 <sup>2</sup></td></tr><tr><td>E21A058</td><td>ISW2-MA</td><td>0.921 <sup>1</sup></td><td>0.282</td><td>3.123 <sup>1</sup></td></tr><tr><td>E21A064</td><td>DAR-MA</td><td>0.744 <sup>1</sup></td><td>0.506</td><td>18.865 <sup>1</sup></td></tr><tr><td>E21A077</td><td>MDAR-MA</td><td>0.817 <sup>1</sup></td><td>0.428</td><td>-11.015 <sup>1</sup></td></tr></tbody></table>

For 2 out of 8 target stations (i.e., E21A002 and E21A022), the geographically nearest and the most physically similar station were the same. For 3 out of the remaining 6 target stations (i.e., D21A169, E21A064, and E21A077), higher NSE values were obtained using DAR with the most physically similar station as the donor station. For 3 out of the target stations (i.e., D21A167, D21A213, and E21A058), higher NSE values were obtained using DAR with the geographically nearest station as the donor station. According to these results, for half of the target stations in the study area, the geographical distance seems to be a good selection criterion as the donor station; however, for the remaining half of the target stations, geographical distance cannot identify the best donor station. Therefore, donor station selection criteria can provide different estimated results that vary from basin to basin.

As can be seen in Table 7, for all target stations other than E21A064, higher NSE values were obtained with MDAR compared with DAR. Especially for D21A167, E21A002, and E21A022, negative NSE values obtained with DAR using the most physically similar donor station improved considerably when the three most physically similar donor stations were used with MDAR. The performance of the DAR method was unsatisfactory for D21A167, E21A002, and E21A022. This was mostly due to the significant increase in the drainage area ratio between the donor and target stations. D21A213 has the smallest drainage area (65.3 km²) in the Middle Euphrates Basin and was determined as the most physically similar donor station for both D21A167 (250 km²) and E21A022 (5882.4 km²). Moreover, E21A002 has the largest drainage area (25,515.6 km²) in the Middle Euphrates Basin. Its drainage area is more than four times the next largest station. For all target stations other than D21A169 and E21A058, MDAR using of the three most physically similar donor stations produced better NSE values than that using the two most physically similar donor stations. For D21A169, the NSE value decreased from 0.852 to 0.595 when the three most physically similar donor stations were used instead of the two most physically similar donor stations (see Table 7). In case of the use of the three most physically

<footnote n="1">¹ Very good, ² good, ³ satisfactory, and ⁴ unsatisfactory.</footnote>

<footnote n="1">¹ Very good, ² good, ³ satisfactory, and ⁴ unsatisfactory.</footnote>
similar donor stations, the third most physically similar donor station for D21A169 was determined as D21A167. The NSE value obtained for D21A169 using the DAR method and utilizing D21A167 was lower than the NSE values obtained from the other two donor stations (i.e., E21A058 and E21A077). The drainage area of the donor station D21A167 is very close to the target station D21A169. On the other hand, the drainage areas of the other two donor stations are much larger than D21A169. Hence, the weight of donor station D21A167 for streamflow estimation of D21A169 is significantly larger compared to the other two. Consequently, the NSE value obtained for D21A169 using MDAR with the three most physically similar donor stations is predominantly influenced by donor station D21A167. Similarly, for E21A058, the decrease in the NSE (i.e., from 0.894 to 0.839) was due to the same reason as for D21A169.

As can be seen in Table 8, in case of the use of the two most physically similar donor stations, the best performance results were obtained with ISW1 for 5 out of 8 target stations (i.e., D21A167, D21A213, E21A002, E21A022, and E21A058). On the other hand, in case of the use of the three most physically similar donor stations, the best performance results were obtained with ISW1 for 5 out of 8 target stations (i.e., D21A167, D21A213, E21A002, E21A022, and E21A077). In both cases, the most reasonable estimation results were mostly obtained when ISW1 was applied instead of ISW2 and ISW3. Moreover, the NSE values mostly improved when the three most physically similar donor stations were used instead of the two most physically similar donor stations.

As can be seen in Table 9, for all target stations other than for E21A064, the MDAR and the ISW methods resulted in higher NSEs compared to the DAR method. For 6 out of 8 target stations, the results can be rated as “very good” for NSE according to the performance ratings in Table 5. For 7 out of 8 target stations, the RSR values were considered satisfactory (i.e., less than 0.7) according to the performance ratings recommended by Moriasi et al. [47]. The negative PBIAS values for D21A213, E21A002, E21A022, and E21A077 demonstrate that the method overestimated daily streamflow, while positive PBIAS values for D21A167, D21A169, E21A058, and E21A064 demonstrate underestimation. For all target stations, the statistical methods with the MA tend to achieve slightly higher NSE values. However, the PBIAS values of the target stations did not change when the statistical methods with the MA were used.

For the target station E21A058 as the example, the estimated daily streamflow values from the statistical methods without MA were compared to the observed (original) daily streamflow values in the hydrograph and scatter plots in Figure 6. The remarkably better agreement between observed and estimated daily streamflow values by three statistical methods (i.e, DAR, MDAR, and ISW) was obtained for E21A058 compared to the other target stations in the Middle Euphrates Basin. ISW2 gave a coefficient of determination (R²) of 0.91, which was higher than the R² values of 0.87 and 0.90 obtained by using the DAR and MDAR, respectively. The NSE values for these methods ranged from 0.814 to 0.907, and the best NSE value was achieved by ISW2. The best NSE performance for E21A058 was obtained using ISW2 with the three most physically similar donor stations.

On the other hand, for the target station E21A058 as the example, the estimated daily streamflow values from the statistical methods with MA were compared with observed-MA (smoothed) daily streamflow values in the hydrograph and scatter plots in Figure 7. The statistical methods with MA performed slightly better for E21A058.
<figure><img src="image_72.png" /><figcaption>**Figure 6.** Comparison of daily streamflow estimated using DAR (top row), MDAR (mid row), and ISW2 (bottom row) with observed (original) daily streamflow at station E21A058 in the Middle Euphrates Basin.</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>**Figure 7.** Comparison of daily streamflow estimated using DAR-MA (top row), MDAR-MA (mid row), and ISW2-MA (bottom row) with observed-MA (smoothed) daily streamflow at station E21A058 in the Middle Euphrates Basin.</figcaption></figure>
## 4.2. Upper Euphrates Basin

Using the same procedure applied for the Middle Euphrates Basin, the statistical methods were applied on each of seven streamflow gauging stations in the Upper Euphrates Basin for the purpose of estimating daily streamflow. Table 11 shows the sequence of the donor stations for each target station in the Upper Euphrates Basin which was determined according to the similarity coefficient.

<table><caption>Table 11. Physically similar donor stations for each target station in the Upper Euphrates Basin.</caption><thead><tr><th rowspan="2">Target Station</th><th colspan="6">Donor Station</th></tr><tr><th>1st</th><th>2nd</th><th>3rd</th><th>4th</th><th>5th</th><th>6th</th></tr></thead><tbody><tr><td>D21A001</td><td>E21A054</td><td>E21A051</td><td>D21A193</td><td>E21A033</td><td>E21A033</td><td>E21A066</td></tr><tr><td>D21A193</td><td>E21A033</td><td>E21A066</td><td>E21A056</td><td>E21A051</td><td>E21A051</td><td>E21A054</td></tr><tr><td>E21A033</td><td>E21A066</td><td>D21A193</td><td>E21A051</td><td>E21A056</td><td>E21A056</td><td>D21A001</td></tr><tr><td>E21A051</td><td>E21A054</td><td>E21A056</td><td>E21A066</td><td>E21A033</td><td>E21A033</td><td>D21A193</td></tr><tr><td>E21A054</td><td>E21A051</td><td>D21A001</td><td>E21A033</td><td>E21A056</td><td>E21A056</td><td>D21A193</td></tr><tr><td>E21A056</td><td>E21A051</td><td>E21A066</td><td>E21A033</td><td>D21A193</td><td>D21A193</td><td>D21A001</td></tr><tr><td>E21A066</td><td>E21A033</td><td>E21A051</td><td>E21A056</td><td>D21A193</td><td>D21A193</td><td>D21A001</td></tr></tbody></table>

According to the selection criteria of the donor stations, the NSE values obtained for the target stations in the Upper Euphrates Basin were given in Table 12 for DAR and MDAR and in Table 13 for ISW with three different power parameters (1, 2, and 3). Performance evaluations of the best statistical method without and with the MA for each target station in the Upper Euphrates Basin are presented in Table 14 and in Table 15, respectively.

For 2 out of the 3 target stations (i.e., E21A054 and E21A056) for which the geographically nearest and the most physically similar stations were not the same, higher NSE values were obtained using DAR with the most physically similar station (see Table 12). For E21A054 and E21A056, the NSE values improved considerably when the most physically similar station was used. These results indicate that the physical similarity may be a better selection criterion for the donor station in the study area compared to the geographical distance between the stations.

<table><caption>Table 12. NSE values for the DAR and MDAR methods in the Upper Euphrates Basin.</caption><thead><tr><th rowspan="2">Target Station</th><th>The Geographically Nearest Donor Station</th><th>The Most Physically Similar Donor Station</th><th>Two Most Physically Similar Donor Stations</th><th>Three Most Physically Similar Donor Stations</th></tr><tr><th>DAR</th><th>DAR</th><th>MDAR</th><th>MDAR</th></tr></thead><tbody><tr><td>D21A001</td><td>0.619</td><td>0.619</td><td>0.630</td><td>-0.328</td></tr><tr><td>D21A193</td><td>0.299</td><td>-0.646</td><td>-0.053</td><td>0.089</td></tr><tr><td>E21A033</td><td>0.446</td><td>0.446</td><td>0.377</td><td>0.319</td></tr><tr><td>E21A051</td><td>0.883</td><td>0.883</td><td>0.932</td><td>0.362</td></tr><tr><td>E21A054</td><td>0.159</td><td>0.885</td><td>0.585</td><td>-5.128</td></tr><tr><td>E21A056</td><td>-3.199</td><td>0.769</td><td>0.415</td><td>-0.586</td></tr><tr><td>E21A066</td><td>-0.067</td><td>-0.067</td><td>0.698</td><td>0.732</td></tr></tbody></table>

<table><caption>Table 13. NSE values for ISW with powers of 1, 2, and 3 in the Upper Euphrates Basin.</caption><thead><tr><th rowspan="2">Target Station</th><th colspan="3">Two Most Physically Similar Donor Stations</th><th colspan="3">Three Most Physically Similar Donor Stations</th></tr><tr><th>ISW1</th><th>ISW2</th><th>ISW3</th><th>ISW1</th><th>ISW2</th><th>ISW3</th></tr></thead><tbody><tr><td>D21A001</td><td>0.633</td><td>0.629</td><td>0.626</td><td>0.597</td><td>0.626</td><td>0.631</td></tr><tr><td>D21A193</td><td>0.015</td><td>-0.085</td><td>-0.185</td><td>0.307</td><td>0.217</td><td>0.108</td></tr><tr><td>E21A033</td><td>0.416</td><td>0.439</td><td>0.444</td><td>0.353</td><td>0.412</td><td>0.435</td></tr><tr><td>E21A051</td><td>0.931</td><td>0.928</td><td>0.921</td><td>0.811</td><td>0.877</td><td>0.905</td></tr><tr><td>E21A054</td><td>0.753</td><td>0.775</td><td>0.794</td><td>0.244</td><td>0.580</td><td>0.719</td></tr><tr><td>E21A056</td><td>0.404</td><td>0.541</td><td>0.640</td><td>-0.788</td><td>-0.306</td><td>0.083</td></tr><tr><td>E21A066</td><td>0.479</td><td>0.191</td><td>0.033</td><td>0.645</td><td>0.331</td><td>0.097</td></tr></tbody></table>
<table><caption>**Table 14.** Performance evaluation of the best statistical method without the MA for each target station in the Upper Euphrates Basin.</caption><thead><tr><th>Target Station</th><th>Method</th><th>NSE</th><th>RSR</th><th>PBIAS</th></tr></thead><tbody><tr><td>D21A001</td><td>ISW1</td><td>0.633<sup>2</sup></td><td>0.606</td><td>40.537<sup>2</sup></td></tr><tr><td>D21A193</td><td>ISW1</td><td>0.307<sup>3</sup></td><td>0.832</td><td>-51.370<sup>3</sup></td></tr><tr><td>E21A033</td><td>DAR</td><td>0.446<sup>3</sup></td><td>0.744</td><td>47.035<sup>2</sup></td></tr><tr><td>E21A051</td><td>MDAR</td><td>0.932<sup>1</sup></td><td>0.261</td><td>-10.006<sup>1</sup></td></tr><tr><td>E21A054</td><td>DAR</td><td>0.885<sup>1</sup></td><td>0.339</td><td>-7.892<sup>1</sup></td></tr><tr><td>E21A056</td><td>DAR</td><td>0.769<sup>1</sup></td><td>0.480</td><td>25.436<sup>2</sup></td></tr><tr><td>E21A066</td><td>MDAR</td><td>0.732<sup>1</sup></td><td>0.518</td><td>-22.121<sup>1</sup></td></tr></tbody></table>

<table><caption>**Table 15.** Performance evaluation of the best statistical method with the MA for each target station in the Upper Euphrates Basin.</caption><thead><tr><th>Target Station</th><th>Method</th><th>NSE</th><th>RSR</th><th>PBIAS</th></tr></thead><tbody><tr><td>D21A001</td><td>ISW1-MA</td><td>0.646<sup>2</sup></td><td>0.595</td><td>40.537<sup>2</sup></td></tr><tr><td>D21A193</td><td>ISW1-MA</td><td>0.305<sup>3</sup></td><td>0.834</td><td>-51.371<sup>3</sup></td></tr><tr><td>E21A033</td><td>DAR-MA</td><td>0.453<sup>3</sup></td><td>0.740</td><td>47.035<sup>2</sup></td></tr><tr><td>E21A051</td><td>MDAR-MA</td><td>0.939<sup>1</sup></td><td>0.248</td><td>-10.006<sup>1</sup></td></tr><tr><td>E21A054</td><td>DAR-MA</td><td>0.897<sup>1</sup></td><td>0.322</td><td>-7.892<sup>1</sup></td></tr><tr><td>E21A056</td><td>DAR-MA</td><td>0.780<sup>1</sup></td><td>0.469</td><td>25.436<sup>2</sup></td></tr><tr><td>E21A066</td><td>MDAR-MA</td><td>0.749<sup>1</sup></td><td>0.501</td><td>-22.120<sup>1</sup></td></tr></tbody></table>

As can be seen in Table 12, for all target stations other than E21A033, E21A054, and E21A056, higher NSE values were obtained with MDAR compared with DAR. Especially for E21A066, the NSE value improved considerably with MDAR as compared with DAR. For D21A001, E21A051, E21A054, and E21A056, the NSE values decreased when MDAR was applied using the three most physically similar donor stations instead of the two most physically similar donor stations (see Table 12). In case of use of the three most physically similar donor stations, the third most physically similar donor stations for the target stations D21A001, E21A051, and E21A054 were determined as D21A193, E21A066, and E21A033, respectively. All NSE values obtained for these target stations using the DAR method utilizing their third most physically similar donor stations were negative. The drainage areas of these target stations and their third most physically similar donor stations are very close to each other. Therefore, the weight of their third most physically similar donor stations for streamflow estimation of these target stations is significantly larger. As a result, the NSE values obtained for the target stations D21A001, E21A051, and E21A054 using MDAR with three most physically similar donor stations are predominantly influenced by their third most physically similar donor stations. On the other hand, for E21A056, the reason is slightly different from the others. The NSE value obtained for E21A056 using the DAR method utilizing its third most physically similar donor station E21A033 was too low (i.e., -10.029). Although the contribution of the donor station E21A033 is not much more than the other two, this leads to poor estimation performance for E21A056 when MDAR was applied using the three most physically similar donor stations.

As can be seen in Table 13, in case of the use of the two most physically similar donor stations, the best performance results were obtained with ISW1 for 4 out of 7 target stations (i.e., D21A001, D21A193, E21A051, and E21A066). On the other hand, in case of the use of the three most physically similar donor stations, the best performance results were obtained with ISW3 for all target stations other than D21A193 and E21A066. As the power parameter increased from 1 to 3, the NSE values mostly improved when the three most physically similar donor stations were used, whereas the NSE values mostly decreased when the two most physically similar donor stations were used. Moreover,

<footnote n="1">¹ Very good, ² good, ³ satisfactory, and ⁴ unsatisfactory.</footnote>

<footnote n="1">¹ Very good, ² good, ³ satisfactory, and ⁴ unsatisfactory.</footnote>
the NSE values mostly improved when the two most physically similar donor stations were used instead of the three most physically similar donor stations.

As can be seen in Table 14, for all target stations other than for E21A033, E21A054, and E21A056, the MDAR and the ISW methods resulted in higher NSE compared to the DAR method. For 4 out of 7 target stations, the results can be rated as “very good” for the NSE according to the performance ratings in Table 5. For 5 out of 7 target stations, the RSR values were considered satisfactory (i.e., less than 0.7) according to the performance ratings recommended by Moriasi et al. [47]. The negative PBIAS values for D21A193, E21A051, E21A054, and E21A066 demonstrate that the method overestimated daily streamflow, while positive PBIAS values for D21A001, D21A033, and E21A056 demonstrate underestimation. For all target stations other than D21A193, the statistical methods with MA tend to achieve slightly higher NSE values. However, the PBIAS values of the target stations did not change when the statistical methods with MA were used.

For the target station E21A051 as the example, the estimated streamflow values from the statistical methods without the MA were compared to observed (original) streamflow values in the hydrograph and scatter plots in Figure 8. The remarkably better agreement between observed and estimated streamflow values by three statistical methods was obtained for E21A051 compared to the other stations in the Upper Euphrates Basin. Both MDAR and ISW1 gave a coefficient of determination (R²) of 0.94, which was higher than the R² values of 0.89 by using DAR. The NSE values for these methods ranged from 0.883 to 0.932, and the best NSE value was achieved by MDAR. The best NSE performance for E21A051 was obtained using MDAR with the two most physically similar donor stations.

<figure><img src="image_74.png" /><figcaption>**Figure 8.** Comparison of daily streamflow estimated using DAR (top row), MDAR (mid row), and ISW1 (bottom row) with observed (original) daily streamflow at station E21A051 in the Upper Euphrates Basin.</figcaption></figure>

On the other hand, for the target station E21A051 as the example, the estimated streamflow values from the statistical methods with the MA were compared to observed-MA (smoothed) streamflow values in the hydrograph and scatter plots in Figure 9. The statistical methods with the MA performed slightly better for E21A051.
<figure><img src="image_75.png" /><figcaption>**Figure 9.** Comparison of daily streamflow estimated using DAR-MA (top row), MDAR-MA (mid row), and ISW1-MA (bottom row) with observed-MA (smoothed) daily streamflow at station E21A051 in the Upper Euphrates Basin.</figcaption></figure>

**5. Conclusions**

This study provided a comparative evaluation of three statistical methods, DAR, MDAR, and ISW, which estimate daily streamflow at ungauged basins. These statistical methods were applied to two study basins: The Middle and Upper Euphrates basins in Turkey. DAR was implemented with the most physically similar donor station determined using the similarity coefficient. On the other hand, the two and the three most physically similar donor stations were used with both MDAR and ISW. By using three different power parameters (1, 2, and 3) in ISW, the effect of the selection of different power parameters on the accuracy of the daily streamflow estimation was tested. In addition, this study investigated the effects of the statistical methods using the smoothed data by the MA on the accuracy and reliability of daily streamflow estimation. Three statistical methods using the smoothed data by the MA, named DAR-MA, MDAR-MA, and ISW-MA, were proposed. The performance of each statistical method was evaluated in terms of the NSE, RSR, PBIAS, and $R^2$ between the observed and estimated daily streamflow. When the estimated daily streamflow values at the target station were obtained from the statistical methods using the observed (original) daily streamflow values at the donor station(s), they were compared to the observed (original) daily streamflow values at the target station. On the other hand, when the estimated daily streamflow values at the target station were obtained from the statistical methods using the observed-MA (smoothed) daily streamflow values at the donor station(s), they were compared to the observed-MA (smoothed) daily streamflow values at the target station. These two approaches were presented to estimate the daily streamflow values with and without MA. It is believed that the results will help decision makers choose the best one for their objectives.

In the Middle Euphrates Basin, the DAR method resulted in negative NSE values, indicating unsatisfactory performance for 3 out of 8 target stations when the most physically similar donor station was used. These negative NSE values obtained with DAR improved considerably when the three most
physically similar donor stations were used with MDAR. Higher NSE values were mostly obtained from both MDAR and ISW used with the three most physically similar donor stations instead of the two most physically similar donor stations. ISW with a power parameter of 1 (i.e., ISW1) mostly outperformed compared to ISW2 and ISW3, when both the two and the three most physically similar donor stations were used. The results obtained for 8 target stations in the Middle Euphrates Basin indicated that ISW for 4 stations, MDAR for 3 stations, and DAR for 1 station performed best in estimating daily streamflow. For all but one target station, the NSE values obtained were greater than 0.6, indicating good or very good performance. For all target stations, the performance of the best statistical method for each target station slightly improved when the smoothed data by the MA was used.

In the Upper Euphrates Basin, for one target station, the NSE value improved from a negative to over 0.7 when MDAR was applied instead of DAR. Higher NSE values were mostly obtained from both MDAR and ISW used with the two most physically similar donor stations instead of the three most physically similar donor stations. ISW1 used with the two most physically similar donor stations and ISW3 used with the three most physically similar donor stations gave better performance than the others. The results obtained for 7 target stations in the Upper Euphrates Basin indicated that DAR for 3 stations, MDAR for 2 stations, and ISW for 2 stations performed best in estimating daily streamflow. For all but two target stations, the NSE values obtained were greater than 0.6, indicating good or very good performance. For 6 out of 7 target stations, the performance of the best statistical method for each target station slightly improved when the smoothed data by the MA was used.

The overall results suggest that, besides the statistical method selection, the selection of appropriate donor stations is an important step to achieve better streamflow estimates at target stations. Also important is that increasing the number of donor stations can also improve or decrease estimation performance. Besides, the estimation performance of the statistical methods can vary from basin to basin. Moreover, data preprocessing can have a positive effect on the estimation performance of statistical methods.

Finally, since obtaining reliable and accurate streamflow estimations is very important in water resource studies, the statistical methods used in the study can be easily applied for decision making and design in many water resources projects that have difficulty in obtaining data.

**Author Contributions:** Conceptualization, M.U.Y. and B.O.; methodology, M.U.Y. and B.O.; software, M.U.Y.; formal analysis, M.U.Y. and B.O.; investigation, M.U.Y. and B.O.; writing—original draft preparation, M.U.Y.; writing—review and editing, M.U.Y. and B.O.; visualization, M.U.Y. and B.O.; supervision, B.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by Research Fund of Istanbul Technical University. Project Number: 40717.

**Acknowledgments:** This study was supported by the Research Fund of Istanbul Technical University under the project "Improvement of Streamflow Estimation in Ungauged Basins" (Project number: 40717). The authors thank the General Directorate of State Hydraulic Works (DSI) and the Turkish State Meteorological Service (DMI) for providing data.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Sivapalan, M.; Takeuchi, K.; Franks, S.W.; Gupta, V.K.; Karambiri, H.; Lakshmi, V.; Liang, X.; McDonnell, J.J.; Mendiondo, E.M.; O'Connell, P.E.; et al. IAHS decade on predictions in ungauged basins (PUB), 2003–2012: Shaping an exciting future for the hydrological sciences. *Hydrol. Sci. J.* **2003**, *48*, 857–880. [CrossRef]

2. Ergen, K.; Kentel, E. An integrated map correlation method and multiple-source sites drainage-area ratio method for estimating streamflows at ungauged catchments: A case study of the Western Black Sea Region, Turkey. *J. Environ. Manage.* **2016**, *166*, 309–320. [CrossRef] [PubMed]

3. Yilmaz, M.U.; Onoz, B. Evaluation of statistical methods for estimating missing daily streamflow data. *Tek. Derg.* **2019**, *30*, 9597–9620. [CrossRef]

4. Besaw, L.E.; Rizzo, D.M.; Bierman, P.R.; Hackett, W.R. Advances in ungauged streamflow prediction using artificial neural networks. *J. Hydrol.* **2010**, *386*, 27–37. [CrossRef]
5. Huo, Z.; Feng, S.; Kang, S.; Huang, G.; Wang, F.; Guo, P. Integrated neural networks for monthly river flow estimation in arid inland basin of Northwest China. *J. Hydrol.* **2012**, *420*, 159–170. [CrossRef]

6. Noori, N.; Kalin, L. Coupling SWAT and ANN models for enhanced daily streamflow prediction. *J. Hydrol.* **2016**, *533*, 141–151. [CrossRef]

7. Chang, F.J.; Chen, Y.C. A counterpropagation fuzzy-neural network modeling approach to real time stream-flow prediction. *J. Hydrol.* **2001**, *245*, 153–164. [CrossRef]

8. Ozger, M. Comparison of fuzzy inference systems for streamflow prediction. *Hydrol. Sci. J.* **2009**, *54*, 261–273. [CrossRef]

9. Toprak, Z.F.; Eris, E.; Agiralioglu, N.; Cigizoglu, H.K.; Yilmaz, L.; Aksoy, H.; Coskun, H.G.; Andic, G.; Alganci, U. Modeling monthly mean flow in a poorly gauged basin by fuzzy logic. *CLEAN-Soil Air Water* **2009**, *37*, 555–564. [CrossRef]

10. Khu, S.T.; Liong, S.Y.; Babovic, V.; Madsen, H.; Muttil, N. Genetic programming and its application in real-time runoff forecasting. *J. Am. Water Resour. Assoc.* **2001**, *37*, 439–451. [CrossRef]

11. Maity, R.; Kashid, S.S. Hydroclimatological approach for monthly streamflow prediction using genetic programming ISH. *J. Hydraul. Eng.* **2009**, *15*, 89–107. [CrossRef]

12. Mehr, A.D.; Kahya, E.; Olyaie, E. Streamflow prediction using linear genetic programming in comparison with a neuro-wavelet technique. *J. Hydrol.* **2013**, *505*, 240–249. [CrossRef]

13. Lin, J.Y.; Cheng, C.T.; Chau, K.W. Using support vector machines for long-term discharge prediction. *Hydrol. Sci. J.* **2006**, *51*, 599–612. [CrossRef]

14. Solomatine, D.P.; Maskey, M.; Shrestha, D.L. Instance-based learning compared to other data-driven methods in hydrological forecasting. *Hydrol. Process.* **2008**, *22*, 275–287. [CrossRef]

15. Yilmaz, A.; Muttil, N. Runoff estimation by machine learning methods and application to Euphrates Basin in Turkey. *J. Hydrol. Eng.* **2014**, *19*, 1015–1025. [CrossRef]

16. Wu, C.L.; Chau, K.W.; Li, Y.S. Methods to improve neural network performance in daily flows prediction. *J. Hydrol.* **2009**, *372*, 80–93. [CrossRef]

17. Di, C.; Yang, X.; Wang, X. A Four-stage hybrid model for hydrological time series forecasting. *PLoS ONE* **2014**, *9*, e104663. [CrossRef]

18. Mehr, A.D.; Kahya, E. A Pareto-optimal moving average multigene genetic programming model for daily streamflow prediction. *J. Hydrol.* **2017**, *549*, 603–615. [CrossRef]

19. Zhou, J.; Peng, T.; Zhang, C.; Sun, N. Data pre-analysis and ensemble of various artificial neural networks for monthly streamflow forecasting. *Water* **2018**, *10*, 628. [CrossRef]

20. Wu, C.L.; Chau, K.W. Prediction of rainfall time series using modular soft computing methods. *Eng. Appl. Artif. Intell.* **2013**, *26*, 997–1007. [CrossRef]

21. Waseem, M.; Ajmal, M.; Kim, U.; Kim, T.W. Development and evaluation of an extended inverse distance weighting method for streamflow estimation at an ungauged site. *Hydrol. Res.* **2016**, *47*, 333–343. [CrossRef]

22. Razavi, T.; Coulibaly, P. Improving streamflow estimation in ungauged basins using a multi-modelling approach. *Hydrol. Sci. J.* **2016**, *61*, 2668–2679. [CrossRef]

23. Huang, W.C.; Yang, F.T. Streamflow estimation using kriging. *Water Resour. Res.* **1998**, *34*, 1599–1608. [CrossRef]

24. Farmer, W.H. Ordinary kriging as a tool to estimate historical daily streamflow records. *Hydrol. Earth Syst. Sci.* **2016**, *20*, 2721–2735. [CrossRef]

25. Tencaliec, P.; Favre, A.C.; Prieur, C.; Mathevet, T. Reconstruction of missing daily streamflow data using dynamic regression models. *Water Resour. Res.* **2015**, *51*, 9447–9463. [CrossRef]

26. Masselot, P.; Dabo-Niang, S.; Chebana, F.; Ouarda, T.B. Streamflow forecasting using functional regression. *J. Hydrol.* **2016**, *538*, 754–766. [CrossRef]

27. Swain, J.B.; Patra, K.C. Streamflow estimation ungauged catchments using regional flow duration curve: comparative study. *J. Hydrol. Eng.* **2017**, *22*, 04017010. [CrossRef]

28. Burgan, H.I.; Aksoy, H. Annual flow duration curve model for ungauged basins. *Hydrol. Res.* **2018**, *49*, 1684–1695. [CrossRef]

29. Farmer, W.H.; Vogel, R. Performance-weighted methods for estimating monthly streamflow at ungauged sites. *J. Hydrol.* **2013**, *477*, 240–250. [CrossRef]
30. Yilmaz, M.U. Performance-Weighted Methods for Estimating Monthly Streamflow: An Application for Middle Part of Euphrates Basin. Master’s Thesis, Istanbul Technical University, Institute of Science and Technology, Istanbul, Turkey, 2014.

31. Merz, R.; Blöschl, G. Regionalisation of catchment model parameters. *J. Hydrol.* **2004**, *287*, 95–123. [CrossRef]

32. Shu, C.; Ouarda, T.B.M.J. Improved methods for daily streamflow estimates at ungauged sites. *Water Resour. Res.* **2012**, *48*. [CrossRef]

33. HRachowitz, M.; Savenije, H.H.G.; Blöschl, G.; McDonnell, J.J.; Sivapalan, M.; Pomeroy, J.W.; Arheimer, B.; Blume, T.; Clark, M.P.; Ehret, U.; et al. A decade of Predictions in Ungauged Basins (PUB)—A review. *Hydrol. Sci. J.* **2013**, *58*, 1198–1255. [CrossRef]

34. Alipour, M.H.; Kibler, K. A framework for streamflow prediction in the world's most severely data-limited regions: test of applicability and performance in a poorly-gauged region of China. *J. Hydrol.* **2018**, *557*, 41–54. [CrossRef]

35. Alipour, M.H.; Kibler, K.M. Streamflow prediction under extreme data scarcity: a step toward hydrologic process understanding within severely data-limited regions. *Hydrol. Sci. J.* **2019**, *64*, 1038–1055. [CrossRef]

36. Patil, S.; Stieglitz, M. Controls on hydrologic similarity: role of nearby gauged catchments for prediction at an ungauged catchment. *Hydrol. Earth Syst. Sci.* **2012**, *16*, 551–562. [CrossRef]

37. Heng, S.; Suetsugi, T. Comparison of regionalization approaches in parameterizing sediment rating curve in ungauged catchments for subsequent instantaneous sediment yield prediction. *J. Hydrol.* **2014**, *512*, 240–253. [CrossRef]

38. Yang, X.; Magnusson, J.; Rizzi, J.; Xu, C.Y. Runoff prediction in ungauged catchments in Norway: comparison of regionalization approaches. *Hydrol. Res.* **2017**, *49*, 487–505. [CrossRef]

39. Emerson, D.G.; Vecchia, A.V.; Dahl, A.L. Evaluation of Drainage-Area Ratio Method Used to Estimate Streamflow for the Red River of the North Basin, North Dakota and Minnesota; U.S. Geological Survey Scientific Investigations Report 2005-5017; U.S. Geological Survey: Reston, VA, USA, 2005.

40. Asquith, W.H.; Roussel, M.C.; Vrabel, J. Statewide Analysis of the Drainage-Area Ratio Method for 34 Streamflow Percentile Ranges in Texas; U.S. Geological Survey Scientific Investigations Report 2006-5286; U.S. Geological Survey: Reston, VA, USA, 2006.

41. Oudin, L.; Andréassian, V.; Perrin, C.; Michel, C.; Le Moine, N. Spatial proximity, physical similarity, regression and ungaged catchments: A comparison of regionalization approaches based on 913 French catchments. *Water Resour. Res.* **2008**, *44*. [CrossRef]

42. Zelelew, M.B.; Alfredsen, K. Transferability of hydrological model parameter spaces in the estimation of runoff in ungauged catchments. *Hydrol. Sci. J.* **2014**, *59*, 1470–1490. [CrossRef]

43. Arsenault, R.; Brissette, F.P. Continuous streamflow prediction in ungauged basins: the effects of equifinality and parameter set selection on uncertainty in regionalization approaches. *Water Resour. Res.* **2014**, *50*, 6135–6153. [CrossRef]

44. He, Y.; Bárdossy, A.; Zehe, E. A review of regionalisation for continuous streamflow simulation. *Hydrol. Earth Syst. Sci.* **2011**, *15*, 3539–3553. [CrossRef]

45. Burn, D.H.; Boorman, D.B. Estimation of hydrological parameters at ungauged catchments. *J. Hydrol.* **1993**, *143*, 429–454. [CrossRef]

46. Nash, J.E.; Sutcliffe, J.V. River flow forecasting through conceptual models part I—A discussion of principles. *J. Hydrol.* **1970**, *10*, 282–290. [CrossRef]

47. Moriasi, D.N.; Arnold, J.G.; Van Liew, M.W.; Bingner, R.L.; Harmel, R.D.; Veith, T.L. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. *Trans. ASABE* **2007**, *50*, 885–900. [CrossRef]

48. Christiansen, D.E.; Haj, A.E.; Risley, J.C. Simulation of Daily Streamflow for 12 River Basins in Western Iowa Using the Precipitation-Runoff Modeling System; U.S. Geological Survey Scientific Investigations Report 2017-5091; U.S. Geological Survey: Reston, VA, USA, 2017.

49. Kalin, L.; Isik, S.; Schoonover, J.E.; Lockaby, B.G. Predicting water quality in unmonitored watersheds using artificial neural networks. *JEQ* **2010**, *39*, 1429–1440. [CrossRef]

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

© 2020 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 (http://creativecommons.org/licenses/by/4.0/).
Article

# UAV and LiDAR Data in the Service of Bank Gully Erosion Measurement in Rambla de Algeciras Lakeshore

Radouane Hout<sup>1,*</sup>, Véronique Maleval<sup>1</sup>, Gil Mahe<sup>2</sup>, Eric Rouvellac<sup>1</sup>, Rémi Crouzevialle<sup>1</sup> and Fabien Cerbelaud<sup>1</sup>

<sup>1</sup> Department of Geography, Faculté des Lettres et Sciences Humaines, Limoges University, 39E rue Camille Guérin, 87036 Limoges CEDEX, France; veronique.maleval@unilim.fr (V.M.); eric.rouvellac@unilim.fr (E.R.); remi.crouzevialle@unilim.fr (R.C.); fabien.cerbelaud@unilim.fr (F.C.)

<sup>2</sup> IRD, UMR HSM IRD, University Montpellier, 34090 Montpellier, France; gil.mahe@ird.fr

* Correspondence: radouane.hout@unilim.fr; Tel.: +33-7-69-21-47-85

Received: 30 July 2020; Accepted: 24 September 2020; Published: 1 October 2020

**Abstract:** The Rambla de Algeciras lake in Murcia is a reservoir for drinking water and contributes to the reduction of flooding. With a semi-arid climate and a very friable nature of the geological formations at the lakeshore level, the emergence and development of bank gullies is favored and poses a problem of silting of the dam. A study was conducted on these lakeshores to estimate the sediment input from the bank gullies. In 2018, three gullies of different types were the subject of three UAV photography missions to model in high resolution their low topographic change, using the SfM-MVS photogrammetry method. The combination of two configurations of nadir and oblique photography allowed us to obtain a complete high-resolution modeling of complex bank gullies with overhangs, as it was the case in site 3. To study annual lakeshore variability and sediment dynamics we used LiDAR data from the PNOA project taken in 2009 and 2016. For a better error analysis of UAV photogrammetry data we compared spatially variable and uniform uncertainty models, while taking into account the different sources of error. For LiDAR data, on the other hand, we used a spatially uniform error model. Depending on the geomorphology of the gullies and the configuration of the data capture, we chose the most appropriate method to detect geomorphological changes on the surfaces of the bank gullies. At site 3 the gully topography is complex, so we performed a 3D distance calculation between point clouds using the M3C2 algorithm to estimate the sediment budget. On sites 1 and 2 we used the DoD technique to estimate the sediment budget as it was the case for the LiDAR data. The results of the LiDAR and UAV data reveal significant lakeshore erosion activity by bank gullies since the annual inflow from the banks is estimated at 39 T/ha/year.

**Keywords:** Rambla de Algeciras; semi-arid; lake; lakeshores; silting; bank gullies; UAV; LiDAR; DoD; M3C2

## 1. Introduction

Gully erosion is a very common form of water erosion, this geomorphological process is triggered after the concentration of surface or groundwater in narrow areas. This drags soil particles and rock fragments from these narrow flow paths to considerable depths. The steep slopes, the low density of plant cover, the overexploitation of arable lands, overgrazing as well as the torrential flow regime due to irregular rainfall, are all factors favorable to the initiation and even the acceleration of the gully erosion phenomenon [1–3].
This process of global erosion occurs in urban, forest and even rural areas [4]. That is why we find the term “gully” depicted in several languages, it is called, “Ravine” in French, “wadi” in Arabic, “'cárcava” in Spain, “lavaka” in Madagascar, “donga” in South Africa, “voçoroca” in Brazil and “barranco” in Argentina [4].

Globally, gullying accounts for 10% up 94% of the total sediment yield produced by water erosion [5], with an average contribution of 50% to 80% in arid areas [6]. In Western European countries, the same process generates between 30 and 80% [5] of the sediment.

Bank gullies are themselves a very particular form of gullies which are formed after a height drop from terraces or river banks, the lake banks, and which develop by headward retreat in erodible hillslopes. Bank gullies are themselves a very particular form of gullies which are formed after a height drop from terraces or river banks, the lake banks, and which develop by headward retreat in erodible hillslopes.

Furthermore, Vanmaercke et al (2016) [7] have demonstrated that gullying and bank gullies are the main sources of sediment loss. It is also important to note that they also have off-site effects [8] such as degradation of water quality, siltation of dams, pollution, degradation of agricultural land, reduction of soil moisture, hydrological dysfunction, and even fertility transfers [9]. Vanmaercke et al (2016) [7] confirmed that the recession rate of bank gullies heads varies between 0.01 m/year and 135.2 m/year with an average of 5 m/year, and a median of 0.89 m/year in Spain. In most landscapes that are subject to different climatic conditions and different anthropogenic activities, we can observe the presence of various forms of gully erosion, namely ephemeral gullies, permanent or classic gullies and bank gullies [5].

In the literature, researchers address the question of gullying according to different approaches:

* Approach based on the quantitative estimation of the sediment balance, using direct measurements in the field, namely, measuring tape, total station (TS), DGPS (differential GPS), drone (UAV), LiDAR, etc. [10];

* Indirect approach based on the modeling of the gully's sub-processes, by numerical or statistical gully forecasting models such as the SCS (Soil Conservation Service) model [11], the Thompson models [12,13], AnnAGNPS [14], EGEM (Ephemeral Gully Erosion Model) [15], GULTEM (model to Predict Gully Thermoerosion and Erosion) [16], or CHILD (Channel-Hillslope Integrated Landscape Development) [17].

The work of Evans [18] has shown that the activity of bank gullies lead to the evolution of hillslope-channel coupling related to flow erosion. In Morocco, the work of Maleval [19] demonstrated the importance of sedimentary inputs from the bank gullies in the silting of the SMBA dam. Monitoring and mapping of gullies for lakeshores remain one of the very important objectives in the management of drinking water reservoirs. Very high resolution topographic acquisitions have considerably improved our ability to better understand the functioning and evolution of geomorphological changes in the shoreline. The emergence of new drone photography techniques in recent years has greatly facilitated digital terrain modeling and geomorphological analysis. This ease is mainly due to the appearance of civil drones (Unmanned Airborne Vehicles, UAVs) and to SfM-MVS 3D photogrammetry algorithms, but also to the development of the computational capacity of computers. This method is increasingly being used for modeling and quantifying gullying in several regions of the world and in various climatic contexts, with work carried out in Morocco [20], Spain [21], China [22], Italy [23], Ethiopia [24], the United States [25], Australia [26], etc.
In this paper, the gullies for lakeshores of the Rambla de Algaciras dam have been mapped using UAV and LiDAR. Our objective was to understand the influence of gullies for lakeshores in the silting of artificial dams using 4D models and photo-interpretation of aerial photos. A detailed methodology for the collection, processing and validation of very high resolution 4D models has been developed to quantify sediment budgets and geomorphological changes of lake bank gullies in the short and medium term. The topographic complexity of certain bank gullies was taken into account in this work through a combination of oblique and vertical acquisition of aerial photos.

In this article we also measured the accuracy of digital terrain models taken by UAVs with the k-fold cross-validation method, to estimate the sediment budget on bank gullies using the DoD method and the M3C2-PM algorithm. The conclusions of this study will guide dam managers to implement several strategies to take into account the contributions of bank gullies which are generally underestimated in the hydro-sediment balances.

## 2. Methodology

### 2.1. Study Site

In the south of Spain, in Alhama de Murcia, the Rambla de Algeciras reservoir was built in 1995 on the Rambla de Algeciras river. It has a torrential character, it rises at the gully of Valdelaparra and flows into the Guadalentín river. The dam is 80 m high and the area of the lake is 235.50 ha and 42.13 $hm^3$ (Figure 1). The dam was built under the project "Plan General de Defensa contra Avenidas de la Cuenca del Segura", intended to develop and regularize water resources in the Segura watershed in south-eastern Spain. The Rambla de Algeciras dam is located in the catchment area by the same name, southwest of the city of Murcia at the northern limit of the depression of the Guadalentín river (left bank). Surrounded by las sierras del Cura, de Espuña and La Muela. This watershed covers an area of 44.91 $km^2$, and 15 km, with an altitude that varies between 200 and 1320 m above sea level. It undergoes a high rate of erosion reaching 100 t/ha/year (Global Atlas of the Murcia region 2002), linked to the torrential nature of precipitation and the detrital nature of the soil composed of marls and polygenic conglomerates.

Locally, Lake Rambla de Algeciras has a semi-arid climate, with short episodes of intense precipitation. Annual rainfall 393 mm, 60% of which is concentrated over 3 to 4 days/year.

The lakeshore and substrate of the Rambla de Algeciras dam are made up of different friable formations of the Tortonian period, on the banks of the lake we find blue loams interspersed with whitish calcareous loams and gypsiferous loams in discordance with calcarenite and molasses formations, on the right bank of the dam we note the presence of alternating soft and consolidated layers of polygenic conglomerates in a marly sand matrix.

The choice of gullies aimed to represent as faithfully as possible the erosion activity on the lakeshore, in a semi-arid climatic context. This choice takes into account several parameters such as the shape, the exposure, the head of the gullies, the density of the plant cover, and the lithology. Based on the diachronic photo-interpretation of aerial photos (PNOA), three gullies models were selected (Figure 1).
<figure><img src="image_77.png" /><figcaption>Figure 1. Location of the study area.</figcaption></figure>

2.1.1. Site 1

The bank gully (B.G-1) (Figure 2C), extends over approximately 130 m. It is V-shaped in the upper and middle part, then takes a U-shape downstream from the gully, where the slopes are steep and exposed to the East (low grass cover) or to the West (bare soil). Upstream of the gully (head) is covered by tree layers, planted as part of the restoration project (Proyecto de correccion y repoblacion forestal de las cuencas de las Ramblas de Belen, Librillas y Algeciras. Provincia de Murcia). The North-south oriented gullies have been reforested and are difficult to access. Thus we studied the B.G-1 which is located 250 m from the reservoir. Access to the bank gully of the lakeshore protected by the forest cover and is very difficult. It is for this reason that we chose site B.G-1 which is located at 250 m from the lake.

2.1.2. Site 2

The second bank gully (B.G-2)(B.G-1) (Figure 2D) is U-shaped. It is oriented East-west with very steep slopes facing North (herbaceous vegetation) and South (bare soil). It is 110 m long, and receives tributary gullies.

2.1.3. Site 3

The third bank gully (B.G-3)(B.G-1) (Figure 2B) has a complex shape and is the result of a combination of tunnel gullying (pipig). The slopes of the gully (B.G-3) covered by bare soil. According to photo-interpretation, the B.G-3 gully began to develop after the lake level receded. It's likely to turn into badlands.
<figure><img src="image_78.png" /><figcaption>**Figure 2.** Photos of the lakeshore of Rambla de Algeciras. (A) Lakeshore. (B) Third bank gully. (C) Bank gully (B.G.-1). (D) Second bank gully.</figcaption></figure>

## 2.2. Precipitation Analysis

### 2.2.1. Rainfall Extremes

Rare and very rare extreme events were identified using the prec90p index, of the STARDEX project (Statistical and Regional Dynamical downscaling of Extremes for European regions). The prec90p index is used to determine heavy rain events and intense rain when the cumulative daily precipitation values are greater than or equal to 90 and 95 percentiles, respectively. Between 2009 and 2016 there were two heavy showers (Figure 3).

* 7 June 2011, with cumulative precipitation of 81 mm in 2 h and a high intensity of precipitation which exceeds 40 mm/h.

* 28 September 2012, with cumulative precipitation was 84 mm produced in 4 h and an intensity of 21 mm/h.

These two rain events are considered as extreme events and automatically play an important role in erosion. The modification of the precipitation regime with the decrease in the amount of precipitation and rainy days as well as the greater frequency of rare and very rare events during the eight years of topographic surveys can be explained by the effect of climate change in the Mediterranean area [27].
<figure><img src="image_79.png" /><figcaption>Figure 3. Rain event between 2009 and 2019 in the station of Rambla de Algeciras.</figcaption></figure>

### 2.2.2. Rainfall Erosivity Factor “R-Factor”

The R-factor was estimated using equation (Equation (1)) (zone-2) proposed by ICONA for the region of Murcia. The monthly analysis of the factor R traces the evolution of the aggressiveness of the rains on the gullies between the study period of 2009 and 2016 as well as for the year 2018. Between 2009 and 2016 the devastating activity of the rains is concentrated annually in March and September. This activity is linked to the increase in intense rain events during this period. Throughout the period 2009 to 2016, the R-factor peaks are recorded in June 2011 and September 2012, which coincides with an anomaly in the kinetic energy explained by the recording of two extreme phenomena. The year 2018 is relatively wet therefore implying a natural increase in the factor R. From the graph (Figure 4) we can see that the average value of R, calculated for the wet season, is quite low compared to that of the dry season. During the wet season the values of (R) are moderate with an average of 8.8 MJ cm ha⁻¹ h⁻¹ year⁻¹, through contribution of the season which recorded an average of 17 MJ cm ha⁻¹ h⁻¹ year⁻¹. The lowest values are recorded between February and April. Conversely, during the dry season, the highest cumulative value of R is 102 MJ cm ha⁻¹ h⁻¹ year⁻¹, and it is associated with very intense precipitation in August and November. R: Climate aggressiveness (MJ cm ha⁻¹ h⁻¹ year⁻¹) PMEX: The maximum monthly value (mm) MR: Cumulative precipitation between October–May (mm) MV: Cumulative precipitation between June–September (mm):

$$R = e^{-1.235} \times \text{PMEX}^{1.297} \times \text{MR}^{-0.511} \times \text{MV}^{0.366} \times F_{24}^{0.414}; \quad (1)$$

where

R: Rainfall erosivity factor (MJ cm ha⁻¹ h⁻¹ year⁻¹)

PMEX: Maximum monthly pre-cipitation (mm)

MR: Total rainfall from October to May (mm)

MV: Total rainfall from June–September (mm)

F₂₄: Ratio of the square of the maximum annual rainfall in 24 hours (mm) to the sum of the maximum

$$F_{24} = \frac{(P_{24h,annual})^2}{\sum_{i=1}^{12} P_{24h,i}}.$$
<figure><img src="image_80.png" /><figcaption>**Figure 4.** R-factor.</figcaption></figure>

## 2.3. UAV Data Processing

### 2.3.1. UAV Data Collection

The aerial photography missions were carried out easily and at low cost by the quadrirotor drone DJI Inspire 1 v2 (Figure 5A), equipped with the metric and multispectral camera ZENMUSE X3 without distortion or aspect and with a focal distance of 20 mm (Figure 5D). The DJI Inspire1 V2 is sold "Ready-To-Fly", and is designed to have a maximum speed of 22 m/s (ATTI mode, no wind), and a maximum flight time of approximately 18 min. It can reach 2500 m altitude in flight. In addition, it is equipped with a satellite positioning system which uses GLONASS and GPS system.

<figure><img src="image_81.png" /><figcaption>**Figure 5.** (A) Drone DJI Inspire 1 v2. (B) GPS Trimble Geo 7x. (C) Ground control points. (D) Camera ZENMUSE X3.</figcaption></figure>
Flight mission planning was done in navigation mode via the Pix4Dcapture mobile application. After the definition of some flight parameters (flight plan, flight altitude, recovery rate, shooting angle, starting and arrival point of the UAV) the application automatically pilots the UAV. The acquisitions of aerial photos by drone were carried out according to two configurations, for the two gullies (B.G 1 and 2) the photos were taken nadir acquisitions. However, by geomorphological constraints of having gully (B.G-3) we combined the two modes of acquisition nadir and oblique at (30°) in order to model in 3D the overhanging areas of the gully (B.G-3). To ensure a successful photographic reconstruction we have chosen a longitudinal overlap of 80% and a lateral overlap of 60% to avoid gaps between the different photographs. A minimum of 10 ground control points (GCPs) (Figure 5C) (Table 1) is required on each site in the middle and at the stable edge, and the same targets were used in each drone mission (except for those which were covered by the sediments), in order to obtain optimal results. The surveys were carried out by the GPS Trimble Geo 7x (Figure 5B), which includes the following parameters: the ellipsoid GRS 1980, the datum: European Terrestrial Reference System 1989, coordinate system: ETRS89 / UTM zone 30N, with a horizontal accuracy of 1 cm + 1 ppm, and vertical of 1.5 cm + 2 ppm. We used the Stop&Go kinematic approach for the collection of GCPs points. We then made subsequent corrections using RINEX data from a reference station located 30 km from our study area.

<table><caption>**Table 1.** Characteristics of UAV photo missions.</caption><thead><tr><th>Bank Gullies</th><th>Missions</th><th>Local Time</th><th>Flight Height (m)</th><th>Cloud Cover (Octas)</th><th>Number of Images</th><th>Density (pts/m²)</th><th>Coverage (Ha)</th><th>Number of GCPs</th><th>Angle of Acquisition (°)</th></tr></thead><tbody><tr><td rowspan="3">B.G-1</td><td>22/01/18</td><td>11 h</td><td>27.6</td><td>02/08</td><td>436</td><td>517</td><td>2.95</td><td>7</td><td>90</td></tr><tr><td>04/06/18</td><td>15 h</td><td>36.4</td><td>01/08</td><td>304</td><td>408</td><td>2.34</td><td>7</td><td>90</td></tr><tr><td>18/12/18</td><td>16 h</td><td>27.5</td><td>03/08</td><td>422</td><td>524</td><td>2.7</td><td>7</td><td>90</td></tr><tr><td rowspan="3">B.G-2</td><td>22/01/18</td><td>15 h</td><td>31.3</td><td>02/08</td><td>646</td><td>425</td><td>4.3</td><td>13</td><td>90</td></tr><tr><td>04/06/18</td><td>11 h</td><td>29.4</td><td>01/08</td><td>465</td><td>455</td><td>1.79</td><td>13</td><td>90</td></tr><tr><td>18/12/18</td><td>15 h</td><td>30</td><td>03/08</td><td>384</td><td>397</td><td>2.02</td><td>10</td><td>90</td></tr><tr><td rowspan="3">B.G-3</td><td>22/01/18</td><td>17 h</td><td>37.4</td><td>02/08</td><td>741</td><td>394</td><td>2.28</td><td>16</td><td>90 + 30</td></tr><tr><td>04/06/18</td><td>10 h</td><td>28.9</td><td>01/08</td><td>520</td><td>433</td><td>1.28</td><td>14</td><td>90 + 30</td></tr><tr><td>18/12/18</td><td>10 h</td><td>33.3</td><td>03/08</td><td>602</td><td>328</td><td>2.5</td><td>11</td><td>90 + 30</td></tr></tbody></table>

### 2.3.2. Processing

The aerial photo processing is based on the SfM-MVS photogrammetry technique. We have chosen the Agisoft PhotoScan commercial software to process all the photos taken by the drone, eventually generating dense point clouds and orthomosaics. This software generates exportable datasets and thus has the most comprehensive and easy to use workflow compared to other commercial and open source photogrammetry software's. The following steps on PhotoScan consisted of establishing dense point clouds which we will then export in ASCII or PTS "XYZ RGB" format (Figure 6).
<figure><img src="image_82.png" /><figcaption>Figure 6. A schematic of a workflow that was applied to UAV data.</figcaption></figure>

### 2.3.3. Post-Processing

The SfM-MVS technique is recognized as the fastest and inexpensive technique for 3D modeling. Post-processing was applied to the raw point clouds obtained by the SfM-MVS technique before generating the digital terrain model (Figure 6). The post-processing step was carried out by three open source external software: CloudCompare software (http://www.cloudcompare.org/), ToPCAT (https://code.google.com/archive/p/point-cloud-tools/downloads), and SAGA GIS software

#### Filtering and Coregistration of Data

CloudCompare open source software is designed for processing and comparing 3D point clouds. The dense 3D point clouds obtained were imported in ASCII XYZRGB formats into the CloudCompare software. We then filtered the point clouds using the SOR filter. Once the point clouds were filtered, we applied the ICP coregistration algorithm to match the homologous points of the multi-temporal point clouds, so that each dense point cloud would refer to the same area [28–32].

#### The Decimation and Interpolation of Point Clouds

Generating a digital terrain model from dense point clouds directly derived from MVS algorithms often requires large computing resources. The interpolation of the points for this kind of processing drains the machines in terms of memory, hence the interest for the decimation of the clouds of dense points which has the potential to reduce the number of points while preserving the topographic complexity of the soil surface. Persistent in the point cloud, we used the ToPCAT geospatial toolbox (Topographic Point Cloud Analysis Toolkit). ToPCAT is an intelligent analysis and decimation algorithm that breaks down the very dense point cloud into a set of non-overlapping meshes, suitable for detecting morphological changes. It then calculates the minimum elevation ($Z_{min}$) of each mesh. ($Z_{min}$) is considered to be the elevation of bare ground points [33–36]. The ($Z_{min}$) elevation obtained
was then converted to TIN (Triangular Irregular Networks) in order to generate the digital terrain model under the SAGA GIS software.

## 2.4. LiDAR-PNOA Data Processing

In the Murcia region, two ALS airborne LiDAR missions were conducted, on 24 November 2009 and between 20 and 25 August 2016 (Figure 7). During the first mission, data were collected using a Leica ALS50-II sensor, and for the second mission, data were collected using a Leica ALS60-II sensor.

<figure><img src="image_83.png" /><figcaption>Figure 7. LiDAR flight turnings.</figcaption></figure>

The LAZ files (binary file compression format LAS version 1.1) that we used come from the ETRS89 UTM Zone 30 geodetic reference system for X, Y and Z coordinates in the EGM08-REDNAP system. Before the classification of the LiDAR point cloud, we prepared the LiDAR raw data in the open source software LASTools. We first converted the LAZ files into LAS files using the las2las algorithm, then we applied a vertical and a horizontal alignment between the flight lines using the lasoverlap algorithm, and finally we transformed the point cloud into EGM08-REDNAP orthometric coordinates using the ascii file <EGM08-REDNAP.asc> and the lasheight algorithm.

Subsequently, we removed outliers and noise using the SOR filter in the CloudCompare software. After pre-processing the point cloud, we chose the MCC-LiDAR v2.15 (Multiscale Curvature Classification) application, an open source command line tool for processing LiDAR data. It classifies data points as ground or non-ground, using the MCC algorithm, to obtain a digital model of bare terrain. Due to the low density of PNOA LiDAR point clouds ranging from 0.5 to 1 point/m², we did not apply a mesh to the cloud.

Finally, the ground points were triangulated and then interpolated with the linear interpolation algorithm under the SAGA GIS software, in order to obtain two ground DTMs of 1 m resolution for 2009 and 2016.

### 2.4.1. Error Modeling

#### Cross Validation: K-fold Cross-Validation

It is likely that the roughness-based estimates are incorrect, For quantify spatially variable SfM-DEM uncertainty we used cross-validation, then the results of this technique will be integrated
from a FIS error model more representative of the different sources of uncertainty. This statistical technique of resampling without replacement consists in dividing the decimated point cloud into validation (20%) and training data (80%) (Figure 8). The clouds are divided repeatedly with the k-folds (k-fold cross validation) algorithm with a k = 5 (5 subsets of validation point cloud). The training points are interpolated to build a DEM, which compared the validation points in order to quantify the uncertainty of DEM-SfM. We chose the Scikit-Learn (v-0.22.0) free software machine learning library to apply k-fold cross validation on the point cloud decimated by the ToPCAT algorithm. 80% of the (Z<sub>min</sub>) points were randomly selected to generate a DEM with a resolution of 5 cm. Then we used the remaining 20% of (Z<sub>min</sub>) to compare them with the DEM. The average absolute difference is taken as an indicator of the altimetric uncertainty of the data. The cross-validation test was repeated five times to check the consistency of the results.

<figure><img src="image_84.png" /><figcaption>Figure 8. K-fold cross-validation model.</figcaption></figure>

Fuzzy Inference System (FIS)

Within the SfM-MVS data, the sources of error vary according to the type of site. In the case of bank gullies, the main sources of error are roughness, slope of the terrain, camera geometry, interpolation error and point cloud density. In this case we were interested in a more comprehensive error modeling that takes into account the different sources of uncertainty. To meet the modelling requirements we adopted a spatially variable FIS (Fuzzy Inference System) error model, which is more reliable model errors in the SfM-MVS-4D data [35,37].

We defined four variables as the error sources of the SfM-MVS data to model with the FIS, namely slope, point density, interpolation error, and 3D point cloud error. Slope is the parameter most representative of the topographic complexity of the terrain, and it is the easiest parameter to derive from the DTM or point clouds. The slope (degrees) was generated directly from the decimated 'Zmin' point clouds with a resolution of 5 cm.

The modeling of the bare ground roughness from the value of the standard deviation of the decimated dense point clouds σ<sub>Zmin</sub> can be considered to take into account the main source of error of the SfM-MVS-4D data. Terrain roughness affects the ability of the SfM-MVS algorithm to reconstruct the dense point cloud, since the error due to roughness depends on the position of the scan, which is never the same for multi-date missions, as well as the nature and occupancy of the terrain. We recall that vegetation is absent on our study gullies. This facilitates the calculation of roughness.

Then we also generated the interpolation error from a simple absolute difference between the altitude of the TIN node and the DTM-SfM Z<sub>DTM</sub> - Z<sub>NodesTIN</sub>. Finally, for the error variable of the 3D point clouds we used the results of cross-validation.

The matrix data of the different variables were reclassified into four classes ('low', 'medium', 'high' and 'extreme') (Table 2).
<table><caption>**Table 2.** A four input fuzzy inference system for elevation uncertainty.</caption><thead><tr><th rowspan="2">Rule</th><th colspan="4">Input</th><th>Output</th></tr><tr><th>Slope (%)</th><th>Point Density (Pts/m²)</th><th>3-D Point Quality</th><th>Interpolation Error (m)</th><th>Elevation Uncertainty (m)</th></tr></thead><tbody><tr><td>1</td><td>Low</td><td>Low</td><td>NA</td><td>NA</td><td>High</td></tr><tr><td>2</td><td>Low</td><td>Medium</td><td>NA</td><td>NA</td><td>Medium</td></tr><tr><td>3</td><td>Low</td><td>High</td><td>NA</td><td>NA</td><td>Low</td></tr><tr><td>4</td><td>Medium</td><td>Low</td><td>NA</td><td>NA</td><td>High</td></tr><tr><td>5</td><td>Medium</td><td>Medium</td><td>NA</td><td>NA</td><td>Medium</td></tr><tr><td>6</td><td>Medium</td><td>High</td><td>NA</td><td>NA</td><td>Low</td></tr><tr><td>7</td><td>High</td><td>Low</td><td>NA</td><td>NA</td><td>High</td></tr><tr><td>8</td><td>High</td><td>Medium</td><td>NA</td><td>NA</td><td>High</td></tr><tr><td>9</td><td>High</td><td>High</td><td>NA</td><td>NA</td><td>Medium</td></tr><tr><td>10</td><td>Extreme</td><td>NA</td><td>NA</td><td>NA</td><td>Extreme</td></tr><tr><td>11</td><td>NA</td><td>NA</td><td>High</td><td>High</td><td>Extreme</td></tr></tbody></table>

## 2.5. Geomorphic Change Detection

An error model makes it possible to associate the control measurements made in the field with an accuracy estimate (example standard deviation in X, Y and Z). When applying the SfM-MVS method to create a digital terrain model, the errors generated are related to several factors such as the position and precision of the control points, camera distortions, resolution, number of photos, processing software, topographic complexity, interpolation methods, surface humidity, vegetation, etc., [38–40]. Several sources of error also affect the quality of DTMs for LiDAR data. These sources of error include terrain complexity, vegetation, weather, variation in point density, data filtering process and interpolation technique [41].

### 2.5.1. DEM of Difference (DoD)

This approach measures the vertical distance between two DEMs, through a simple subtraction between the pixels of these two DEMs according to the Equation (2). But before calculating the difference between the DEMs, a threshold value must be set for the minLOD (Equation (2)) [42,43]. This threshold means that subtraction values lower than minLOD are revealed as being uncertainties of the two numerical models. We can thus assign a value for the minLOD with confidence intervals that we can set. The user must define this threshold (statistic) according to the probability that he wishes to have so that the detection of changes can be real (Equation (2)) [44]. This fast method allows to determine the uncertainties related to the roughness of the terrain, the quality of the point clouds and the acquisition of the point clouds. However, great uncertainty in the quantification of sediment budgets occurs when the DEMs are not precise (since the minLOD threshold will be higher).

$$minLoD = \sqrt{\varepsilon_{2i}^2 + \varepsilon_{1i}^2} \quad \text{or} \quad minLoD = t\sqrt{\varepsilon_{2i}^2 + \varepsilon_{1i}^2}; \qquad (2)$$

where

* $t = \frac{Z_{2i} - Z_{1i}}{\sqrt{\varepsilon_{2i}^2 + \varepsilon_{1i}^2}}$

* $\varepsilon_1$ = The individual error of pixel i of first DTM

* $\varepsilon_2$ = The individual error of pixel i of second DTM

* $t$ = is the critical value for a given confidence interval (1.28 for 80% CI, 1.96 for 95% CI); and $z_2$ and $z_1$ are the elevations in a given cell of DEMs.
2.5.2. Multi-scale Model to Model Cloud Comparison-Precision Maps (M3C2-PM)

M3C2-PM is an integrated approach in the CloudCompare software. It is used to directly compare
dense point clouds in 3D instead of comparing DTMs or Mesh. The technique calculates the distance in
the best perpendicular direction to the surface, with a confidence interval normally set at 95% (to assess
the statistical significance of small topographic changes) and a measurement accuracy estimated from
the point cloud roughness, allowing us to detect real topographic changes in elevation and planimetry
between two degraded surfaces by gullying (erosion or sedimentation). The combination of three
factors makes this the preferred method for determining the uncertainty of our digital terrain models:
It can work perfectly even with changes in point density and even in the case of missing data. It
calculates the distance between point clouds along the normal direction of the surface (the calculation
of uncertainty on quasi-vertical areas). It estimates a spatially variable uncertainty for each point cloud
and in accordance with terrain roughness and coregistration error [31,45]. To apply the M3C2 method,
two parameters must be defined:

* 'D': Scale and orientation of the normal are used to calculate the normal area of each point which depends on the roughness and geometry of the point clouds. In the case of bank gullies we first applied a local model of <Quadric >shape along the preferential axis to orient the reference point cloud towards the comparison point cloud. Then we chose the value of D/2 = 10 cm to look for neighbouring points in the reference point cloud while taking into account the terrain roughness.

* 'd': Projection scale, or the diameter of the cylinder in which the average distance of each surface of the point clouds is calculated. In this study we used a value of d = 20cm.

3. Results

3.1. *Data Processing*

3.1.1. Data Precision

In this section we will focus mainly on the application and comparison of different error
quantification techniques related to the collection and processing of SfM-MVS-4D data. Working
on gullies of complex topography, the quantification of errors related to SfM-MVS data along three
types of bank gullies allows on the one hand to analyze the reliability of the methodology adopted for
data collection, and on the other hand to quantify with great precision the sediment budget over a
relatively wet year.

Georeferencing Errors

The mean uncertainty values of the GCPs in XYZ were 1.5 cm (X), 2.4 cm (Y) and 5 cm (Z) (Table 3).
These DGPS data processing results highlighted a flaw in the methodology for collecting GCP point's
network. The points installed at the bottom of the gullies have a low accuracy compared to the others.
They don't exceed 4 cm at best. This can be explained by the poor reception of the GNSS signal at
the bottom of the gullies. The low precision of certain points reduced the general network precision
of the GCPs, but this accuracy of the GCPs does not call into question the precision of the results
since the three study sites were covered by a network of photos of high geometry that testify, by the
total absence of systematic errors, to the shape of the surfaces such as the appearance of domes. In
our case study, to georeference a network of photos with a high overlap, we used a network of GCPs
well distributed spatially but of low precision. This greatly improved the georeferencing of the point
clouds, where the value of the mean vertical RMES XY error is less than 2 cm for the 1st and 2nd sites
and does not exceed 1.3 cm for the third (Table 3). The recessing of the geometry of the aerial photo
network, either using a high rate of overlap, or with the combination of two modes of oblique and
nadir acquisition, allowed us to obtain very precise internal measurements but weakly georeferenced
in the external coordinate system.