Driving Factors of the Hydrological Response of a Tropical Watershed: The Ankavia River Basin in Madagascar

Abstract

Understanding the hydrological behavior of watersheds (WS) and their driving factors is crucial for sustainable water resources management. However, at large scales, this task remains challenging due to the spatial heterogeneity in landscapes, topography and morphology (T), land cover (LC), geology (G), and soil properties (S). In this context, the aim of this study was to identify the key factors that influence the hydrological signatures of four watersheds: Ankavia (WS1: 55% forest cover), Ankaviabe (WS2: 77% forest cover), Sahafihitry (WS3: 41% forest cover), and Antsahovy (WS4: 48% forest cover), over a 10-month study period. These catchments are located within the SAVA region of northeastern Madagascar and have a humid tropical climate. We investigated the relationship between selected catchment descriptors and hydrological signatures by using a Pearson coefficient-based correlation matrix. More specifically, catchment descriptors (extracted from T, LC, G, and S) were correlated with the following hydrological signatures: base flow index (BFI), mean runoff coefficient (rc), mean peak flow (Qp), mean runoff event time scales (ts), high flows (Q5), low flows (Q95), and mean discharge (q_mean). The analysis revealed that land cover, soil properties, and geology seem to be the best predictors for BFI and Q95, while soil properties mainly govern rc, Qp, Q5, ts, and q_mean. These findings provide valuable insights into the key drivers of hydrological behavior that can inform water resource management strategies. In particular, WS2 exhibits better flood buffering capacity but also experiences lower base flows in the dry season, potentially due to higher evapotranspiration. Conversely, WS3 and WS4 (and to a lesser extent WS1) have lower flood buffering capacity, but these watersheds encounter less pronounced low flows in the dry season due to higher BFIs, possibly attributable to lower evapotranspiration rates. The results underscore the importance of responsible land use practices and conservation efforts, which are essential for the sustainable development of the region. By incorporating these practices into water management strategies, we can help ensure a more stable and reliable water supply for communities and ecosystems within the region.

1. Introduction

Global changes, such as population growth, industrial development, agricultural expansion, and climate change, are exerting strong pressures on natural resources in general, and on water resources in particular [1,2,3]. In turn, the resulting stresses on water resources are having a significant impact on production, economic growth, health, livelihoods, and national security [4]. As a result, achieving sustainable management of water resources is a pressing issue, and it will be even more so in the future.

In order to support the sustainable use of water in the face of socio-economic and environmental constraints, understanding the impact of global changes on hydrology is essential [5,6]. This includes an in-depth understanding of the physical drivers of the hydrological response over time and space within the regions to be managed [7]. As water management is often operated across large areas, this represents a significant challenge considering the complex interactions of natural and human processes that may operate concurrently and affect the spatial and temporal variability of the water system [8,9]. This is particularly true in large tropical watersheds because of the rapid pace at which these processes change in response to the growing population pressure, economic development, and climate change. Furthermore, local economies and livelihoods across much of the tropics heavily rely on water availability [10], while the limited economic and development capabilities enhance water-related risks (limitations in water supply, floods, and droughts), making countries from the Global South among the most vulnerable to climate change impacts [11].

Most of our current understandings of hydrological responses are derived from small-scale catchment studies [12,13]. However, in tropical areas, the study of watershed behavior at a small- and meso-scale remains problematic due to a lack of monitoring infrastructure [14,15,16]. As a result, in recent decades, most efforts in tropical areas have been invested in studies of large watersheds rather than small-scale or meso-scale catchments [17,18]. However, large tropical watersheds may be composed of a wide variety of soil types, geological formations, and land covers, making the scaling of hydrological behavior problematic and resulting sometimes in conflicting results [19]. Indeed, whereas some studies conducted in tropical areas have reported considerable hydrological alterations in response to land cover changes in large watersheds [20,21,22], others have shown that the influence of changing forest cover on river hydrology was insignificant [23,24,25]. These inconsistencies may be related to the specificities of the study sites, but also to the large scale of the study of the watersheds and the resulting complexity [22]. More studies at the small- and meso-scale are thus needed to help better understand the hydrological behavior of tropical watersheds. Besides, the challenge to understand catchment response at smaller scales to changing environmental conditions in tropical regions is still very much a topical issue since many climate mitigation strategies are implemented at these smaller scales [7,26].

Madagascar—a tropical island state in southern Africa—is faced with various water-related management issues [27,28,29]. About 80% of the active population’s livelihoods, including agriculture, livestock farming, fishing, and fluvial transportation, are water dependent [27]. The impacts of global changes on hydrology have already been observed, particularly due to the location of the country in the Indian Ocean and its 5000 km coastline, which makes the island extremely vulnerable to climate change and extreme events [30,31]. Floods often occur in northern and central areas, while droughts are recurrent in certain southern regions [28,31]. Over the past decades, to cater to the growing population needs, a significant amount of rainforest cover has been converted into a mosaic of land uses representing various stages of the swidden farming cycle, including severely degraded grasslands [17,32]. Such changes can potentially affect the soil infiltration capacity, leading to increased surface runoff and peak discharge of rivers [33,34]. Unfortunately, the country has very few monitored watersheds, which is a prerequisite to better understanding the hydrological functioning at different scales and designing sustainable water and land use management strategies. Also, creating a credible database for hydrological studies is jeopardized by the highly incomplete and dispersed data from the national weather service stations [31,35].

As a contribution to filling this current knowledge gap, this study aimed to identify the key factors governing hydrological response at the meso-scale using the humid tropical Ankavia catchment in Madagascar as a case study. The data collected will be used in a recently installed rainfall and discharge monitoring network. To reach this overall objective, we address the following specific objectives: (i) to document the Ankavia catchment characteristics; (ii) to analyze the hydrological responses in the catchment; and (iii) to identify the groups of factors (topographic, morphometric, pedo-geological, and land cover properties) that are the most informative to understand these hydrological responses. The study should serve as a basis for the future implementation of an Integrated Water Resources Management (IWRM) system in the SAVA region (northeastern Madagascar).

