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

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,5,6,7,8]. 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,9]. This includes an in-depth understanding of the physical drivers of the hydrological response over time and space within the regions to be managed [10]. 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 [11,12]. 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 [13], while limited economic and development capabilities enhance water related risks (water supply, floods and droughts), making countries from the Global South among the most vulnerable to negative climate change impacts [14,15,16].

Most of our current understanding of hydrological responses are derived from small-scale catchment studies [17,18]. However, in tropical areas, the study of watershed behavior at small- and meso-scales remains problematic due to a lack of monitoring infrastructure [19,20,21]. 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 [22,23]. However, large tropical watersheds may be composed of a wide variety of soil types, geological formations and land uses, making the scaling of hydrological behavior problematic and resulting sometimes in conflicting results [24]. Indeed, whereas some studies conducted in tropical areas have reported considerable hydrological alterations in response to land cover changes in large watersheds [25,26,27], others have shown that the influence of changing forest cover on the river hydrology was insignificant [28,29,30]. These inconsistencies may be related to the specificities of the study sites, but also the large scale of study of the watersheds and the resulting complexity [27]. More studies at the small- and meso-scale are thus needed to help understand better the hydrological behavior of tropical watersheds. In addition, the challenge to understand catchment response at smaller scales to changing environmental conditions in tropical regions is still very much a tropical issue since many climate mitigation strategies are implemented at these smaller scales [31,10].

Madagascar – a tropical island state in southern Africa – is faced with various water-related management issues. About 80% of the active population’s livelihoods, including agriculture, livestock farming, fishing and fluvial transportation, are water-dependent [32]. 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 [33]. Floods often occur in northern and central areas, while droughts are recurrent in certain southern regions [16,33]. Over the past decades, to cater for 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 [22,34]. Such changes can potentially affect the soil infiltration capacity leading to increased surface runoff and peak discharge of rivers [35,36]. Unfortunately, the country has very few monitored watersheds, which is a prerequisite to better understand the hydrological functioning at different scales and to design 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 [37,33].

To contribute to fill this current knowledge gap, this study aimed to identify the key factors governing the hydrological response at the meso-scale using the humid tropical Ankavia catchment in Madagascar as a case study. Use will be made of data collected 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 region of northeastern Madagascar (Figure 1). It covers an area of 1,116 km², and 55% of the basin is occupied by forest. The watershed lies between 14°90' and 15°20' South Latitude, and 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 1,469 m a.m.s.l. at the highest point on the escarpment (Figure 1). In 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 [38]. Moreover, the water-related stresses experienced in the Ankavia basin are similar to those experienced in other basins in the region [39].

According to the meteorological station in Antalaha, the average annual rainfall in the area is 2,241 mm (1960–2017). There is wet rainy season between December and April, while a dry (but still rainy) season occurs from May to November. The maximum temperature of 31.4°C is observed in January, while the minimum temperature of 18.9°C is observed in August [40]. 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 [40].

The geology of Madagascar consist of a patchwork of Precambrian rocks, mainly high-grade metamorphic basement domains overlain by unmetamorphosed sedimentary rocks and poorly weathered low-grade metasediments [41]. 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 [41]. The pedology of the basin is dominated by Ferrasols upstream and Acrisols downstream [42].

2.2. Watershed attributes

Several satellite products and thematic maps (Table 1) 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 Shuttle Radar Topography Mission (SRTM; [43]) was used to derive the topographic and morphometric attributes (Figure 1). The 2019 land cover of the Ankavia catchment was extracted from GlobCover v2 [44], 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 soil class distribution based on the World Reference Base (WRB) (Figure 2c) [45]. Hydrologic soil groups were extracted using HYSOGs250m (Figure 2d) [46], which were created to assist with USDA-based (United states Department of Agriculture) curve-number runoff modeling at regional and continental. The four standard hydrological soil groups (HSG-A, HSG-B, HSG-C, and HSG-D) correspond to soils with varying 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 [47], 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 [48]; (Figure 1). Ref. [47] 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. Starting from September 2018, all four stations were operational and monitored water levels using Solinst® LTC Levelogger Edge M5 C80 pressure sensors 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-hour 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 use (LU), soil properties (S), and geology (G) of the Ankavia catchment were acquired and calculated. Table 2 lists the metadata of catchment descriptors.

2.5. Streamflow analysis

At first, since the atmospheric pressure varies with elevation [54], the atmospheric pressure recorded at the Marofinaritra weather station were corrected using Equation 1 for an isothermal atmosphere to provide atmospheric pressures 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 Marofinaritra 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^-1.mol^-1],

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^-1],

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 A]. During flow gauging, water velocity for low and medium flows (Q) was measured using an OTT MF PRO magnetic induction current meter, which has a measuring range of 0 to 6 m.s^-1 and an accuracy of ± 2% of the measured value ± 0.015 m.s^-1 (between 0 and 3 m.s^-1) and ± 4% of the measured value ± 0.015 m.s^-1 (between 3 and 5 m.s^-1). 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^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 Equation 3 and Equation 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 is the slope of the river (%) at the gauging station, estimated from 30m 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 [55]. These signatures can be computed using various data sources but, in practice, they are most often calculated using rainfall-runoff time series and referred to as event indices or event characteristics. In this study, the HS indicators were derived from rainfall-runoff datasets for the 10-month study period at event time scale.

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 hours [47].

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; [56]) 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 B]. 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 [56].

