Original Research

Hydrological response in a savanna hillslope under different rainfall regimes

Johan van Tol, Stefan Julich, Darren Bouwer, Edward S. Riddell

Received: 11 Oct. 2019; Accepted: 22 Apr. 2020; Published: 29 Oct. 2020

Copyright: © 2020. The Author(s). Licensee: AOSIS.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Soil water is a link between precipitation and the functioning of ecological systems. It is therefore critical to understand exactly how soil water regimes are affected by changes in precipitation. This is especially true for the variable water regimes of savanna ecosystems. Therefore, understanding the effects of precipitation on soil water was the central goal of this article. The hydropedological behaviour of a catena in the Stevenson Hamilton Research Supersite of the Kruger National Park was configured as a conceptual model of catchment modelling framework, a toolbox of various model structures and processes. The model was parameterised using measured hydraulic properties of the soils, and calibrated and validated using measured soil matric potentials and derived actual evapotranspiration (aET) data. The model was then used to simulate hydrological response under five different rainfall scenarios, ranging from 30% drier than the normal rainfall to 30% wetter than the normal rainfall. The scenarios also included rainfall years with fewer but larger rain events, that is, more intense rainfall events. In general, the model performed well with Pearson’s correlation coefficient (R) values ranging between 0.66 and 0.87 and between 0.58 and 0.69 for correlations with daily soil matric potential and daily aET, respectively. Scenario analysis indicates non-linearity in the response of hydrological processes to changes in precipitation. This is especially evident in a seven-fold increase in the duration of saturation at the seepage surface associated with a 30% increase in rainfall. In general, the impact of drying conditions (30% below average rain) has a greater influence on soil water contents, overland flow and percolation from the riparian zone to bedrock than a 30% increase in rainfall would have on the same process.

Conservation implications: This article presents realistic predictions of the potential impact of changes in precipitation on hydrological processes in an important area of the Kruger National Park. These predictions would enable decision-makers to be prepared for the anticipated changes in near-surface hydrological processes associated with climate changes.

Keywords: Hydropedology; Hydrological modelling; Soil water; Catchment modelling framework; Kruger National Park.

Introduction

Savannas are currently experiencing extensive population and agricultural water resource pressure globally (Rockström et al. 2014). Because we are now moving into an era of adaptive management of ecosystems where information on state change thresholds is a prerequisite for their sustained management, it is necessary to invest in understanding the biophysical processes that maintain these systems. In savanna landscapes, soil water is the direct link between precipitation and ecological patterns (Weltzin et al. 2003). The co-domination by trees and grasses is bound by soil water availability, a key factor in establishing form and function. Therefore, spatial soil moisture dynamics is a crucial link in the equilibrium between climate, soil and vegetation in these systems (Rodriguez-Iturbe et al. 1999). The resulting hydrological processes contribute to the biophysical template of these semi-arid systems, controlling the distribution of water and other resources along a continuum within the landscape, which often comprises compound effects of non-linear relationships and threshold-triggered responses. Given this complexity, novel interdisciplinary approaches should be sought to understand hydrological processes in such heterogeneous landscapes (Troch et al. 2008). Moreover, inter-disciplinarity is increasingly valuable for successful landscape management, given the emphasis on hydrological connectivity at landscape scale (Michaelides & Chappell 2009).

The general framework that process hydrologists have used is to garner knowledge of water flow pathways and residence times, as this is essential for predicting the response of a catchment to rainfall (Uhlenbrook, Wenninger & Lorentz 2005). It is suggested that hydropedology is a vital interdisciplinary science for contributions to the earth’s critical zone (Lin 2009). Whilst the fields of pedology and hydrology have long been mutually exclusive, there is a need for the soil scientist to benefit from flow theory when transcribing qualitative descriptions into quantitative expressions, and vice versa, for the hydrologist to develop representative pedotransfer functions in hydrological modelling (Bouma 2006). For instance, Ticehurst et al. (2007) demonstrated that marrying these two disciplines, most notably using soil colour and presence of redox concretions, was effective in indicating depths of saturation and lateral flow occurring on hillslopes in New South Wales, Australia. Testing their conceptual flow path model against hillslope hydrometric observations yielded a satisfactory agreement, but they warned that further catchment information should be sought to reduce model uncertainty.