2. Materials and Methods

2.1. Study Area

The Ankavia watershed (WS1) is situated in the district of Antalaha, a city in the SAVA (Sambava–Andapa–Vohémar–Antalaha) region of northeastern Madagascar (Figure 1). It covers an area of 1116 km² and 55% of the basin is occupied by forest. The watershed lies between 14°90′ and 15°20′ South Latitude and between 49°45′ and 50°15′ East Longitude. The Ankavia river is formed by three tributaries before flowing into the Ankavanana river, which reaches the Indian Ocean. The three upstream sub-catchments are Ankaviabe (WS2; 624 km²; 77% forest cover), Sahafihitry (WS3; 197 km²; 41% forest cover), and Antsahovy (WS4; 70 km²; 48% forest cover). The basin has a total length of 54.6 km, and its elevation varies from 14 m a.m.s.l. at the outlet to 1469 m a.m.s.l. at the highest point on the escarpment (Figure 1). For the purpose of this study, the outlet of the Ankavia was located at the Antsahamanenona Drinking Water Supply Station of Antalaha, which provides drinking water for the city’s 280,000 inhabitants [36]. Additionally, the water-related stresses experienced in the Ankavia basin are similar to those experienced in other basins in the region [37].

According to the meteorological station in Antalaha, the average annual rainfall in the area is 2241 mm (1960–2017). There is a wet rainy season between December and April, while a dry (but still rainy) season occurs from May to November. The daily maximum temperature of 31.4 °C is observed in January, while the minimum temperature of 18.9 °C is observed in August [38]. The western part of the Ankavia watershed is dominated by primary forests, while the eastern part is occupied by a mosaic of natural vegetation composed of shrubs and herbaceous cover [38].

The geology of Madagascar consists of a patchwork of Precambrian rocks, mainly high-grade metamorphic basement domains overlain by unmetamorphosed sedimentary rocks and poorly weathered low-grade metasediments [39]. In the Antalaha area, two different cover sequences lay on high-grade metamorphic and igneous foundation rocks of the Archaean ‘Antongil’ craton and the Neoproterozoic ‘Bemarivo’ belt. The elder of these two cover sequences, the ‘Andrarona’ Group, is composed of low-grade metasedimentary rocks, while the younger sequence, the ‘Ampohafana’ formation, is made up of unmetamorphosed sedimentary rocks [39]. The pedology of the basin is dominated by Ferrasols upstream and Acrisols downstream [40].

2.2. Watershed Attributes

Several satellite products and thematic maps (

Table 1. Summary of the gridded products and maps used to characterize the Ankavia catchment.

Attributes	Data	Spatial Resolution/Map Scale	Sources (Accessed on January 2022)	References
Topography	DEM (Digital Elevation Model)	30 m	https://earthexplorer.usgs.gov/	[41]
Land Cover	GlobCover v2	300 m	https://maps.elie.ucl.ac.be/CCI/viewer/	[42]
Soil Texture	SoilGRIDs	250 m	https://soilgrids.org/	[43]
Dominant Soil Orders	Soil World Reference Base (WRB)	250 m	https://soilgrids.org/	[43]
Hydrological Soil Group	HYSOGs	250 m	https://daac.ornl.gov/SOILS/guides/Global_Hydrologic_Soil_Group.html	[44]
Geology	PRGM	1/100,000	https://www.brgm.fr/fr/reference-projet-acheve/synthese-geologique-miniere-madagascar

) were used to determine the characteristics of the Ankavia catchment and sub-catchments in terms of topography/morphometry, land cover, soil, and geology (Figure 2). The 30 m digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM; [41]) was used to derive the topographic and morphometric attributes (Figure 1). The 2019 land cover of the Ankavia catchment was extracted from GlobCover v2 [42], which has a spatial resolution of 300 m (Figure 2a). The 250 m SoilGRIDs (global gridded soil information) was used to retrieve the soil attributes (Figure 2b), including the organic carbon data at three depths (0–5, 5–15, 15–30 cm), and the dominant soil order based on the World Reference Base (WRB) (Figure 2c) [43]. Hydrologic soil groups were extracted using HYSOGs250m (Figure 2d) [44], which were created to assist with USDA-based (United States Department of Agriculture) curve number runoff modeling at regional and continental scales. The four standard hydrological soil groups (HSG-A, HSG-B, HSG-C, and HSG-D) correspond to soils with increasing levels of runoff potential (low, moderately low, moderately high, and high, respectively). Lastly, the modified geological map (1/100,000 scale) was provided by the Ministry of Energy and Mines through the Mineral Resources Governance Program (PGRM) (Figure 2e).

2.3. Hydrometeorological Datasets

Even though rain gauges were installed as part of the project [45], their spatial coverage was uneven across the watershed, with no gauges in the western part and at the higher altitudes due to accessibility issues (Figure 1). Consequently, the gridded precipitation product IMERG V06 Final was used in the present study [46] (Figure 1). Ref. [45] demonstrated that the IMERG v06 Final data is a reliable alternative to ground-based measurements for event and daily time scales or greater. We used the IMERG product at the event time scale for the time period covered by the discharge measurements, as described below.

To measure the discharge, we installed four gauging stations at the outlet of each catchment (Figure 1). These locations were chosen for their accessibility and the stability of the riverbed sections (predominantly composed of rocks) (Appendix A). Starting on September 2018, all four stations were operational. Water levels were monitored using Solinst^® LTC Levelogger Edge M5 C80 pressure sensors (Level measurement device, Solinst Canada Ltd., Georgetown, ON, Canada) with an accuracy of ±0.3 cm. The atmospheric pressure was recorded using a Solinst^® sensor 3001 LT Barologger Edge, M1.5/F5, with a precision of ±0.05 kPa located at the climatic station of Marofinaritra (Figure 1).