$$Q{t}_{}=\frac{\left(1-BFImax\right)\ast \mathrm{k}\ast Qb,t-1+\left(1-\mathrm{k}\right)\ast BFImax\ast Qs,t}{1-\mathrm{k}\ast BFImax}$$

(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) [57] 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. [58] 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 [59].

2.7. Hydrological signatures (HS) analysis

A large number of hydrological signatures HS have been proposed in the literature but there are some general selection criteria. According to [60], 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 time scale ts is the ratio between quickflow volume and the peak discharge. 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. Q5 and Q95 were derived from the Flow Duration Curves (FDC) which are cumulative frequency curve that show the percentage of the time that a specified discharge is equaled or exceeded during a given period [61]. The FDC combines the flow characteristics of a stream throughout the range of observed discharges, without regard to the sequence of occurrence. Q5 (the 5^th percentile of flow) and Q95 (the 95^th percentile) correspond to the discharge that is 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.

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, LU, 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. 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.

$$DFFITSi=\frac{\widehat{yi}-\widehat{y\left(i\right)}}{\sqrt{MSE\left(i\right)hii}}$$

(6)

The numerator of Equation 6 measures the difference in the predicted responses obtained when the i^th 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 h_ii is a measure of the distance between the y value for the i^th 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 i^th 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 i^th observation to see how much influence that observation has on the analysis. If the predicted responses and estimated slope coefficients are roughly the same, then there is no influential data point (high leverage), and the correlation is considered robust. Otherwise, if the results are significantly different, 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 cumulative area, and another at 70–80% of 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 mature/old stage [50].

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 use (LU) 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 (LU) and negatively correlated with Ferra (S), S-C-L (S), and Gneiss (G). Also, 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 (LU) 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, 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, 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, WS1, WS4, and WS3 had similar values (42–46% of total rainfall), while WS2 had noticeably lower values (31%).

Figure 7. Streamflow characteristics of sub-catchments (%): cumulative streamflow, cumulative runoff (quickflow), cumulative baseflow from September 2018 to June 2019. All values are expressed relative to the total rainfall of each sub-catchment (Fig.6). Date in X-axis in MM/DD/YYYY format.

Figure 7. Streamflow characteristics of sub-catchments (%): cumulative streamflow, cumulative runoff (quickflow), cumulative baseflow from September 2018 to June 2019. All values are expressed relative to the total rainfall of each sub-catchment (Fig.6). Date in X-axis in MM/DD/YYYY format.

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 9). The mean rc values 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. The runoff event time scales range between 0.11 and 0.22 h/km, with WS4 having the highest value, followed by WS1, WS3, and WS2 in descending order. Also, 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 WS2 has the lowest value.

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 ‘outlier or high leverage’ point. Accordingly, we added the ‘*’ indication to display robust correlations. 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 (LU), 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 year 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 from December, notable changes in the cumulative streamflow happen only as from February, which is likely due to a slow response time of the ground water 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) [64]. 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 [65,66].

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% [67,68]. In contrast, in arid and semi-arid regions with poorly drained soils and steep topography, like those found in Africa and Australia, the streamflow to rainfall ratio ranges from approximately 20% to 60% [69,70]. In tropical regions with poorly drained soils and rugged topography, like 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 [71,72]. 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 [73]. Therefore, the discharge to rainfall ratio observed in the Ankavia basin appears consistent with observations from tropical regions with poorly 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 [74] and Indonesia [75] in similar environments. However, other regions in Indonesia have shown smaller BFI values (ranging from 0.49 to 0.56) compared to Ankavia [76]. In contrast, higher BFI values were reported for watersheds in Thailand [77] and Malaysia [78].

The study findings show that the mean runoff coefficient values range from 0.09 to 0.34, which is consistent with similar research conducted in Brazil [74] and in West Africa [79]. However, the observed rc values are slightly higher than the value reported in other parts of Brazil [80], ranging from approximately 0.008 to 0.12, and lower than the value estimated in Malaysia [81], ranging from 0.32 to 0.92. Also, 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) [47]. 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 wide range of values for these two hydrological signatures are found in the literature, which mean that there are more dependent on basin characteristics rather than the specific tropical regions themselves [82,83,84].

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 [71], the Mae Klong River in Thailand [85], and the Mekong River in Vietnam [86] 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 [66]. Hence the observed range of values in this study falls within the ‘normal’ range for rivers in the SAVA region, Madagascar [66].

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 [87]. 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 contributes probably 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 [45].

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 [88,89]. This is consistent with other studies conducted in Brazil [74] and Indonesia [75]. However, some studies contradict this finding, such as a study conducted in the southern Amazon [74], which reported a higher BFI in forested areas than in grazing areas, supporting the "sponge effect theory" [24] that forested watersheds produce more baseflow during the dry season than watersheds with disturbed vegetation [90].

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 [91]. 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 reservoir with high capacities [58].

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" [92]. 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 [11]. Other studies have also demonstrated the strong correlation between runoff coefficient, soil properties and land cover characteristics, such as those conducted by [74] in Brazil and by [79] 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 [49]. 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). Also, 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 base flow and low flows (Q95). In contrast, soils 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 at maintaining sufficient discharge during the dry season. Furthermore, the notably higher peak discharges in WS3 and WS4 vs WS2 highlights the impact of land use 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 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 in Madagascar that are under similar 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.