Meanwhile, Sivapalan (2003) asked why watershed hydrological responses were seemingly simple and hillslope processes comparatively complex, suggesting that by aggregating hillslope complexities into dominant processes, one could parameterise the hillslope as the basic unit within a catchment model. Ridolfi et al. (2003) highlighted these hillslope complexities, including inter alia the spatial heterogeneity of soil properties; climatic variability (although generally uniform at the hillslope scale, may trigger mechanisms within the hillslope to alter its spatial dynamics); lateral redistribution of water along a hillslope because of the formation of a saturated zone within the soil; lateral sub-surface flow in the unsaturated zone; the longitudinal form and profile of a hillslope; and boundary conditions at the valley bottom of the hillslope. It has been advocated that the use of perceptual hydrological models, whilst being largely qualitative conceptualisations, offers good potential based on process understanding of key zones or ‘reservoirs’, and that these ‘soft’ data could be married to ‘hard’ hydrological observations (streamflow and soil water content) to facilitate the reduction in parameter uncertainty (Bouwer et al. 2015; Lorentz et al. 2003; Seibert & McDonnell 2002).

Arid and semiarid regions of the world are highly dependent on the availability of water, affecting recruitment, growth and reproduction, nutrient cycling and net ecosystem productivity. Quantifying the potential consequences of changing precipitation patterns is important in these regions. It is expected that varying regional precipitation regimes will have important ramifications for the distribution, structure, composition and diversity of plant, animal and microbe populations and communities (Easterling et al. 2000; Houghton et al. 2001; Weltzin & McPherson 2003). Although some effects have been estimated, for instance, increasing precipitation variability – but not necessarily the total amount of precipitation – may reduce grassland productivity and livestock-carrying capacity, exacerbate overgrazing and increase rangeland susceptibility to invasive weed species (Knapp et al. 2002).

Given the above rationale, the following three primary objectives are addressed in this article: (1) To configure and parameterise a hillslope-scale hydrological model using soil and environmental data by Bouwer et al. (2020); (2) to calibrate and/or validate the model using hydrometrical measurements presented in Riddell et al. (2020); and (3) to use the calibrated model to simulate hillslope hydropedological processes under different rainfall scenarios.

Site description and methodology

Site characteristics and soil information

A granitic catena (hillslope) of a third-order stream in the Southern Granite Supersite (Smit et al. 2013) in the Kruger National Park (KNP), South Africa, was the focus area in this article. The mean annual precipitation is 560 mm/a (Smit et al. 2013), and the geological formation is granite and gneiss of the Nelspruit suite (Venter 1990).

Vegetation on the crest to upper-midslope positions is dominated by woody vegetation such as Combretum apiculatum and Combretum zeyheri. A prominent seep line exists on the transition between the midst to footslope, which is marked by Terminalia sericea. Below the seep line, there is a long sodic site frequented by Euclea divinorium. Along the footslope, fine-leaved woody species such as Vachellia nilotica dominates (Smit et al. 2013).

Bouwer et al. (2020) conducted a hydropedological survey where 49 soil profiles were classified and interpreted in terms of their hydropedological behaviour (Figure 1). This conceptual model of hillslope hydropedological behaviour served as a basis for both the model set-up and identification of dominant soils and horizons for undisturbed sampling for soil physical measurements. Watermark sensors were installed in representative soil horizons (see Riddell et al. 2020).

FIGURE 1: Hydropedological behaviour inferred from soil distribution pattern, different landscape units and location of watermark sensors on the studied hillslope.

Four modal profiles were identified based on terrain position and soil type. These modal profiles were sampled at 100-mm depth intervals for particle size distribution using the pipette method and chemical analysis. Soil hydraulic properties were determined on five of the dominant horizons (orthic A, yellow-brown apedal B, prismacutanic B, melanic A and pedocutanic B). The saturated hydraulic conductivity of each diagnostic horizon in the modal profiles was measured in duplicate using a modified Bouwer and Rice (1976) double-ring falling-head method:

where K_s is the saturated hydraulic conductivity (mm hr⁻¹), L is the thickness of the horizon, t is the time until constant infiltration and h₀ and h₁ are the head of water above the surface before and after time t, respectively. Water retention characteristics were measured for each horizon on undisturbed cores. The samples were de-aired at –70 kPa at room temperature for 24 h. Water was then allowed to gradually fill the cores. The saturated samples were mounted on a hanging water column set-up in accordance with the procedure of Dirksen (1999). The suction level (h) was set at 0, 38, 50, 100, 200, 400, 600, 800 mm to capture important lower suctions of the retention curves. The gravimetric water content was determined at each level. The samples were then oven-dried at 105°C to determine the bulk density to calculate the volumetric water content (θ_v). The measured bulk density, water retention measurements, together with the particle size distribution, was used to derive Van Genuchten parameters in the retention curve (RETC) program (Van Genuchten, Leij & Yates 1991). These parameters were used as inputs in the hydrological model.

Model configuration and parameterisation

The Catchment Modelling Framework (CMF; Kraft et al. 2011) was used to construct the framework for the hydrological modelling of this hillslope. The CMF is essentially a toolbox to configure a wide range of different model structures based on the finite volume approach. Water fluxes through the landscape are represented as a network of storages and boundary conditions in CMF. Flux-governing equations could be assigned to link the storage units with the next one. These equations can either be fairly simple, for instance, linear storage or tipping bucket approaches, or more complex, for instance, solving of kinematic wave or Richard’s equation. The compounds of the model are assembled using the scripting language Python. Here we used the open-source platform, Spyder version 3.3.3, as the environment for scientific programming in Python (Python Software Foundation 2019).

We divided the 480-m long slope into 24 cells, each covering 2 000 m². These cells were 100-m wide and stretched 20 m downslope. The elevation profile was determined with a differential Global Positioning System (GPS; Riddell et al. 2020). The hydraulic parameters of soils (Table 1), and their spatial distribution presented by Bouwer et al. (2020: Figure 1), guided the spatial distribution of the hydraulic parameters (Figure 2). The vertical discretisation within each cell was 10 cm and the Richard’s equation was used to estimate vertical flow through profiles. Overland flow was calculated using the kinematic wave approach. Water ponding on the surface flows to the next downslope cell, from where it is allowed to re-infiltrate the soil. The thickness of horizons, as observed in the soil survey, was included in the set-up. Below the soils, a 2 m thick ‘bedrock’ layer was included in the model set-up. The hydraulic properties below each of the cells were adjusted as part of the calibration, but in general it reflected the conceptual understanding of ‘recharge’ of groundwater in the crest and riparian zone, and limited deep drainage in the midslope positions. Below the bedrock layer, free drainage into a groundwater store was allowed.

FIGURE 2: Hillslope configuration and hydrological parameters/processes employed in catchment modelling framework for the simulations.

An atmospheric boundary condition was applied to the surface of each cell. Measured climatic parameters (Riddell et al. 2014), including rainfall, minimum and maximum temperatures, relative humidity, solar radiation and wind speed, were included in the simulation on a daily time-scale. A summary of these parameters is presented in Table 1. Potential Evapotranspiration (ET₀) was calculated from the climatic parameters using the Penman–Monteith equation, and the Shuttleworth–Wallace approach was used to distinguish between actual evaporation and transpiration (Shuttleworth & Wallace 1985).

The simulation period commenced on 01 May 2011 and ended on 13 January 2013. Evaluation of the model outputs commenced on 03 November 2011 and 10 November 2011 for evaporation and soil water contents, respectively, because of the availability of reliable data. This allowed the model to ‘settle’ for approximately six months prior to evaluation.

Evaluation and calibration of the model

Soil matric pressures were measured at different horizons of the four sites on the hillslope (Figure 1) and converted to daily average values. Details on the instrumentation and calibration of these measurements are discussed in more detail in Riddell et al. (2014) and Riddell et al. (2020). Actual evapotranspiration (aET) was obtained from Van Eekelen et al. (2015) using the surface energy balance algorithm for land (SEBAL) model (Bastiaanssen et al. 1998). The weekly measurements were then converted into daily values for different landscape positions using the approach reported in Riddell et al. (2020).

Measured daily soil matric pressure and daily SEBAL-derived aET values were used in a manual calibration process. In this process, parameters for which no measured values were available have been varied stepwise in a defined range (Table 2). After each simulation, the Pearson’s correlation coefficient (R) values and root mean square errors (RMSE) between measured and simulated values were calculated. The calibrated parameters comprise the conductivity and porosity of the bedrock layers as well as rooting depth and leaf area index (LAI) of the respective vegetation types. For example, if the simulated water contents of sub-soils accumulated more than what was measured, the conductivity and porosity of the bedrock was increased. The aim of the calibration was not simply to improve model outputs but to maintain the measured parameters and make the unmeasured parameters (i.e. bedrock hydraulic parameters, rooting depth and LAI) more realistic and congruent with existing data for the region and reduce parameter uncertainty.