The data were collected with a 1 h timestep (GMT+ 0). The equipment was checked monthly, and data were regularly checked for consistency and completeness. In this study, we used the hydrometric data collected from September 2018 to June 2019, which is the period during which all 4 gauging stations were operational concurrently. After this period, one or more stations were not operational due to damage by extreme floods or technical issues with the probes.

2.4. Catchment Descriptors (CD) Analysis

A wide range (37 in total) of catchment descriptors representing topography/morphometry (T), land cover (LC), soil properties (S), and geology (G) of the Ankavia catchment were acquired and calculated.

Table 2. Description and calculation of catchment descriptors.

Groups	Variables	Label (Units)	Calculation/Description	Reference
Topography-morphometry (T)	Drainage area	A (km2)	GIS calculation 1: area calculation	[47]
	Perimeter	P (km)	GIS calculation: perimeter calculation	[47]
	Basin length	Lb (km)	GIS calculation: spatial analyst tool/hydrology
	Mean slope	Sl (°)	GIS calculation: spatial analyst tool/surface
	Elevation range	Hmin, Hmax (m)	Hmin is the minimum elevation; Hmax is the maximum elevation.
	Hypsometric curve		The shape of a hypsometric curve is an indicator of dominant geomorphic processes at work in a watershed area. A convex curve indicates that a greater part of the watershed’s area is held relatively high in the watershed. A concave curve indicates that the bulk of the watershed’s area resides at a relatively low elevation.	[48]
	Gravelius coefficient	Kg (-)	Kg = P/(2(πA)1/2) 1 < Kg < 1.5: circular shape 1.5 < Kg < 1.9: elongated shape
	Effective basin width	Rb (km)	Rb = A/Lb where A is the watershed area and Lb is the watershed length
	Elongation ratio	Re	Re= 2(A/π)0.5/Lb Re is the ratio of the diameter of a circle having the same area as the watershed to the maximum watershed length. Re < 0.7: more elongated 0.8 < Re <0.9: oval Re > 0.9: circular	[49]
	Unit shape factor	Ru	Ru = Lb/A½ Ru is the ratio of the length of the watershed to the square root of the area.
	Total relief	Rr (km)	Rr = Hmax − Hmin Rr is the maximum vertical distance between the lowest (outlet) and the highest (divide) points on the valley floor of a watershed.
	Relief ratio	Rh	Rh = Rr/Lb where Rr is the total relief and Lb is the watershed length Rh indicates the overall steepness of the drainage watershed and the intensity of erosional processes operating on the slope of the watershed.	[49]
	Total stream length	Lt (km)	GIS calculation: spatial analyst tool/hydrology	[50]
	Length of main stream	Lc (km)	GIS calculation: spatial analyst tool/hydrology	[50]
	Total number of stream segments	Nu	GIS calculation: spatial analyst tool/hydrology
	Drainage density	Dd (km/km2)	Dd = Lt/A Dd is the stream length per unit area in the watershed. Dd < 2: very coarse density 2 < Dd < 4: coarse density 4 < Dd < 6: moderate density 6 < Dd < 8: fine density Dd > 8: very fine density	[47]
	Constant of channel maintenance	Cm (km2/km)	Cm = 1/Dd Cm expresses the watershed surface required to encompass a 1 km long stretch of stream.	[49]
	Ruggedness number	Rn	Rn = Rr × Dd Rn indicates the structural complexity of the terrain in association with relief and drainage density.	[50]
Land cover variables (LC)	Land cover		% of catchment area covered by the following land cover types: mosaic crop (cropland, herbaceous cover) and mosaic natural vegetation (Crop)/mosaic tree and shrub (Mos)/tree cover (Tree).
	Estimated Curve Number (CN)	CN_est	CN_est = Σ (%land_cover(i) × CN(i))/100) where CNi is the respective CN to land cover, which is based on USDA Curve Number. CN_Crop = 77, CN_Mos = 60, CN_Tree = 55.	[51]
Soil variables (S)	Soil texture		% of catchment area covered by soils having the following soil textures: sandy clay (S-C), clay loam (C-L), sandy clay loam (S-C-L).
	Soil types		% of catchment area covered by the following soil types: Acrisols (Acri), Ferrasols (Ferra), Lixisols (Lixi), Cambisols (Cambi).
	Soil hydrological group		% of catchment area covered by the following hydrological soil group: HSG-C (H-C), HSG-D (H-D).
	Organic Carbon density	OC (g/kg)	OC = OC0–5cm × 1/6 + OC5–15cm × 2/6 + OC15–30cm × 1/2
	Grain diameter	Dg (µm)	Dg = %Sand × øsa + %Clay × øsi + %Silt × øc where øsa, øsi, øc are the median particle diameter of sand, silt, clay, which are 900 µm, 25 µm, 1 µm, respectively.
Geological variables (G)	Geological characteristics		% of catchment area covered by the following geological classes: Gneiss (gneiss), Granites (gran.), Shales (shales), Quartzites (quartz.).

Notes: ¹ Geographic Information Tool (GIS) calculations are performed using the ArcToolbox module within the ArcGIS software 10.4 version ^®.

lists the metadata of catchment descriptors.

2.5. Streamflow Analysis

At first, since the atmospheric pressure varies with elevation [52], the atmospheric pressure recorded at the Marofinaritra weather station was corrected using Equation (1) for an isothermal atmosphere to provide atmospheric pressure at the outlet of each catchment.

P(hi) = P(h₀) × e^−Δh/hs,
with hs = R × T/M × g

(1)

where: P(h_i) is the effective atmospheric pressure at elevation h_i [Pa],

P(h₀) is the atmospheric pressure measured at the Marofinaritra station at elevation h₀ [Pa],