Precipitation scenario analysis

The calibrated model was then used to simulate hydrological processes under five different precipitation scenarios. An average rainfall year (583 mm), represented by the 2009–2010 rain season (i.e. August to July), was identified from the Riddell et al. (2014) database. This rainfall year represented the ‘average’ scenario, or the baseline. The rainfall was then manually changed to reflect a ‘dry year’ by reducing the measured daily rainfall by 30%. Similarly, a ‘wet year’ was created by increasing the measured rainfall by 30%. ‘Dry extreme’ and ‘wet extreme’ scenarios were created by lumping the rainfall of consecutive rainfall days (i.e. >1 mm day⁻¹) and using this value as the daily rainfall on the first day of the rainfall sequence whilst deleting rainfall on the other days. This was carried out to reflect an increase in rainfall intensity and also longer dry periods. These adjustments were applied to both ‘dry’ and ‘wet’ years to create ‘dry extreme’ and ‘wet extreme’ scenarios. A summary of this climatic data is presented in Table 3.

The soils were initially dry for the model runs, with a saturated zone 10 m below the lowest node of the bedrock; then, different yearly climate inputs were repeated twice. The first year was to allow a settling period for the model and the second year for evaluation purposes. The simulation results presented are centred on the key hydrological processes as presented in Figure 1, namely, storage of water in the crest zone, seepage above the sodic site, overland flow in the midslope, and recharge and storage in the riparian zone.

Ethical considerations

The full ethical clearance statement was obtained from the Interfaculty Animal Ethics Committee of the University of the Free State – Clearance number UFS-AED2019/0121.

Results and Discussions

Hydraulic parameters of the different soils

The crest is dominated by freely drained Clovelly (Cv) soils consisting of an orthic A (ot) horizon on a yellow-brown apedal (ye) B horizon. These soils are marked by high hydraulic conductivity (> 50 mm h⁻¹; Table 4). In the upper-midslope, an albic horizon (gs) was formed, predominantly because of lateral movement on the soil–bedrock interface. The difference in conductivity between the albic horizon and impermeable bedrock (Table 4) promotes the generation of lateral flow and the formation of albic properties (Van Tol, Hensley & Le Roux 2013).

The conductivity of both the top soil and sub-soils of the sodic site (see the orthic horizon above prismacutanic [pr] B horizon – Table 4) restricts vertical infiltration, especially during high-intensive rain events, and promotes the generation of overland flow. The well-developed structure of the melanic A horizon (ml) of the Bonheim soils in the footslope is likely to allow infiltration and storage of water (Figure 1) with the pedocutanic (vp) also allowing drainage towards the bedrock. In the valley bottom, the Dundee (Du) soils have the same hydropedological behaviour as the Clovelly soil in the crest (Bouwer et al. 2020). The ot and the stratified alluvium (al) horizons of the Dundee form were consequently assigned similar hydraulic properties as that of ot and ye of the Clovelly soil, respectively (Table 4). Infiltration and vertical drainage are the dominant hydrological processes in the valley bottom.

Model validation

Pearson’s correlation coefficients between measured and simulated soil matric potentials were acceptably high (when compared to similar studies in the same area – see Van Tol et al. 2015) for most of the soil horizons (Table 5), especially in the crest and the upper-midslope. This suggests that there is a good (positive) correlation between measured and simulated changes in the soil water content. Pearson’s correlation coefficient values obtained in a previous modelling study for the same area for the riparian soil ranged between 0.42 and 0.75 (Van Tol et al. 2015), rendering >0.7 R-values of this exercise acceptable. The lower correlation coefficients in the Sterkspruit soil of the midslope could be attributed to macropore flow dominant between the structure units of the prismacutanic B horizon. This bypass flow is not simulated well with the model. The RMSE values between measured and simulated matric, however, are also high, suggesting that the absolute measured matric potentials are not predicted very well. Because the measurements of soil matric potential are also an estimation, albeit a calibrated estimation, there is a degree of error in these measurements. Future work should focus on direct measurements of soil water content (gravimetric measurements) to validate modelling and to confirm the accuracy of, or recalibrate, the watermark sensors.