Δh is the difference in elevation between h_i and h₀,

R is the molar gas constant [8.314 J K⁻¹ mol⁻¹],

T is the absolute temperature of the air, taken at the Marofinaritra station [°K],

M is the average molar mass of gas in the atmosphere [0.02896 kg mol⁻¹],

g is the gravitational acceleration [9.81 m·s⁻²].

Then, the water level from the Levelogger data was obtained using Equation (2).

Z = (Pi − P(h_i))/ρ × g

(2)

where Z represents the water level [m],

Pi is the measured water pressure [Pa] at a given outlet,

P(h_i) is the atmospheric pressure [Pa] at each sub-catchment outlet,

and ρ is the density of water (1000 kg/m³).

The discharge was determined by converting water level data time series to discharge through the rating curve equation for each outlet (Appendix B). During flow gauging, water velocity for low and medium flows (Q) was measured using an OTT MF PRO magnetic induction current meter (OTT hydromet GmbH, Kempten, Germany), which has a measuring range of 0 to 6 m·s⁻¹ and an accuracy of ±2% of the measured value ±0.015 m·s⁻¹ (between 0 and 3 m·s⁻¹) and ±4% of the measured value ±0.015 m·s⁻¹ (between 3 and 5 m·s⁻¹). At least 20 gauging operations were carried out at each outlet to establish the rating curve, using the measured pairs of water level and discharge (Q_i, Z_i) for low and medium flow conditions. Then, a power function (Equation (3)) was fitted to the water level discharge data:

Q = a × (Z_i)^b

(3)

where a and b are fitting coefficients,

Z_i = Z + Z₀, Z is the water level [m] measured by the Levelogger,

and Z₀ is the stage at zero discharge [m].

Due to safety concerns during high flow conditions and practical limitations, the rating curves for high flows were extrapolated by using an empirical method based on Manning’s equation (Equation (4)).

Q = (1/n) × C_s × H_r^2/3 × S_l^1/2

(4)

where n is the roughness coefficient of the channel [s·m^−1/3], derived from stream gauging for lower and medium flow conditions by combining Equations (3) and (4), C_s is the cross-sectional area of flow [m²], H_r is the hydraulics radius [m], the ratio of the cross-sectional area of fluid flow to the wetted perimeter, and S_l is the slope of the river (%) at the gauging station, estimated from 30 m resolution DEM.

2.6. Hydrological Signatures (HS) Analysis

Hydrological signatures (HS) are quantitative metrics or indices that define statistical or dynamical properties of hydrologic data series [53]. These signatures can be computed using various data sources but, in practice, they are most often calculated using rainfall-runoff time series. In this study, the HS indicators were derived from rainfall-runoff datasets for the 10-month study period.

2.6.1. Event Identification

Rainfall events were defined as periods of non-zero basin-average rainfall interspersed with periods of zero rainfall, considering a minimum cumulative rainfall amount of 2.5 mm and a minimum inter-event time of 6 h [45].

A runoff event is characterized by a runoff start, peak(s), and end. The start of a direct runoff event was defined as the moment in time closest to the (first) peak when total streamflow equals baseflow (Figure 3). The end is reached as soon as total streamflow falls back to the baseflow level (see further). Rainfall and runoff events were then matched, resulting in three cases: (i) rainfall event without runoff event, (ii) rainfall with runoff event, and (iii) runoff without (measured) rain event. In the study, we took into account only the two first cases.

2.6.2. Hydrograph Separation

The Eckhardt method (Equation (5); [54]) was used to separate the discharge data into quickflow (Q_f; the direct flow in response to a rainfall event) and baseflow (Q_bf: the delayed flow from storage) (Appendix C). This method relies on two parameters: the filter parameter (k) and the maximum baseflow index (BFI_max). The filter parameter, which describes the rate at which the streamflow decreases with time following a recharge event, is estimated through recession analysis, while the BFI_max is the maximum value of the baseflow index that can be modeled through the recursive digital filter algorithm [54].

$$Q_{b,t}=\frac{\left(1-{BFI}_{max}\right)\times k\times Q_{b,t-1}+\left(1-k\right)\times {BFI}_{max}\times Q_{s,t}}{1-k\times {BFI}_{max}}$$

(5)

where Q_b,t and Q_b,t−1 are the baseflow at time steps t and t − 1, respectively; Q_s,t is the total streamflow at time step t; BFI_max refers to maximum baseflow index (ratio of baseflow to the total flow), and k is the filter parameter.

In the present study, we employed the Web-based Hydrograph Analysis Tool (WHAT) [55] to separate baseflow from total streamflow. In Equation (5), baseflow for the first time step, Q_b,t−1, was assumed to be 50% of the total streamflow. Ref. [56] proposed different BFI_max values for various aquifer types. Based on the hydrological and geological characteristics of the studied watersheds, we considered the catchments to have porous aquifers and perennial streams. Therefore, the default BFI_max was set to 0.80 while the filter parameter ‘k’ was optimized through recession analysis of selected events [57].

2.7. Hydrological Signatures Analysis

A large number of hydrological signatures (HS) have been proposed in the literature but there are some general selection criteria. According to [58], HS should meet the following requirements: (i) they can be related to hydrological processes to better understand catchment behavior; (ii) they are sensitive to events occurring over different time periods; and (iii) they are not redundant. Based on these criteria, this study used the following HSs: baseflow index (BFI); runoff coefficient (rc); peak discharge (Qp); runoff event time scale (ts); high flows (Q5); low flows (Q95); and mean discharge (q_mean) (Table 3).

BFI is the ratio of total baseflow to total discharge during the observation period. Runoff coefficient rc is the ratio of the quickflow volume of a specific runoff event to the corresponding rainfall. The peak discharge, Qp, represents the maximum rate of discharge during the period of runoff resulting from a rainfall event. The time scale ts is the ratio between quickflow volume and the peak discharge of a specific runoff event. It is an indicator of catchment flashiness and a measure of the importance of fast runoff generation processes, such as overland flow or fast subsurface storm flow. Fast-reacting catchments will have small ts, while catchments with slower runoff generation processes tend to have larger ts. In this study, ts is normalized by the length of the main river. Rc, Qp, and ts were determined for each individual event. Q5 and Q95 were derived from the flow duration curves (FDC), which are cumulative frequency curves that show the percentage of the time that a specified discharge is equaled or exceeded during a given period [59]. The FDC combines the flow characteristics of a stream throughout the range of observed discharges, without regard to the sequence of occurrence. Q5 (the 5th percentile of flow) and Q95 (the 95th percentile) correspond to the discharge that equaled or exceeded 5% (high flows) or 95% (low flows) of the time, respectively. Finally, q_mean is the mean daily discharge during the 10-month study period.

Table 3. Hydrological signatures of runoff events used in the present study.

HS (Unit)	Descriptions	Calculation	References
BFI (-)	Ratio of total baseflow (Qbf) to total discharge (Q) during the observation period; hydrograph separation performed according to [55].	BFI = $\frac{\sum Q_{bf}}{\sum Q}$	[55]
rc (mm/mm)	The ratio of the quickflow volume (Qf) of a specific runoff event [mm] to the corresponding rainfall (Ri). Mean of all events [mm].	rc = $\frac{Q_f}{R_i}$	[60]
Qp (mm/h)	Highest value of the discharge during an event (mean value).		[61]
ts (h/km)	Ratio between quickflow volume (Qf) [mm] and the peak discharge (Qp) [mm/h] multiplied by the length of the main river (Lc) [km] (mean of all events).	ts = $\frac{Q_f}{Q_p\times L_c}$
Q5 (mm/day)	5th percentile of the flow duration curve (high flows).
Q95 (mm/day)	95th percentile of the flow duration curve (low flows).
q_mean (mm/day)	Mean daily discharge

Table 3. Hydrological signatures of runoff events used in the present study.

HS (Unit)	Descriptions	Calculation	References
BFI (-)	Ratio of total baseflow (Qbf) to total discharge (Q) during the observation period; hydrograph separation performed according to [55].	BFI = $\frac{\sum Q_{bf}}{\sum Q}$	[55]
rc (mm/mm)	The ratio of the quickflow volume (Qf) of a specific runoff event [mm] to the corresponding rainfall (Ri). Mean of all events [mm].	rc = $\frac{Q_f}{R_i}$	[60]
Qp (mm/h)	Highest value of the discharge during an event (mean value).		[61]
ts (h/km)	Ratio between quickflow volume (Qf) [mm] and the peak discharge (Qp) [mm/h] multiplied by the length of the main river (Lc) [km] (mean of all events).	ts = $\frac{Q_f}{Q_p\times L_c}$
Q5 (mm/day)	5th percentile of the flow duration curve (high flows).
Q95 (mm/day)	95th percentile of the flow duration curve (low flows).
q_mean (mm/day)	Mean daily discharge

2.8. Linking Hydrological Signature and Catchment Descriptors

A correlation matrix was used to select the most significant CDs. When two or more CDs within the same group of variables (T, LC, S, or G) were strongly correlated (r > |0.5|), only one was retained. The links between event characteristics and selected landscape features were then analyzed to provide insights into the physical drivers of the event runoff response.

Pairwise correlations between individual HSs and selected CDs were determined using the Pearson correlation coefficient and summarized in a matrix. It was displayed as a heat map to ease the examination of these interactions. Additionally, to assess the robustness of the correlation between CD and HS, we checked for the presence of any influential data points (high leverage), which can unduly influence any part of the regression analysis, such as the predicted responses or the estimated slope coefficients. To identify these points, we used the “Difference in Fits” measure, denoted DFFITS (Equation (6)).

$${DFFITS}_i=\frac{\widehat{yi}-\widehat{y\left(i\right)}}{\sqrt{MSE\left(i\right)\times l_{ii}}}$$

(6)

The numerator of Equation (6) measures the difference in the predicted responses obtained when the ith data point is included and excluded from the analysis. The denominator is the estimated standard deviation of the difference in the predicted responses, where MSE is the mean square error, and the leverage l_ii is a measure of the distance between the y value for the ith data point and the mean of the y values for all n data points. Therefore, the difference in fits quantifies the number of standard deviations that the fitted value changes when the ith data is omitted. An observation is deemed influential if the absolute value of its DFFITS value is greater than: $\sqrt[2]{\frac{k+2}{n-k-2}}$, where n is the number of observations and k is the number of predictor terms.

If there are no particular influential data points, a second approach is to delete one observation at a time and refit the regression model on the remaining n−1 observations. We then compare the results obtained using all n-observations with the results obtained by deleting the ith observation to see how much influence that observation has on the analysis. If the predicted responses and estimated slope coefficients are roughly the same, it suggests the absence of influential data points (high leverage), indicating a robust correlation. Otherwise, if the results are significantly different, it implies that the correlation is not robust.

3. Results

3.1. Catchment Descriptor Characteristics

3.1.1. Hypsometric Curves

The hypsometric slope–area curves for the catchments are quite similar, with 80% of the catchment areas having slopes between 5 and 25 degrees, and 20% of the catchment areas with slopes above 25 degrees (Figure 4). The elevation–area hypsometric curves for each catchment are more diverse. WS1 and WS4 have two inflection points, one at 20–30% of the cumulative area, and another at 70–80% of the cumulative area. WS3′s curve is initially fairly linear but becomes concave at higher elevations. WS2′s curve shows a marked S-shape, which seems to indicate that WS2 is in a younger geomorphic stage compared to the other three catchments, which are approaching the mature/old stage [48].

3.1.2. Correlation between CDs