The SEBAL-derived aET and simulated aET values correspond relatively well for all landscape units during the first part of the simulations (Figure 3) (Note that the cumulative graphs in Figure 3 are not of consecutive days, as there were gaps in the SEBAL-derived aET values. The corresponding simulated values were also omitted from the statistics). A relative under-simulation of SEBAL-derived aET compared to simulated aET is observed towards the end of the simulation period. The differences could be explained in the dynamic nature of vegetation, which flourishes after wet periods. This is not simulated properly by the model, hence the underestimation. In theory, it implies that parameters such as LAI and root depth should also be dynamic and change with the changes of actual vegetation characteristics during a simulation. For this modelling exercise, the relatively high R and low RMSE values were deemed acceptable for establishing the impacts of changing precipitation patterns in the next section.

FIGURE 3: Cumulative surface energy balance algorithm for land-derived actual evapotranspiration and cumulative-simulated actual evapotranspiration for: (a) the crest, (b) the footslope and (c) the valley bottom of the studied hillslope together with (d) corresponding correlation coefficients and root mean square error of daily actual evapotranspiration (mm day−1).

Precipitation scenarios

A summary of the impacts of the different precipitation scenarios is presented in Table 6. It is important to note that there is non-linearity between different processes and properties. For example, a 30% decrease in rainfall will result in an average decrease of around 12% in the water content of the A horizon of the soils of the crest, but a 30% increase in precipitation will only yield an average increase of around 5% in the water content of the same horizon. This non-linearity is more pronounced in complex processes such as overland flow. Here, a 30% decrease in rainfall will result in a 67% decrease in overland flow volume, whereas a 30% increase will only cause a 17% increase in the volume of overland flow during a rainfall year (Table 6).

Soil water contents of selected soil horizons

Figure 4 presents exceedance curves of soil water contents of various soil horizons under different precipitation scenarios. Sub-soil horizons of both the Clovelly and Dundee soils have higher water contents, for longer, than surface horizons. This is presumably because of the impact of evaporation on surface horizons. In general, the water content is considerably lower for the dry and dry extreme scenarios when compared with average or wet scenarios. The duration curves of the ‘average’ scenario are more closely associated with wet scenarios than dry ones. This could imply that a drying climate will have a greater impact on soil water contents than a wetting climate (see also differences in the average water content presented in Table 6). In the simulations, soil water content affects evapotranspiration; therefore, higher soil water content has a higher potential to produce evapotranspiration. For this reason, the impact of an increase in rainfall patterns on soil water contents is not as clear as the impact of a decrease in rainfall.

FIGURE 4: Duration of exceedance curves of soil water contents of (a) orthic A horizon and (b) yellow-brown apedal B horizon of the Clovelly soil on the crest position; and (c) orthic A horizon and (d) stratified alluvium B horizon of the Dundee soil in the valley bottom.

There is very little difference between the ‘dry’ and ‘dry extreme’, and ‘wet’ and ‘wet extreme’ scenarios in all the simulations. This is likely because of relatively high hydraulic conductivity of all horizons in the Clovelly and Dundee soils (Table 4). An increase in rainfall intensity will therefore have little impact on infiltration in the soils when simulated with a daily time-step.

Simulations indicate that the Dundee soil will be wetter than the Clovelly soil (Figure 4). This despite higher simulated evapotranspiration from the Dundee soil (Figure 3) and similar hydraulic properties of the two soils (Table 4). Overland flow from the midslope to the Dundee soil contributes to higher water contents in this soil.

Seepage line functioning

Under the ‘dry’ and ‘dry extreme’ scenarios, a positive matric pressure (i.e. seepage) was simulated for one day only: a decrease of 85% compared to the ‘average’ (Table 6). The ‘wet’ and ‘wet extreme’ scenarios could result in fivefold and sevenfold increase, respectively, in the duration of saturation at seepage line. The lateral flux volume in the profile under the ‘dry’ and ‘dry extreme’ simulations will decrease by 45% and 31%, respectively, whilst the ‘wet’ year will increase lateral flux by 66% and a ‘wet extreme’ scenario by 122%. Lower fluxes associated with a drying climate could result in the accumulation of salts at the seepage line and reduce the spreading of salts over the sodic site. A wetting climate, on the other hand, could result in excessive leaching of salts and nutrients from the slope, and hinder the normal functioning of the sodic site.