In the first step, the least correlated variables within each attribute were selected based on Figure 5a. As a result, only Slope, Kg, and Dd were retained out of the 14 topographic (T) attributes. Among the four-land cover (LC) attributes, only the CN_est was selected. Ferra, S-C, S-C-L, and OC were chosen from the 11 soil attributes. Finally, Gran and Gneiss were kept among the geology (G) attributes.

Figure 5b was used to assess the between-group correlations, i.e., correlations between CDs from different groups. Kg (T) was found to be positively correlated with CN_est (LC) and negatively correlated with Ferra (S), S-C-L (S), and Gneiss (G). Additionally, there was a strong correlation between Slope (T) and OC (S), as well as between Dd (T) and S-C (S) and between Dd and Gran (G). CN_est (LC) is correlated negatively with both Ferra (S) and Gneiss (G). Finally, S-C (S) shows a positive correlation with Gran (G), and S-C-L (S) is positively correlated with Gneiss (G).

3.2. Hydrological Signatures and Their Spatial Patterns

3.2.1. Rainfall Characteristics

Most rainfall fell between mid-December and mid-April (wet rainy season, Figure 6). The cumulative rainfall amount varied across the catchments during the study period. Specifically, WS4 in the western part received the highest amount of rainfall (~1800 mm), whereas WS3 in the eastern part received the lowest amount of rainfall (~1670 mm).

3.2.2. Streamflow Characteristics

In terms of flow distribution (Figure 7a), WS1 had the highest cumulative flow during the study period, accounting for 66% of total rainfall, while WS2 had the lowest (46%). When considering cumulative quickflow (Figure 7b), WS1 had the highest value (20% of total rainfall), whereas WS2, WS3, and WS4 had fairly similar values (12–15%). In terms of cumulative baseflow (Figure 7c), WS1, WS4, and WS3 had similar values (42–46% of total rainfall), while WS2 had noticeably lower values (31%).

3.2.3. Streamflow Characteristics

The flow duration curve (Figure 8) reveals that the ‘high flows’ (0–5% exceedance probability) for all the watersheds are relatively similar. The curves differ more strongly in the medium flow range (20–80 percent), with WS1, WS3, and WS4 displaying the highest values and WS2 having the lowest. In terms of low flows (95–100% exceedance probability), the WS2 and WS1, in that order, exhibit the lowest, while the WS3 and WS4 have the highest.

3.2.4. Hydrological Signatures

The BFI values range from 0.60 to 0.68, with the watersheds ranked in ascending order, as follows: WS2, WS1, WS3, and WS4 (Figure 9a). Values of rc, Qp, and ts are highly variable across events in all four watersheds (Figure 9b–d). Considering the mean values, we can observe that the runoff coefficient range between 0.09 and 0.34, with WS2 displaying the lowest value, followed by WS3, WS4, and WS1. Mean peak discharge Qp ranges from 0.34 to 1.38 mm/h, with WS2 having the lowest value, followed by WS1 and WS3, while WS4 shows the highest value. Regarding the runoff event time scales, they range between 0.11 and 0.22 h/km, with WS4 having the highest value, followed by WS1, WS3, and WS2 in descending order. Additionally, the low flows Q95 follow a similar trend to the BFI. Finally, the high flows Q5 and mean discharge generated by WS1 are the highest, followed by WS4 and WS3, while the lowest value pertains to WS2 (Figure 9e).

3.2.5. Correlation between Hydrological Signatures

According to Figure 10, there are strong correlations between the hydrological variables. Specifically, BFI exhibits a strong positive correlation with Q95 (low flows). rc is positively correlated with Q5 and q_mean, while Qp is related to ts. Furthermore, there is a strong correlation between Q5 and q_mean.

3.3. Correlation between Selected CDs and HSs

The correlation plot shows links between the selected catchment descriptors (CDs) and the hydrological signatures (HSs). However, not all the strong correlations (indicated by red and green circles) are reliable due to the strong dependence on a single ‘high leverage’ point. Accordingly, we added an ‘*’ to indicate a robust correlation. As a result, each HS was found to be correlated with at least two or more CDs (Figure 11).

Specifically, BFI and Q95 both exhibit a positive correlation with Kg (T) and CN_est (LC), and a negative correlation with Ferra (S), S-C-L (S), and Gneiss (G). rc is negatively correlated with Slope (T) and Ferra (S), while Qp is positively correlated with Kg (T) and negatively with Ferra (S) and Gneiss (G). ts is negatively correlated with S-C (S) and Gran (G). Q5 is negatively related to Slope (T) and OC (S), while q_mean is negatively correlated with Ferra (S).

4. Discussion

4.1. General Hydrological Characteristics

The study period (September 2018–June 2019) can be qualified as a normal year. The cumulative rainfall curves clearly show a distinct wet rainy and dry (but still rainy) season. The hot and wet rainy season occurs from December to April, while the cool and dry season occurs from May to November. During the rainy season, the region receives a large amount of precipitation, but this does not immediately translate into increased river discharge. Indeed, although rainfall increases sharply as of December, notable changes in the cumulative streamflow happen only as of February, which is likely due to a slow response time of the groundwater table. Conversely, during the dry season, the region receives little precipitation, resulting in a period of low flows in the rivers. This strong seasonality of rainfall and river discharge in the northeastern region of Madagascar is the result of a complex interplay between several atmospheric and oceanic processes, such as convective precipitation and trade winds, as well as the movement of the intertropical convergence zone (ITCZ) [62]. Specifically, the position of the ITCZ shifts northward from March to August and southward between September and December, resulting in distinct rainy and dry seasons in tropical regions [63,64].