Overland flow in midslope

The functioning of the seepage line impacts the soils downslope. Owing to drier soils and lower than average rain events, overland flow is likely to decrease between 65% and 67% for ‘dry extreme’ and ‘dry’ scenarios, respectively. It appears that under dry conditions, the soil water content has a greater influence on overland flow generation than the nature of precipitation events (i.e. only a difference of 2%). This assumption is, however, only applicable to the daily time-step, which was used in these simulations; using a finer time-step (e.g. hourly), the rainfall intensity might be higher than the infiltration rates of the soils and could produce significant overland flow. Under wet conditions, the rainfall intensity has a much stronger influence (18% increase from ‘wet’ to ‘wet extreme’ scenario). Under ‘wet’ and ‘wet extreme’ scenarios, the seepage line produces more return flow (66% and 122%, respectively), and the surface horizons of the soils lower in the hillslope transect are expected to be close to saturation. Large rain events are then more likely to generate significant overland flow when compared to the drier antecedent soil water conditions expected under the dry climate.

Percolation from riparian zone soils

The assumption is that large rain events are more likely to generate significant groundwater recharge than a number of smaller events (Riddell et al. 2020). Our simulations however indicate that the highest volume of percolation from the Dundee soil of the riparian zone will occur during the ‘wet’ scenario (33% increase) (note that percolation in this sense does not necessarily imply recharge into the groundwater but merely exfiltration into the bedrock layer below the soil). The ‘wet extreme’ scenario only resulted in an increase of 17% (Table 6). A possible reason for this is that the simulations based on matrix flow do not take macropore flow into account; the latter is often dominant during extreme events. This is also evident in the percolation duration curves (Figure 5). Percolation from the ‘dry’ and ‘dry extreme’ scenarios will decrease by 49% and 53%, respectively, which will impact groundwater levels negatively.

FIGURE 5: Duration of exceedance curves of percolation (mm day−1) from the lowest node of the Dundee soil in the riparian zone under different precipitation scenarios.

Conclusion

In this article, we applied pedology and soil physics to quantify the hydrological response of a hillslope under different rainfall scenarios. The quantification of hydrological processes based on this hydropedological approach is important for an improved understanding of the functioning of the hydrological cycle and contribution towards better resource management in the savanna ecosystems. The CMF model was relatively accurate in predicting soil matric potentials and aET. The model structure could, however, be improved by including a dual porosity module as well as time-dependent vegetation parameters (e.g. LAI and rooting depth).

The findings presented suggest that hillslope processes in savanna ecosystems are especially vulnerable to a drying climate. It also highlights non-linearity in the response of hillslope hydrological processes to changes in rainfall. This non-linearity is caused not only by non-linearity in the soil-water retention function of soils but also by the connectivity of landscape. For example, the seepage line would remain saturated seven times longer, with a 30% increase in rainfall because of increased lateral flow from upslope positions. Similarly, drainage from the riparian zone decreases by more than 50%, with a 30% decrease in rainfall because of lower overland flow volumes received from upslope. It is important that resource managers and policymakers in the savanna ecosystems understand and quantify this non-linearity and the influence of hydrological connectivity at different scales. Higher or lower rainfall cannot be simplistically translated as higher or lower recharge, streamflow, runoff, soil water storage and seepage.

Future work should investigate whether these hillslope-scale processes could be extrapolated to catchment and landscape scale. Different processes (e.g. groundwater contributions and return flow from streams) are associated with higher-order streams, which influence the sensitivity of hydrological processes to changing rainfall regimes. It would also be important to establish the impact of changes in hydrological processes on other ecosystem components through interdisciplinary research.

Acknowledgements

The authors declare that they have no financial or personal relationships that may have inappropriately influenced them in writing this article.

J.v.T. was responsible for modelling and compiling the first draft of the manuscript. D.B. and E.S.R. conducted the fieldwork and improved the first draft. S.J. assisted with the modelling and the final write-up of the manuscript.

This work was partially funded by the University of the Free State under the ‘Multi-disciplinary Program’ and The Water Research Commission (WRC) of South Africa under project K5/2051.

Data are available and may be provided, on request, by the corresponding author. Data from all research done within Kruger National Park is placed within the SANParks repository (not for free, open access).

The views and opinions expressed in this article are those of the authors and do not necessarily reflect the official policy or position of any affiliated agency of the authors.

References