The proportion of precipitation that turns into streamflow within the Ankavia basin is somewhat variable, ranging from 40% to 60% depending on the watershed. Generally, in humid regions with well-drained soils and gentle topography, such as parts of the eastern United States, Europe, and Southeast Asia, the proportion of rainfall that becomes streamflow is lower, ranging from about 10% to 40% [65,66]. On the other hand, in arid and semi-arid regions characterized by limited vegetation cover and poorly drained soil, the streamflow proportion ranges from approximately 20% to 60% [67]. In tropical regions with highly drained soils and rugged topography, such as those in Asia and Brazil, the proportion of precipitation that flows as streamflow ranges from 30% to 60%, and sometimes can exceed 90% during extreme events [68,69]. Finally, in mountainous areas with high precipitation and steep topography, such as the Andes, the Alps, and the Himalayas, this rate can reach up to 80% or more [70]. Therefore, the discharge-to-rainfall ratio observed in the Ankavia basin appears consistent with observations from tropical regions with highly drained soils and rugged topography.

The calculated BFI values in this study range from 0.60 to 0.68, which is consistent with studies conducted in Brazil [71] and Indonesia [72] in similar environments. However, other regions in Indonesia have shown smaller BFI values (ranging from 0.49 to 0.56) compared to Ankavia [73]. In contrast, higher BFI values were reported for watersheds in Thailand [74] and Malaysia [75].

The study findings show that the mean runoff coefficient values range from 0.09 to 0.34, which is consistent with similar research conducted on meso-scale watersheds in Brazil [71] and in West Africa [76]. However, the observed rc values are slightly higher than the value reported in other parts of Brazil [77], ranging from approximately 0.008 to 0.12, and lower than the value estimated in Malaysia [78], ranging from 0.32 to 0.92. Additionally, some runoff coefficient values greater than 1 were observed in this study, especially during some extreme events, which is theoretically impossible since the amount of rainfall should always be greater than the volume of runoff. This could be due to several factors. Firstly, the rainfall data used in the study were obtained from IMERG and are known to underestimate some classes of rainfall, in particular, extreme rainfall events (>100 mm/d) [45]. Additionally, the method used to separate the hydrographs could also contribute to uncertainty in the estimation of runoff. Moreover, since the rating curve was extrapolated for high water, the runoff during extreme events may be biased (possibly overestimated).

The mean peak discharge and runoff event time scales in the Ankavia watershed range from 0.34 to 1.38 mm/h and 0.11 to 0.22 h/km, respectively. It is important to note that a wide range of values for these two hydrological signatures are found in the literature, which means that there are more dependent on basin characteristics, rather than the specific tropical regions themselves [79,80,81].

When compared to other tropical regions, the Ankavia River in Madagascar has a slightly higher range of flow rate values. For instance, the Sapucaí River in Brazil [68], the Mae Klong River in Thailand [82], and the Mekong River in Vietnam [83] have flow rates (Q5, q_mean, Q95) of (5.48, 2.45, 1.33), (5.6, 2.19, 1.11), and (2.50, 1.74, 1.06) [mm/d], respectively. The region where the Ankavia River is located receives high levels of precipitation, which contributes to its relatively high flow rate even in normal years such as 2019 (~1800 mm). Overall, the Ankavia River exhibits similar discharge values when compared to other basins studied in the SAVA region. For instance, the Lokoho River, which is located to the north of the Ankavia basin, has (Q5, q_mean, Q95) of (7.83, 4.5, 1.94), as reported by ref. [37]. Similarly, the Sambirano River, situated further north, has a Q95 flow of 1.29 mm/d and an average flow of 4.06 mm/d [64]. Hence, the observed range of values in this study falls within the ‘normal’ range for rivers in the SAVA region, Madagascar [64].

4.2. Catchment Descriptors

We examined all the catchment attributes using the correlation matrix based on Pearson’s coefficient (Figure 5a). A critical analysis of these correlations is crucial before deriving any linkages; otherwise, potentially misleading relationships may be inferred [84]. In addition, strong correlations between two or more attributes belonging to the same group may arise because these parameters roughly evaluate the same characteristics. This is for instance the case with the attributes related to catchment area and perimeter (Kg, Cc, Re, etc.) which all determine the watershed morphometry.

The strong relationships that have been found between CDs of different groups in the present study can be explained by the spatial distribution of these CDs in the Ankavia watershed (Figure 2 and Table S1). For instance, WS2 has a more compact and rougher landscape, lower CN (i.e., higher forest cover), more Ferrasols, more SCL soils, and more Gneiss than the other WSs. In contrast, WS3 and WS4 are more elongated, less rough, have higher CN (i.e., lower forest cover), and more Acrisols and CL soils than WS2. The differences in land cover characteristics are likely due to the absence of secondary roads and navigable rivers in WS2, making it less accessible and, therefore, less prone to human activity. Moreover, the presence of natural parks and protected areas in WS2 also probably contributes to its higher forest cover.

The distribution of soils in the watershed depends on various factors such as topography, climate, and geology. In our case study, Acrisols dominate in the downstream (northeast) part, while Ferrasols dominate in the upstream (southwest) part of the watershed. Typically, Acrisols are found on erosive slopes of low hills while Ferrasols are found on adjacent stable pediments or high hills [43].

4.3. Which Catchment Descriptors Are the Best Predictors for Each Hydrological Signature?

The BFI calculated in this study ranges from 0.60 (WS2) to 0.68 (WS4) and is positively correlated with CN_est, the latter being negatively correlated with the forest cover. Land cover attributes thus appear to be an important predictor of baseflow (BFI). Indeed, forest cover has a higher evapotranspiration rate than cropland, degraded land, and tree fallow cover, which may explain why BFI declines as forests cover increases [85,86]. This is consistent with other studies conducted in Brazil [71] and Indonesia [72]. However, some studies contradict this finding, such as a study conducted in the southern Amazon [71], which reported a higher BFI in forested areas than in grazing areas, supporting the “sponge effect theory” [19] that forested watersheds produce more baseflow during the dry season than watersheds with disturbed vegetation [87].

BFI also appears to be correlated with Kg and the relative area occupied by Ferralsols (and sandy clay soils). The positive correlation with Kg is likely the result of a confounding effect since Kg is highly correlated with CN_est (the most forested watershed is also the most compact). The negative correlation with Ferralsols may also be spurious, as Ferralsols are negatively correlated with CN_est. Indeed, Ferralsols are known to have favorable physical properties, especially high infiltration capabilities [88]. A higher BFI would, therefore, be expected in Ferralsol-dominated watersheds, which is not observed here. Similarly, the negative correlation with gneiss may also be spurious, especially since fractured hard rocks such as quartzites and gneisses more generally lead to aquifers reservoirs with high capacities [56].

The mean runoff coefficient values in our study range from 0.09 to 0.34, depending on the watershed. The WS2 watershed, which has a higher proportion of Ferralsols and forest compared to the other basins (Figure 2c), has the lowest runoff coefficient. Ferralsols are known for their excellent physical properties, including high permeability due to their well-developed microaggregate structure, sometimes referred to as “pseudo-sands” [89]. Conversely, Acrisols, present mainly downstream of the basin, have a poorer structure and horizons with clay accumulation, which may explain higher runoff, particularly if the soil is cultivated. While the correlation is not strong, we observe that the runoff coefficient is lower in WS2 (Figure 9b), which has the most forested cover, and higher in the other less forested WSs. This observation is consistent with the well-known role of forest cover in intercepting rainwater (litter and canopy) and protecting and improving the hydrological conditions of the soil thanks to the litter layer that accumulates on the surface. The average slope is negatively correlated with the runoff coefficient, but this relationship is most likely spurious since the range of average slopes across the watersheds is small (~2°) and topographic parameters may not significantly affect the runoff coefficient [8]. Other studies have also demonstrated a strong correlation between runoff coefficient, soil properties, and land cover characteristics, such as those conducted by ref. [71] in Brazil and by ref. [76] in West Africa.

The mean peak discharge in the Ankavia catchment ranges from 0.34 mm/h to 1.38 mm/h. The peak discharge of WS2 is found to be the lowest compared to the other three WSs. In this case study, there appears to be a positive correlation between Qp and Kg (T); however, this relationship may be spurious. Indeed, more compact watersheds are expected to lead to shorter response times and higher peak discharges. Instead, the relationship between Qp and runoff/quickflow (rc) seems more relevant. Indeed, the positive relationship between these two parameters is well-established in the literature, indicating that if runoff increases (rc), Qp also increases [47]. Therefore, predictors of runoff coefficients, such as soil type, are also predictors for peak discharge.

The WS1, WS2, and WS3 catchments exhibit faster hydrologic responses (lower runoff event time scale ts) compared to WS4, which is located towards the west. Our correlation analysis suggests a strong positive relationship between Qp and ts (Figure 10). Meanwhile, soil attributes, which are critical factors for runoff, also appear to contribute to the hydrologic responsiveness of the catchments (Figure 11). Additionally, we observed that WS4 has the lowest proportion of granite and the highest proportion of shale.

The correlation analysis reveals that Q5 and q_mean are positively correlated with rc and that Q95 is positively correlated with BFI (Figure 10). Q5 (high flows) is the result of runoff, and factors that are good predictors of rc, such as soil attributes, are therefore also good at predicting Q5. It is worthing that the mean discharge values are heavily influenced by extreme values. Hence, soil attributes that are predictors of rc and Q5 also serve as good predictors for q_mean. Q95 (low flows) originate from soil releases, groundwater, or channel storage; hence, catchment descriptors for predicting the BFI and Q95 are expected to be roughly the same (Figure 11).

In summary, the more forested WS2 watershed has a better flood buffering capacity (lower Q5 and Qp values), but lower BFI (potentially due to higher ET). Consequently, this watershed experiences more pronounced low flows during the dry season (lower Q95 values). In contrast, the less forested WS3 and WS4 watersheds (and to a lesser extent WS1) have higher Q5 and Qp values, which indicates a lower flood buffering capacity than WS2. Nevertheless, they have higher BFI than WS2, resulting in less pronounced low flows (higher Q95) during the dry season.

5. Conclusions

In this study, we aimed at gaining a better understanding of the factors driving the hydrological response of the Ankavia watershed by investigating the linkages between watershed descriptors and hydrological signatures using a fairly unique dataset for northeastern Madagascar. Our findings indicate that all watershed attributes play a role in governing hydrological responses to a greater or lesser extent. Land cover, soil properties, and geology are the best predictors of baseflow and low flows (Q95). In contrast, soil properties mainly govern the runoff coefficient (rc), mean discharge (q_mean), peak flow (Qp), and runoff event time scales (ts).

The findings indicate that the WS2 watershed, primarily covered by forests, has the lowest runoff coefficients and peak flows. This emphasizes the importance of protecting the forest cover in this area, not just to protect biodiversity but also to mitigate downstream flooding risks. However, this watershed is also characterized by the lowest low flows (Q95), reflecting a limited capacity for maintaining sufficient discharge during the dry season. Furthermore, the notably higher peak discharges in WS3 and WS4 vs. WS2 highlight the impact of land cover on the hydrological characteristics of the river basin, which suggest that human activities, such as deforestation, shifting cultivation, as well as agriculture, can have a considerable impact on the hydrological behavior, and must be considered in water resource management in the SAVA region.

This study highlights the complexity of the interactions between the factors governing hydrological systems in a large tropical watershed. Nevertheless, the findings can serve as a basis for interpreting hydrological responses in other ungauged basins within the SAVA region but also more generally in the northeastern part of Madagascar that are under similar hydrological conditions. Ultimately, the study provides valuable insights into the hydrological characteristics of the sub-basins within the Ankavia watershed and underscores the significance of sustainable land use practices and conservation efforts for effective watershed management.