The hydrological bene ts of restoration: A modelling study of alien tree clearing in four mountain catchments in South Africa

Alanna J. Rebelo (  alanna.rebelo@gmail.com ) Stellenbosch University https://orcid.org/0000-0002-7544-9895 Petra B. Holden University of Cape Town https://orcid.org/0000-0002-3047-1407 Jason Hallowes Ekosource Bruce Eady Ekosource James Cullis Zutari Karen J. Esler Stellenbosch University https://orcid.org/0000-0001-6510-727X Mark New University of Cape Town https://orcid.org/0000-0001-6082-8879


Introduction
Ecosystem degradation has undermined the well-being of approximately 3.2 billion people globally, with signi cant economic impacts (Abhilash, 2021). To prevent, halt and reverse such degradation, the United Nations Decade of Ecosystem Restoration, launching in 2021, aims to accelerate and up-scale restoration efforts globally (UNEP & FAO, 2020). Increases in tree cover in certain ecosystems (for example alien tree invasion of grasslands or shrublands, or native bush encroachment) are key global drivers of degradation causing biodiversity and ecosystem service loss as well as economic and human health impacts (IPBES, 2019;Rai & Singh, 2020;Stafford et al., 2017;Stevens et al., 2016). It is estimated that one fth of the planet is at risk due to biotic invaders (IPBES, 2019). Many countries are investing in ecological restoration of ecosystems that are degraded by these trees, for example through alien clearing interventions coupled with either passive or active vegetation recovery (Dimson & Gillespie, 2020;Fenouillas et al., 2021;Holmes et al., 2020;Nuñez et al., 2017;Stafford et al., 2017), or bush thinning (Hare et al., 2020;Smit, 2004).
Restoration, including clearing of invasive alien trees, has also been termed ecological infrastructure investment. Ecological infrastructure is de ned as "the underlying framework of natural elements, ecosystems, and functions and processes that are spatially and temporally connected to supply ecosystem services" (Dominati, 2013). The concept was developed to mainstream the value of biodiversity, or in other words to communicate the bene ts of natural ecosystems to other sectors, in a way that may appeal to planners, engineers and investors (DEA & SANBI, 2016;Pringle et al., 2015). Investments into ecological infrastructure via speci c interventions have been de ned as "arti cial or natural actions that aim to enhance chosen ecosystem services in intact to transformed landscapes, informed by an understanding of ecology" . To upscale these ecological infrastructure investments (i.e. undertake restoration at a larger scale), it is necessary to leverage previously untapped sources of funds, such as from the private sector (Beck et al., 2019;Trinomics & IUCN, 2019). However to access these funding streams, there is a need for a clear evidence-base of the bene ts of investment into ecological infrastructure, and these bene ts are often not quanti ed . Monitoring and evaluation of success, and measurements of the bene ts of previous investments into ecological infrastructure has often been coarse-scale, conceptual or lacking Dimson & Gillespie, 2020;Halme et al., 2013;Kumar et al., 2021;Ntshotsho et al., 2011;Rebelo et al., 2021;Seddon et al., 2020).
The water towers of the world, often mountainous areas, have been so named for the disproportionate amount of run-off they produce globally (Vivirolli et al., 2007). In semi-arid, water-scarce countries, these regions are especially important for water supply and security of cities and irrigated agriculture downstream, especially in the context of a changing climate . Trees invading these water towers constitute a threat to water security, especially in shrubland, savanna and grassland ecosystems. This is because tree species are estimated to have high water-use relative to native vegetation for these ecosystem types (Meijninger & Jarmain, 2014), reducing available water resources Skowno et al., 2019). Previous research modelling the hydrological bene ts of investments into ecological infrastructure has been relatively coarse, at large spatial scales (4 400-5 610km 2 ), using lumped hydrological models, and where there is a general lack of adequate, long-term stream ow data for validation (Hughes et al., 2018b(Hughes et al., , 2018aMander et al., 2017;Rebelo et al., 2015).
Measuring the hydrological bene ts of investments in ecological infrastructure applications at coarse spatial and temporal scales makes it di cult to disentangle the bene ts of investments into ecological infrastructure from other activities in complex catchments (e.g. with high degrees of human in uence, such as groundwater abstraction, irrigation extent and methods, soil degradation (Jiang & Wang, 2019)), or background processes and variations (e.g. re and vegetation succession). Termed scaling-issues (Blöschl & Sivapalan, 1995), these may also mask site-level effects, downstream effects, and effects on downslope processes, particularly at scales that link to local risks and livelihoods (Warburton et al., 2012). One of the twenty-three unsolved problems in hydrology has been framed as how model structural, parameter or input uncertainty can be disentangled and reduced in hydrological prediction (Blöschl et al., 2019). Another that is of applicability here is that of hydrological laws at the catchment scale, and how these change with scale (Blöschl et al., 2019).
Therefore ne-scale hydrological modelling could present an opportunity to isolate the bene ts of investments into ecological infrastructure as well as to understand and compare spatial variability in modelling results at appropriate scales.
This study aims to address these gaps by: (1) bolstering the evidence base of the water-related bene ts of investment into ecological infrastructure using a ne-scale modelling approach (with eld work and remote-sensing inputs, and well-validated models); and (2) explicitly modelling and quantifying uncertainties, for example those linked to spatial variation in geology, using a novel approach which involves the independent development and calibration of three hydrological models for different catchments in water towers within the same region. Therefore, this study aims to address two research questions: 1. What are the hydrological bene ts of investing in ecological infrastructure, and are the results from this ne-scale modelling approach congruent with the international literature? 2. What insights do the detailed exploration of uncertainty yield, and what are the implications for decision-making?

Methods
We modelled the impacts of clearing invasive alien trees on stream ow at a ne-scale in four strategic water providing catchments (Le Maitre, Walsdorff, et al., 2018), or water towers, using the fully distributed MIKE-SHE modelling tool. To maximise certainty and to improve on previous modelling studies, we chose headwater catchments in water tower areas as their strategic upstream position renders them important for regional water supply, there are reduced levels of complexity due to limited land-use in the mountains, as well as due to the strategic location of reasonable quality, long-term observed ow records.

Study Region: Catchment descriptions
Four catchments ranging from 46 to 78 km 2 in size were selected from two strategic water providing catchments in the Western Cape of South Africa ( Figure   1, Table 1), the Upper Berg and Dwars catchments from the Berg Secondary Catchment (8 821 km 2 ), and the Du Toits and the Elands catchments from the Breede Secondary Catchment (12 562 km 2 ). All four catchments were selected because they are high-altitude (elevation ranging from 543.1 to 1 085.2 mamsl, slope ranging from 20.8 to 31.3°), upper/watershed catchments draining into the Berg and Breede Rivers, and are thereby least transformed, allowing for relatively accurate gauging and therefore modelling exercises. These four catchments also form part of critical water towers of the region, supplying the Western Cape Water Supply System which supports the metropolitan city of Cape Town and surrounding agriculture (LaVanchy et al., 2019). These catchments also fall within strategic ground-and surface water source areas for the region . Strategic Water Source Areas are regions within South Africa that provide a disproportionate amount of water resources relative to their size.
The Berg and Breede Catchments are located within the re-adapted Fynbos Biome (a biodiverse shrubland), and have a mild Mediterranean climate characterised by winter rainfall resulting from the passage of mid-latitude cold fronts (Midgley et al., 2003). Soils are mainly nutrient poor (highly leached) dystrophic lithosols associated with the quarzitic sandstones of the Cape Supergroup (Midgley et al., 2003). Mean annual precipitation is estimated to range between 1648.5 to 2553.3 mm for the catchments (Table 1). For all the catchments, with the exception of the Dwars, urban and agricultural development has been minimal. The Upper Berg catchment was impounded in 2009 by the Berg River Dam. There are various transfers between the Berg and Breede Rivers, which affect the Upper Berg and Dwars catchments. In the case of the Dwars catchment, the transfer is upstream of gauging weirs, and therefore there is no suitable gauge for validation, whereas in the Upper Berg, gauging is above the inter-basin transfer. Despite limited development in three of the catchments, all of them have been invaded by alien trees to differing degrees, through riparian infestations (e.g. Acacia mearnsii -Black Wattle) and high-altitude terrestrial infestations (e.g. Pinus spp, Acacia longifolia -Longleaf Wattle), making them suitable candidate catchments for investments into ecological infrastructure through alien plant clearing.

Model Choice and Con guration
A spatially explicit, ne-scale modelling tool was required for this study that could process ne-scale data inputs, and generate results at the grid cell scale, on a daily time-step, both for stream ow but also for other water balance components.
To achieve this, we selected the MIKE-SHE modelling tool. MIKE-SHE is an advanced, exible framework for hydrological modelling (DHI, 2017) and is a fully integrated surface and groundwater model. It may be con gured as either a lumped or distributed model, and has various simple and advanced algorithm options (numeric engines) to represent the different hydrological processes. The more advanced numeric engines are deterministic, physics-based (i.e. solves the partial differential equations describing mass ow and momentum transfer) and distributed model codes. Therefore the parameters used in these equations can be measured and used in the model, however this makes data acquisition intensive, and therefore costly. MIKE-SHE covers the major hydrological processes, including process models for evapotranspiration, overland ow, unsaturated ow, groundwater ow and channel ow, as well as their interactions. To simulate channel ow, MIKE-Hydro River is used and coupled to the MIKE-SHE model. MIKE-Hydro River is able to model complex channel networks, lakes, reservoirs and river structures. Although not well-calibrated for South Africa, MIKE-SHE has been used in a number of related applications globally, including investigations of anthropogenic impacts and land-use/land-cover change (Butts & Graham, 2005).
For this research, we used a fully-distributed version of MIKE-SHE coupled with a MIKE-Hydro channel network. As a trade-off between resolution (the nest resolution would have been 12.5 m digital elevation model) and run-times, we selected a grid size of 60 m for each catchment and simulations covered a period of 15 years (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) at a daily timestep (and 2-hour timestep for overland ow and the unsaturated zone ow), with varying validation periods based on the available observed ow datasets (Table 1). For the evapotranspiration process model, we used the Kirsten and Jensen numeric engine, for overland ow: the fully-distributed nite difference numeric engine, for the soils (unsaturated zone ow) we used the 1D Finite Difference Richards Equation, and for groundwater ow (saturated zone ow) we used the 3D Finite Difference -Darcy Flow numeric engine. For MIKE-Hydro, we used the 'river' model type and the 'hydrodynamic' river module, for a xed time step with a length of 30 to 90 seconds. We used the 1D St Vernant Equations for channel ow (speci cally the fully dynamic equation, high order friction for wave approximation); no ow routing, and the Manning method for water-level. In order to decide how best to con gure the MIKE-SHE model, a conceptual model of each of the catchments was constructed (Figure 2). After a brief description of the conceptual framework used, the con guration of each of the hydrological processes is described in the following sections, after which the model parameters and input data are described.

Catchment Conceptual Model
Using values extracted from literature (see Literature Review and Table S1, Supplementary Material), as well as available geological maps (DWAF, 2007), a conceptual framework for each catchment was constructed, and used to inform model con guration and parametrisation. The Peninsula and Nardouw formations of the Table Mountain Group Aquifer are the major aquifers in the upper Berg and Breede catchments, separated by the Winterhoek Mega-Aquitard. Both the Peninsula and Nardouw formations are composed of highly fractured quarzitic sandstone and are heterogeneous and anisotrophic entities, but are often assumed to be homogeneous and anisotropic at local scales for modelling exercises (Lin, 2007). As most of these rocks have low primary porosities, they only become adequate aquifers where fractured (Pietersen & Parsons, 2002). The hydraulic conductivity of the Table Mountain Group Aquifer is mainly controlled by fractures or fracture networks, and therefore a decrease in the hydraulic conductivity with depth suggests a closing of fractures (Lin, 2007). Therefore hydraulic conductivity is not uniform, but strongly associated with fractures, and very di cult to acquire data for, and accurately represent in hydrological models. In the Upper Berg, in both the Assegaaiboschkloof and Wolwekloof tributaries, there is a strong fracture running along the riverbed, suggesting high connectivity between river ow and groundwater (Lasher, 2011). The thickness of the aquifer is estimated from existing, coarse layer les (DWAF, 2007).
To con gure the MIKE-SHE model in a way to mimic the subsurface piston ow which is characteristic of these catchments (Midgley & Scott, 1991), a separate, perched aquifer was created (Figure 2), which we name the "talus aquifer". Parametrisation with values collected from the literature revealed that the Table Mountain Group Aquifer would never produce the observed recession curve after a rain event (Table 1). With only the unsaturated zone (soils) and Table  Mountain Group Aquifer, a critical part of the hillslope hydrology for the region would be excluded: the intermediate translation of rainfall into stream ow.
Therefore we added the perched aquifer with higher conductivities above the Table Mountain Group Aquifer. This talus aquifer may be considered to be composed of talus material as well as weathered bedrock ( Figure 2). Due to the differing conductivities of the two aquifers, there is a permeability barrier between them, which results in the majority of the water in the perched aquifer moving laterally quite quickly after a rainfall event. Therefore the talus aquifer produces the inter ow portion of stream ow, and the Table Mountain Group Aquifer produces slower moving base ow. The talus aquifer may be transient, only becoming saturated for short durations during the winter. Recharge to the Table Mountain Group Aquifer may come from both the direct exposure to the surface (e.g. at mountain tops and cliffs where the rock is exposed) as well as vertical drainage (percolation) from the perched aquifer. The low conductivities parametrised for the Table Mountain Group Aquifer (either the Peninsula or Nardouw formations), means that it is unnecessary to add more than one aquifer to represent these layers.

Overland Flow
The explicit solver type was used (advised for where over-bank spilling is allowed), as well as the weir-formula to calculate overland-river exchange. For overland ow, uniform values were used for all three parameters over the whole catchment in each case (Table 1), established through calibration (see section on Calibration).

Evapotranspiration
Actual evapotranspiration is calculated for every grid cell on a sub-daily timestep, based on the input land-use/land-cover type, an associated input leaf area index, crop factor, canopy interception factor and root depth. Canopy interception, which is interception and then evaporation by the vegetation canopy, is calculated rst in MIKE-SHE, using a canopy interception factor which is multiplied per computational time-step, which for these four models was on average 1-2 hours. The remaining water reaches the soil surface, either producing runoff, or percolating into the soil. Part of this water is then either evaporated from the soil, from open water bodies, or transpired from the vegetation (based on reference evapotranspiration, the input crop factor and root depth). Finally, the remainder of the in ltrating water recharges the groundwater, where the roots may still extract it directly, if the roots are deep enough, or indirectly through capillary rise. Actual evapotranspiration is calculated using empirically derived equations based on the work of Kristensen & Jensen (1975).

Unsaturated Zone
Based on the model con guration and soil input data, the model was not sensitive to the complexity of the unsaturated zone (see section on Calibration) and therefore a simple set-up was used, with uniform values for the entire catchment in each case. The unsaturated zone was con gured as two layers: a top (0.5 m depth) and a bottom (connection) layer (15 m depth). Using the fully distributed con guration for both unsaturated and saturated zones means that the model pro les for these two zones need to overlap in MKE-SHE to allow the water table to rise and fall over this boundary and simulate capillary rise from an uncon ned aquifer. To allow seamless integration, the lower unsaturated zone was parametrised with the same values as the talus aquifer, rendering the bottom soil layer a purely conceptual layer. This was necessary as a result of the characteristically thin (<0.5 m) soils measured in the eld for these catchments (see section on Fieldwork). Parameter values can be found in Table 1. Lateral ow is not calculated in unsaturated layers, only vertical redistribution. No macropore ow was modelled in the unsaturated zone.

Saturated Zone
The saturated zone was also con gured as two layers (see Catchment Conceptual Framework and Figure 2), with a pre-conditioned conjugate gradient, transient solver-type and no under-relaxation. The connection layer (talus aquifer) was parametrised with higher hydraulic conductivities than the Table  Mountain Group Aquifer (Peninsula Formation) and has the same parameters as the lower soil layer. For parametrisation values, see Table 1. It was assumed that there is little or no ux into or out of the catchments in terms of the regional aquifer, purely because there are no data available on this. However a gradient was set out of the catchment in each case, at the catchment outlet. Initial potential head was set to 0.9 m below the ground for the talus aquifer for the Elands and Du Toits catchments, and was hot-started for the Berg and Dwars. The initial potential head for the Peninsula Aquifer was hot-started for all models. Hot-starting involves running the model for a full cycle (e.g. 14 years), and outputting the values for potential head of the aquifers at the end of this cycle to start the actual model run. This allows the potential head of the aquifers to reach equilibrium.

River and Channel Components
For each catchment, simple river branches were de ned, along which cross sections were automatically inserted, using the input topography, 500 m apart and 50 m wide ( Figure S1). Fieldwork informed the morphology of these cross sections (see section on Fieldwork), but only the relative elevation values were used to de ne the shape of the cross sections, not the absolute elevation values, to avoid scale-issues due to differences in resolution (i.e. eld GPS values compared to a digital elevation model resampled to 60 m). Where we did not have eld-measured cross section data (e.g. between measured points), the width and depth of cross sections were inferred from upstream and downstream measurements, and Google Earth imagery. To improve the connection of the MIKE-SHE and MIKE-Hydro topography and that of the cross sections, topography was lifted along the river channel in MIKE-SHE such that bank minus ground between the topography and cross-sections was never greater than zero, and sinks in the topography due to the resampling of the digital elevation model were lled in MIKE-SHE (sinks can cause ponding which may increase model run times).

Model Parameters and Input Data
The input data for various model parameters for each of the four study catchments are summarised in Table 1, and the details given in four subsequent sections: metrological input data, land-use/land-cover input data, eldwork and labwork, and literature review sections.  Table 1). For an overview of the selected rainfall stations, see Figure S2, Supplementary Material. Where there were periods of missing data (Table S2, Supplementary Material), rainfall records were patched using a linear regression approach, using nearby stations where relationships were very strong, i.e. correlation coe cient greater than 0.8, but mostly greater than 0.9 (Pearson Correlation). The more complex set-up was used for rainfall in MIKE-SHE, whereby results were spatially distributed using the Thiessen polygons, each assigned a rainfall station ( Figure S2, Supplementary Material) for catchments with more than one station (Berg, Dwars), and for where there was only one station driving the catchment, no polygon was used (Du Toits and Elands).
Due to the mountainous nature of the catchments, and the fact that most gauges were from low altitude rainfall stations (with the exception of some of the SAEON gauges), we corrected precipitation for elevation using a lapse rate. To create a lapse rate layer, we used a 12.5 m digital elevation model (ASF DAAC, 2015), a rainfall surface (DWAF, 2007), and the Thiessen polygons previously described for the relevant catchments. Rasters were converted to point les, and all were combined such that each point had an elevation and rainfall value, as well as a catchment number, at the scale of the nest input. Lapse rate was then calculated as follows: Lapse Rate = (cell rainfall/station rainfall -1) / (cell elevation -station elevation) * 10 000 This was done one Thiessen polygon at a time (where relevant), such that the correct station rainfall and elevation were used in each calculation. In some cases, where there were many elevations in the catchment similar to that of the rainfall station, with differing rainfall values, this resulted in extreme lapse rates (either far too high, or far too low), which can result in rainfall that is too high, although MIKE-SHE does not allow negative rainfall. Therefore to stop an overestimation of rainfall, we capped the lapse rates, in various ways, depending on the catchment and the elevation of the rainfall station(s) driving the lapse rate calculation. In one case, the lapse rate was used as is (low elevation station, Elands), in another case, the lapse rate was capped above zero and an additional 3% was added (Du Toits) based on station data near the catchment, and in two cases, the lapse rate was capped by one standard deviation on either side of the mean (Berg, Dwars) ( Table 1). The resultant point le was then converted to a raster for import into MIKE-SHE. The lapse rates are then converted into rainfall correction factors internally within MIKE-SHE using the following formula: Correction Factor = Lapse Rate * (cell elevation -station elevation) The corrected precipitation for each cell is then calculated within MIKE-SHE as follows: Cell rainfall = station rainfall * (1 + Correction Factor) MIKE-SHE requires Penman-Monteith reference vegetation potential evapotranspiration as well as rainfall as a climatic input data. We used ARC reference evapotranspiration data from automatic weather stations, which was modelled using the FAO Penman-Monteith method (Zotarelli et al., 2016). To try and input distributed evapotranspiration, we considered METREF data ( Figure S3, Supplementary Material), however the ARC reference evapotranspiration had much greater annual amplitude, and therefore using METREF resulted in an underestimation of actual evapotranspiration for key land-use/land-cover types, and this approach was therefore discarded. However the daily time-series of the ARC reference evapotranspiration data had inaccuracies and did not cover the full modelling time period. This resulted in actual evapotranspiration rates for key land-use/land-cover types also being slightly underestimated. Therefore we developed a method to model the average shape of the daily ARC time series using this dataset, considering amplitude and wavelength, to remove the daily noise, that could be extrapolated to cover the full period. To achieve this, we tted a sine curve in r, to predict the daily ARC time series based on the daily reference evapotranspiration input data ( Figure S3, Appendix 1, Supplementary Material). We then exported the predicted values at a daily timestep and input this evapotranspiration at a catchment scale for each of the four catchments. Only two reference evapotranspiration datasets were used for the four catchments: reference evapotranspiration from the Berg catchment (ARC Station 30890), and from the Elands catchment (ARC 30942_1) (Table 1).

Land-use/Land-cover Data
For all four catchments, the native vegetation is fynbos shrublands with pockets of Afromontane Forest in the re-protected ravines (Table S3, Supplementary Material). This native vegetation has been invaded by alien trees to varying degrees . We used a newly generated spatial database of 14 land-use/land-cover classes for the Upper Berg and Breede Catchments, derived using an SVM classi cation of a Sentinel 2 image at 20 m spatial resolution and combined with the 2014 Western Cape Province Land-use/Land-Cover Map for urban and agricultural areas, and a map of res overlain to determine vegetation age   (Smithers & Schulze, 2004). The speci c ACRU parameter is vegetation interception or VEGINT (the potential interception by vegetation input as mm/day on a month-by-month basis), which was then adjusted to a 2-hour timestep (mm/2 hour). The back-calculation is necessary because ACRU vegetation interception is already inclusive of leaf area index, whereas in MIKE-SHE the vegetation interception gets multiplied by leaf area index within the model. The adjustment from a daily to two-hourly value is necessary because MIKE-SHE calculates the vegetation interception at every timestep, which is sub-daily and in our case on average two-hourly. Therefore if this adjustment is not done, the canopy interception is overestimated, and not enough water reaches the soil surface for in ltration and transpiration processes. Vegetation factors (crop factors) were extracted from the ACRU COMPOVEG and FAO56 (FAO, 1998) databases and averaged, except for wetlands where crop factors between these databases and  were averaged. Due to the fact that ACRU crop factors are for A-Pan evaporation, they were corrected for reference evapotranspiration by multiplying by a factor of 1.2 (Schulze, 1995).

Fieldwork and Labwork
During a two-week eldtrip in November 2019, we measured leaf area index, soil properties and cross sections in each of the four catchments in the Berg and Breede Catchments, South Africa. Leaf area index was measured in the eld using a LI-COR LAI-2000 plant canopy analyser for each of the main landuse/land-cover categories used in this modelling study for which leaf area index data are scarce (Alien trees: Black Wattle, Gum, Pine, Other as well as low and high density fynbos, burnt fynbos, indigenous forest and wetland). For each class, at least 20 readings were taken from different points of representative patches under cloudy conditions. In addition to the eld-scale measurements we calculated leaf area indices from Sentinel-2 satellite imagery for pure pixels of the same classes, but also including irrigated and dryland agriculture, to test the ability to upscale. These two independent sets of results were compared with a literature review (Table S5, Table S4, Figure S4, Supplementary Material), showing reasonable congruence, and from these, the nal leaf area index inputs for the model were selected ( density was measured using a Kopecky Ring of known volume, soil depth was measured using an auger (cored to bedrock), and from the soil that was cored, composite samples were taken of each soil layer that emerged (where there was more than one). All soil samples were stored in plastic bags, and bulk density samples in plastic containers for processing.
Bulk density samples were weighed before and after oven drying for 48 hours at 105°C and soil moisture and bulk density calculated, and values for the latter expressed as g/cm 3 . For each soil sample, a soil particle size analysis was performed using a hydrometer and sieve according to the classi cation of DoAD (1991) (Gee & Bauder, 1986). Determination of water retention measurement by controlled out ow cell for 0-1 bar potentials using the provided bulk density was also performed for each sample (Lorentz et al., 1991). Saturated Hydraulic Conductivity was measured by packing to speci ed dry bulk density and conducting hydraulic conductivity testing (Lorentz et al., 2001). The data are available in Table S6, Supplementary Material. Soil particle size, water retention and saturated hydraulic conductivity laboratory work was done by the Soil and Water Laboratory at the Centre for Water Resources Research (CWRR) of the University of KwaZulu-Natal in Pietermaritzburg, South Africa.
Finally, cross sections of each of the main river channels for each of the four catchments were measured ( Figure S1, Supplementary Material). We used a Differential Global Positioning System (DGPS) to map out the shape of each cross section. Data were processed and then manually input into MIKE-HYDRO (see River and Channel Components for more details).

Literature Review
The literature was reviewed for (a) aquifer parameter values (Table S1, Supplementary Material), (b) vegetation leaf area index (Table S4, Table S4

Calibration
In the three catchments with adequate observed ow data, the models were calibrated for the dates available (Table 1). For the Dwars catchment, the Berg calibration parameters were used for the nal model, given similar geologies and close proximities (Blöschl, 2005). Initial exploratory calibration showed that the models were insensitive to soil data (unsaturated zone), and vegetation properties, but extremely sensitive to parametrisation of the saturated zone and slightly sensitive to overland ow parameters. Therefore the focus was initially on calibrating the saturated zone, followed by a more minor calibration for overland ow parameters. Due to the many parameters and the long-run times of a fully distributed MIKE-SHE coupled MIKE-Hydro models, the automatic calibration tool of MIKE-SHE, AUTOCAL, was used. We used the parameter optimisation function for the four parameters of the saturated zone (horizontal hydraulic conductivity, vertical hydraulic conductivity, speci c yield, speci c storage) in a rst phase, followed by one parameter of overland ow (Manning Number) in the second phase. Upper and lower bounds, and initial values for the calibration were taken from the literature (Table S1, Supplementary Material). For the calculation of comparison statistics only one function type was used for the objective function (the weighted sum of absolute values), and thus no transformations were done as there were no aggregations. The statistic type used to compare with the observed ow was RMSE. For parameter optimisation, we used the population simplex evolution with Monte Carlo sampling using initial parameter values. Population size was set to 10, with two points in the simplex (DHI, 2017). Other parameters for parameter optimisation were left at the default (e.g. three loops to convergence). We set the number of simultaneous simulations to four, based on our number of processors, and used a single thread engine.
For the unsaturated zone, we rst attempted a simple set-up of one uniform soil over the whole catchment based on average values. We then tried incorporating more spatial variation, but it had no signi cant effect in improving results, and in some cases signi cantly increased model run times. For the saturated zone, we found that the highest accuracies (relative to observed ow) were obtained when the Peninsula Aquifer is contributing below 0.01% to discharge (Figure 2). This is also congruent with the literature values for the Peninsula Aquifer, which is characterised by very low conductivity values (Table   S1, Supplementary Material). Although there were some variations between the catchments, parametrising the talus aquifer with much higher hydraulic conductivities than the Peninsula Aquifer worked best, thus representing the inter ow in the catchment. For overland ow, we initially tried a range of distributed options (e.g. changing Manning Number for wetlands compared to slopes), however this reduced model performance (from ~80% to ~50%) and therefore we reverted to uniform values for overland ow.

Validation
Following satisfactory calibration, validation was performed for each of the catchments with adequate observed stream ow data (i.e. excluding the Dwars catchment). The relevant periods included in the analysis were based on available observed stream ow data, always allowing at least one year for model equilibration.

Ecological Infrastructure Scenarios
Five ecological infrastructure scenarios were considered for each of the catchments in this study (Figure 3, Table 2), including a baseline scenario (scenario 1), the current state of the catchment (scenario 2), and three theoretical scenarios (scenarios 3-5). The baseline scenario (scenario 1) was considered to be a reference condition for the catchments, such that all land-use remains the same (e.g. any agriculture, urban areas, water impoundments), however the invasive alien trees are all cleared to a maintenance level (i.e. no more adult trees), representing the best (realistic) state of ecological infrastructure for each catchment. Impacts on stream ow were compared to the baseline scenario for subsequent analyses. Results could not be compared to the current state (scenario 2), due to differing levels of invasion in each of the catchments. The current state of ecological infrastructure in the catchments (scenario 2) as of 2019 was based on classi cation of Sentinel-2 imagery at a spatial resolution of 20 m , representing varying levels of alien tree infestation and ages.
Scenarios 3 and 4 are theoretical scenarios of degraded ecological infrastructure based on scenario 1 (i.e. again being realistic by not replacing current agriculture, urban areas and impoundments with invasive alien trees) that are the inverse of each other: scenario 3 has terrestrial infestations, and scenario 4 has infestations in the riparian zones. Riparian zones were de ned using the Advanced Land Observation Satellite (ALOS) dataset, which has 15 different landform classes in the dataset including variations of: peak/ridge, cliff, upper slope, lower slope, and valley classes. We combined the valley classes to de ne the riparian zones. The dominant terrestrial invasions in the four catchments are pine infestations of the mountains, though there are also other invasions of gums and wattles, and therefore to standardise stream ow reduction results, for terrestrial invasions, mature pines (7-15 years old) at 100% density were used (

Analyses
To answer the research questions, ve analyses were performed.

Exceedance Flow Probabilities
Flow duration curves were constructed for each catchment and every scenario using the hydroTSM package in r, and displayed on a log scale to improve low ow separation (Zambrano-Bigiarini, 2020). The percentage of time each of the stream ow magnitudes given as input was equalled or exceeded were extracted for six percentiles: 5th, 10th, 25th, 50th, 75th and 95th using the fasstr package in r (Goetz et al., 2020).

Stream ow Reduction
For each catchment and every scenario, the stream ow reduction was calculated by taking the difference of the stream ow for each scenario relative to that of the baseline (scenario 1) and expressing it as a percentage of that of the baseline at either a seasonal or annual timestep. These results were used as input into the seasonality analysis, as well as the spatial analysis (see below). In addition, to situate the modelling results within the theoretical understanding of the stream ow reduction curves of Scott & Smith (1997), an additional modelling simulation was done, whereby the total invadable area (i.e. excluding water impoundments, urban areas, roads and agriculture) of the catchment was replaced by pine (7-15 years old at 100% density) and the mean annual reduction in stream ow calculated relative to the baseline (scenario 1). These points were then overlain on the new stream ow reduction curves from the South African paired catchment data, which incorporates uncertainty estimates (Moncrieff et al., 2021).

Seasonality
To investigate the relationship between stream ow reduction in wet relative to dry years, we used the annual stream ow reductions for each of the scenarios, and compared these to the total precipitation for each year for each of the four catchments using Pearson Correlations in R. To investigate seasonal variations in stream ow reduction, we used the seasonal stream ow reductions and plotted Box-and-Whisker plots in R for each of the scenarios (including all catchments), and for each of the catchments (including all scenarios), as well as one plot for all catchments and scenarios together.

Spatial Analysis: Riparian versus Terrestrial Impacts
To gain a theoretical understanding of modelled stream ow reduction impacts in different parts of the landscape, the stream ow reduction impacts from scenario 3 (terrestrial invasion) and scenario 4 (riparian invasion) were expressed per unit of area (% per km 2 ). In addition to stream ow reduction (%), the change in mean stream ow (m 3 .s −1 .km −2 ), mean annual volume (Mm 3 .km −2 ) and mean annual runoff (mm.km −2 ) were also calculated. The ratio of stream ow reduction in terrestrial relative to riparian areas was also calculated for each catchment and the overall mean was also calculated. The stream ow reduction (%) for terrestrial compared to riparian areas was also graphed.

Water Balance Extractions
MIKE-SHE explicitly models multiple ow pathways at the scale of grid cells, where overland ow moving between grid cells can in ltrate, inter ow can percolate into an aquifer, and water tables can rise and saturate soils, producing surface ow, amongst others. For these ow pathways, the MIKE-SHE water balance extraction records the pathway by which water entered the river, regardless of whether this was the dominant pathway through the landscape. So for example, water may have travelled overland for 90% of its pathlength, however if it in ltrates metres away from the river channel, it will be classi ed as inter ow. Likewise with inter ow surfacing meters from the river (overland ow), or inter ow percolating into an aquifer just before entering the river (base ow).
The water balance results can be extracted from MIKE-SHE from a ow result catalogue le (.sheres) for any date range, time-step, or area (based on de ned catchments or catchments), and multiple post-processings can be performed on the extraction. We extracted results at a daily timestep (incrementally), from 2005 to 2018 (excluding 2004 to allow the model one year to equilibrate), for the entire catchment area, for each of the scenarios for each of the four catchments. The mean annual values were calculated in each case for the 14-year period. Six major water balance categories were extracted, including: precipitation, evapotranspiration, overland ow, base ow, out ow and storage (Table S9, Supplementary Material). For evapotranspiration, six minor water balance components were extracted (total evapotranspiration, canopy interception, soil evaporation, transpiration, groundwater evapotranspiration and open water evaporation), and two each for boundary out ow (overland and subsurface boundary out ow), and base ow (inter ow and deep groundwater ow).
Inter ow is the portion of water reaching the river channel from the connection layer (talus aquifer), and the deep groundwater ow is the same for the Peninsula Aquifer. For the MIKE-SHE water balance type and item name, see Table S9 in the Supplementary Material.

Results
The results are presented in three sections, rstly validation results, secondly the hydrological bene ts of investing in ecological infrastructure (based on the scenario results) and thirdly theoretical insights into spatial variation in gains and water balance partitioning.

Validation
Model performance was satisfactory in terms of performance measures and evaluation criteria (Moriasi et al., 2015), with less than 9% difference in means (r=0.77-0.79), and similar minimum and maximum values (Table 3). PBIAS is low (between -10 and 2), indicating low percentage deviation from the mean, and in the case of the Upper Berg and the Du Toits catchments, a slight underestimation of rainfall, whereas for the Elands catchment, a slight over-estimation.
Nash-Sutcliffe e ciency is acceptable for a daily timestep (NSE=0.52-0.59), whereas the Nash-Sutcliffe e ciency of the logged discharge suggests that prediction of lower ows is better (NSE=0.74-0.81) (Moriasi et al., 2015). There is slight underestimation of peak ows, and slight overestimation of base ows overall (Figure 4, Figure S5, Figure S6). Table 3 Daily stream ow statistics (m 3 /s) for the simulated compared to the observed for the baseline scenario (scenario 2) for the Upper Berg and Elands catchments and an adjusted scenario the Du Toits catchment to re ect the appropriate historical validation period of the three case studies with accurate gauges, Western Cape, South Africa. For the Dwars River comparison to two catchment gauges, please see The hydrological bene ts of investing in ecological infrastructure The different catchments produce different amounts of discharge in their current state (Scenario 2), ranging from as high as 4.2 m 3 /s for the Upper Berg catchment, to as low as 1.1 m 3 /s for the Du Toits catchment, due to differences in area, mean annual precipitation and levels of alien tree infestation. This translates into volumes ranging from 131.9 Mm 3 to 34.1 Mm 3 respectively (Table 4). The models predicted mean annual runoff (MAR) values ranging from 305 mm (Du Toits -Scenario 5) to 1544 mm (Upper Berg -Scenario 1), with scenario impacts on stream ow ranging from 0.4 (Elands -Scenario 2) to 29.5% (Du Toits -Scenario 5) (Table 4, Figure 5). The increase in stream ow (%) that is predicted for clearing the current levels of infestations (Scenario 2) for all four study catchments ranges from 0.4 to 8.8% (Figure 6). If all catchments were theoretically infested with pine trees in all terrestrial areas (Scenario 3), stream ow bene ts of clearing would range from 12.7-25.3%, whereas if riparian zones were theoretically infested with Black Wattle (Scenario 4), bene ts of clearing would range from 1.5-4.8%. Clearing a theoretically fully infested catchment (Scenario 5), would yield stream ow increases from 15.1-29.5%. This translates to changes in volume of about 1.5 and 2.5 Mm 3 , and in mean annual runoff of 18 and 32 mm respectively.
The largest impacts on exceedance probabilities are in the mid to low ows for all four catchments (Figure 7, Table 4). The mean effect of all scenarios relative to the baseline (scenario 1) across all case studies was 25% on the 5th percentile (lowest ows), and 12% on the 95th percentile (highest ows) (Table   S10, Supplementary Material). Therefore the impacts on stream ow are proportionally greater on mid to low ows compared to high ows.

Impacts of ecological infrastructure investments in wet compared to dry years
For all four catchments, there is a negative relationship between the percentage difference in stream ow of all scenarios relative to the baseline and total annual rainfall (Figure 8). This was the most marked for the Du Toits and Elands catchments. Therefore the impact of the reduction in stream ow is exacerbated in dry years relative to wet years.

Seasonal impacts of ecological infrastructure investments
In summer and autumn, stream ow reductions were more pronounced relative to winter and spring, for all scenarios and all catchments (Figure 9). This trend held over all scenarios ( Figure S8) and all catchments ( Figure S9), although the trend was most pronounced for Scenario 4 (the riparian invasion).

Stream ow reduction curves
Using the same models, but simulating stream ow reductions for total invadable area by pine 7-15 years old at 100% density, yields annual stream ow reductions of 15.4 to 30.1% for the 14-year period from 2005 to 2018, which is within the 50% predictive interval of Moncrieff et al., (2021), and well aligned to the initial curve tted by Scott and Smith 1997 ( Figure 10).

Theoretical insights
Impacts in terrestrial relative to riparian invasions Terrestrial invasions ranged in impact from 0.2-0.6% per km 2 on water resources, compared to 1.7 times greater impact by riparian invasions which ranged from 0.4-1.1% per km 2 depending on the catchment (Figure 11).

Water balance partitioning
The mechanism by which invasive alien trees are decreasing stream ow is primarily transpiration according to all the simulations. Transpiration is the component of the water balance that has the greatest impact on discharge between scenarios, ranging from 33-51% of rainfall for Scenario 1, and 43-65% for Scenario 5 (Figure 12, Figure 13). Base ow also declines with alien tree invasion. According to the way the models were parametrised, in terms of base ow, most of the contribution is predicted to be from the talus aquifer (99.99%), therefore the inter ow component, with very little contribution from the Peninsula Aquifer (<0.001%) (  (Table S14, Supplementary Material), which lies between the range of possible evapotranspiration for mixed terrestrial/riparian fynbos of 520 to 1460 mm/a according to the literature (Table 5).

Discussion
The hydrological bene ts of investing in ecological infrastructure In this section we consolidate our ndings on the hydrological bene ts of investing in ecological infrastructure, speci cally the effects of tree clearing, and discuss whether these ne-scale modelling results are congruent with the international literature for the hydrological impacts of alien tree invasion, bush encroachment, and ecological infrastructure investments in general. We end with a discussion around the feasibility of investment to upscale ecological restoration.

Comparison with hydrological impacts of alien tree invasions and bush encroachment globally
This study demonstrated that clearing a catchment fully infested with invasive alien trees such as pines, increased available water resources from 15.1-29.5% (translating to increases in mean annual runoff of 18-32 mm and volumes of 1.5-2.5 Mm 3 ). The percentages of annual ow augmentation are comparable with estimates from another national-scale study in South Africa, suggesting 7.71% for the entire Berg catchment, and 8.74% for the entire Breede catchment based on estimated infestation levels . Another modelling study of the Kromme River, Eastern Cape of South Africa, found that wetland rehabilitation, including clearing of riparian Black Wattle infestations, could yield water resource augmentation of up to 30%, which is also comparable to our results (Rebelo et al., 2015). According to the theory of limits to evaporation, alien tree invasions may be considered analogous to plantations of similar species (Calder, 1998). Using data from some of the longest and most detailed paired catchment experiments, a series of curves were generated to predict impacts of afforestation on stream ow (Moncrieff et al., 2021;Scott & Smith, 1997). The mean annual stream ow reduction modelled by this study is comparable to that of the Scott & Smith (1997) curves for sub-optimal conditions (i.e. in colder areas).
A review of paired catchment experiments from around the world speci cally considered the impact of afforestation experiments on stream ow (e.g. conversions of shorter vegetation to tree cover) (Brown et al., 2005). For three experiments with afforestation with pine or eucalypt in India, New Zealand and South Africa, the average percentage reduction in monthly stream ow ranged from 20 to 60% depending on the season, which compares well to these results. For New Zealand, with reasonably constant rainfall, this resulted in a constant reduction in stream ow. However both India and South Africa were examples of seasonal rainfall, resulting in high seasonal variation in reductions in yield (Brown et al., 2005). There is evidence of exotic trees decreasing available water resources in countries across the world, for example afforestation of grassland regions of Argentina with plantations has raised concerns of salinization of the soils due to their proportionally higher water-use (Jackson et al., 2005). A systematic review demonstrated that exotic tree plantations in high elevation Andean grasslands may reduce water yield by up to 40% (Bonnesoeur et al., 2019), which also compares well with our ndings. A review of the relationship between forest management and water resources found that for certain climates (e.g. in Australia, Brazil, Chile, China, India, South Africa, Spain and the USA), where grassland or shrublands are replaced with exotic plantations, stream ow and groundwater recharge are signi cantly reduced (Garcia-Chevesich et al., 2017). Invasions by alien trees have been likened to the issue of bush encroachment in terms of potential hydrological impacts. A study on woody encroachment in Cerrado vegetation in Brazil found that 0.9% less rainfall reached the soil with every increase in 1m 2 /ha of tree basal area, which has implications for soil saturation and groundwater recharge (Honda & Durigan, 2016). Conversely a review on the impacts of bush encroachment on groundwater recharge found highly variable results (Acharya et al., 2018).
In this study, the largest impacts of clearing alien trees on exceedance probabilities are in the mid to low ows, and the impact of the augmentation of stream ow is greater in dry years relative to wet years. Likewise, in summer and autumn, stream ow augmentation was more pronounced than in winter and spring. In a global review of paired-catchment experiments, this trend was observed for both summer-and winter-rainfall regions, however in some cases for summer-rainfall regions similar changes are observed over all months (Brown et al., 2005).

Comparison to studies of ecological infrastructure investments across the world
Other studies considering the bene ts of investments into ecological infrastructure internationally have found that interventions can provide signi cant cobene ts. A study on the rehabilitation of the Yitong River in northern China found that riparian vegetation buffering, and river channel enhancement signi cantly reduced non-point source pollution and improved water quality, as well as provided biodiversity and economic bene ts (Mi et al., 2015). An international review of nature-based solutions, a synonym for ecological infrastructure interventions, of over 1700 documents from science and practice concluded that co-bene ts in terms of the environment, social and economic dimensions are apparent (Raymond et al., 2017). Similarly, a review of European nature-based solutions for urban water management found that particular value was derived in terms of drought and ood protection, the water-energy-food nexus, and water puri cation, with co-bene ts for biodiversity, society and urban microclimate (Volkan Oral et al., 2020).
Nature-based solutions also have strong potential for climate change adaptation as well as mitigation, but major barriers include evidence-based implementation, as well as nancial and governance challenges . It is essential that in the age of big data, cloud processing, cloud-based computing and the fourth industrial revolution, water management keeps pace with such developments to ensure more accurate and judicious planning in the face of climate change. Creative solutions are especially necessary for data and resource scarce nations that suffer from water scarcity. For example, in the desert megacity of Lima, Peru, which has almost 10 million inhabitants and less than 10 mm of rainfall per year, a transdisciplinary approach to water management was developed (Schütze et al., 2019). This approach combines adaptation and methods from hydrology, social sciences, water engineering and modelling, as well as input from stakeholders, which encourages ownership and acceptance of solutions (Schütze et al., 2019).

Feasibility of making an investment case to upscale restoration efforts
Other hydrological modelling studies to date, despite technical limitations, have demonstrated that the potential for returns looks promising, with the proposition that protecting and rehabilitating ecological infrastructure could generate meaningful gains in water quantity (Mander et al., 2017). Furthermore, the costs of interventions are comparable with those of traditional built infrastructure solutions, suggesting economic feasibility (Mander et al., 2017). This study, using a ne-scale, physical-based model, nds that the water-related bene ts of clearing mature infestations of alien trees, such as pines, from shrublands, are increased available water resources of 15.1-29.5%. This study also con rms that gains in water resources are particularly during the dry season and drought years, which has been found in many other paired-catchment experiments globally, both in winter-and summer-rainfall regions (Brown et al., 2005). This is an important nding as it suggests that investments into ecological infrastructure could help to improve resilience, by increasing availability of water resources particularly when they are needed the most. This is also signi cant from a water resource management perspective, as for any catchments that are impounded, augmentation of summer ows is likely to be entirely captured, whereas increases in winter ows are likely to be lost through dam over ow, and could also perpetuate ood impacts downstream. Therefore this investment in ecological infrastructure through alien tree clearing is likely to have water regulation bene ts not limited only to improving water supply. Considering the co-bene ts that investments into ecological infrastructure are thought to produce (Pagano et al., 2019), this makes a compelling case for upscaling investment into restoration.
Theoretical insights into spatial variation in gains and water balance partitioning In this section, we discuss the theoretical ndings of this research, as well as insights yielded by the detailed exploration of uncertainty.
The impact of clearing terrestrial compared to riparian invasions Clearing terrestrial invasions ranged in impact from 0.2-0.6% per km 2 on water resources, compared to 1.7 times greater impact when clearing riparian invasions, which ranged from 0.4-1.1% per km 2 . One of the twenty-three unsolved problems in hydrology is what the hydrologic laws are at the catchment scale and how they change with scale (Blöschl et al., 2019). This factor of 1.7 could be tested elsewhere, with different models and at different scales to determine whether it is a consistent hydrologic law. Other studies have found that riparian invasions by wattle, eucalypt, poplar and willows (high water-users) reduce stream ow by a factor of 2 in grasslands and savannas where there is dry-season dormancy, and 1.5 for other biomes with evergreen vegetation (Dye & Jarmain, 2004;Le Maitre et al., 2016). A study on the long-term impact of wattles in grasslands showed that the relative impact in riparian compared to terrestrial areas was 16 mm and 78 mm for areas of 7.5 ha and 65 ha respectively (Everson et al., 2014). This results in a 1.78 factor different between riparian and terrestrial areas, which is congruent with our ndings.
Mechanism underlying changes in water resources following restoration The mechanisms for the changes in stream ow according to the MIKE-SHE model are decreased evapotranspiration following restoration (speci cally the transpiration component thereof) and increased overland ow, with slightly increased inter ow (the portion of the base ow moving through the talus aquifer with shorter residence times). Inter ow was modelled to range from 3.7-9.1% for the four catchments compared to 38.5-61.9% overland ow. It should be clari ed that in these catchments, overland ow is rare, except when severe res have damaged the soils Scott, 1993). Technically there would not be much overland ow after a rainfall event in these catchments but a rapid inter ow response through the soil and regolith following preferential ow paths, more generally known as 'quick ow' . According to a study done in the area using stable isotopes, less than 5% of the storm ow comprised direct runoff (or overland ow) (Midgley & Scott, 1991). Since we know that the MIKE-SHE water balance extraction records only the pathway by which water entered the river, regardless of whether this was the dominant pathway through the landscape, it is possible that a portion of this 'overland ow' is moving mainly below the ground (inter ow), but daylighting at the valley bottom due to saturation, being recorded as 'overland ow' before it enters the stream. Therefore drawing conclusions should be done with caution, but it is likely that restoration is increasing inter ow and base ow given what we know about overland ow in these catchments .
The reason for the decreased evapotranspiration with restoration is because the invasive alien trees have higher evapotranspiration rates than the indigenous vegetation that replaces them (Meijninger & Jarmain, 2014). Since this modelling study is focussed on change in vegetation according to ve ecological infrastructure scenarios, it is critical to ensure that the internal actual evapotranspiration modelled by MIKE-SHE is corroborated by that in the literature.
Terrestrial fynbos evapotranspiration is on average 701 mm/a, compared to 928 mm/a of terrestrial invasive alien trees (difference of 227 mm/a) (Dzikiti et al., 2014;Meijninger & Jarmain, 2014). For riparian systems, riparian fynbos has a mean actual evapotranspiration rate of about 1017 mm/a compared to that of invasive alien trees of 1279 mm/a (difference of 262 mm/a) (Dye & Jarmain, 2004;Dzikiti et al., 2014Dzikiti et al., , 2016. The mean difference in terrestrial and riparian water use is remarkably similar between fynbos and invasive alien trees, with 316 mm/a and 351 mm/a respectively. The question remains whether this is a hydrologic law at the catchment scale, and whether this changes with scale, or whether this is simply a function of the MIKE-SHE model structure (Blöschl et al., 2019). The actual evapotranspiration predicted for fynbos (both riparian and terrestrial) in this study was 821 mm/a (731-866 mm/a) and for invasive alien trees (both riparian and terrestrial), was 1078 mm/a (1031-1121 mm/a). The difference was predicted to be 256.48 mm/a on average (202-323 mm/a). These modelled actual evapotranspiration means and ranges are well supported by the literature, suggesting that the results of this modelling study are reliable. In terms of partitioning, the models predicted that on average 313 mm/a (29%) of the actual evapotranspiration was due to canopy interception and evaporation from the soil, which is slightly higher than (but still comparable to) the 235 mm (23%) estimated from an adapted MOD16 model of a riparian gum (Eucalyptus camaldulensis) infestation in the Berg River, for an annual actual evapotranspiration of 1039 mm/a  .
Insights yielded by the detailed exploration of uncertainty A point of interest is the variation in the water related bene ts among the four water tower catchments, even after standardising for different catchment discharges (expressing bene ts as a percentage change relative to the baseline scenario in each case) as well as standardising per unit area of infestation (which accounts for the effect of different catchment sizes and percentages of infestations). This may be explained to some extent by the indigenous vegetation replacing the alien tree infestations. For example, high-density fynbos and indigenous forest have higher actual evapotranspiration compared to low-density fynbos, due to differences in characteristics such as biomass and root depth. However there must be something else over-riding this effect, given that the Du Toits catchment has the greatest bene ts in terms of increase in water resources following restoration, and yet has the highest coverage of highdensity fynbos in the baseline scenario (scenario 1). This should result in a smaller proportion of gains, as the difference between alien and indigenous vegetation water use is smaller.
Similarly, the unsaturated zone (soils) and overland ow parameters cannot account for these differences, since they are very similar between catchments.
This leaves only one factor, and this is the in uence of the saturated zone. Though we know that there are limitations to how the saturated zone is modelled in this study (i.e. fractures with higher hydraulic conductivities and connectivities cannot be captured), it does suggest that underlying geology plays a role in determining the scope of the bene ts of ecological infrastructure investments. In particular, the differences between the Berg and Breede catchments are marked, the Breede apparently yielding higher returns than the Berg. This is a signi cant nding that improves our understanding of uncertainties in hydrological modelling of the bene ts of ecological infrastructure investments. In understanding impacts of invasive alien trees on ow reduction, not only tree type, density and age is important (Moncrieff et al., 2021;Scott & Smith, 1997), but also underlying geological characteristics. Likewise, impacts on inter ow augmentation seem to be different for different saturated zone parameters, ranging from 0.4-0.7% for the Berg catchments, and 1-1.7% for the Breede catchments (almost 2.5 times greater). Other studies have also found that the relative contribution of ecological restoration is strongly linked to both geology and vegetation changes (Lian et al., 2020).
This raises the question of validity of the saturated zone parametrisation and conceptualisation, given that it has important implications for the results. The main question is whether the Peninsula Aquifer really plays such a minor role in these catchment hydrological processes (i.e. contributing practically no base ow). If this is the case, what are the implications for management and planning, and the proposed groundwater abstraction from the system? If it is not the case, and the Peninsula Aquifer is contributing signi cantly to base ow, will it ever be possible to adequately represent a fractured aquifer system in distributed models like MIKE-SHE, in data scarce regions such as these? In a study of the Berg River, hydrograph separation analysis indicated that the Berg River is 38% dependent on subsurface water discharges annually (Madlala et al., 2018). This study estimated, on average, 5.4% of inter ow and base ow combined for all four catchments, with overland ow playing a more important role overall (46.2%). This may be in part because the water balance extraction tool of MIKE-SHE does not consider the dominant pathway through the landscape, but the pathway by which water enters the river. If this were the case, it would require that water be resurfacing through the unsaturated zone, above the ground, becoming overland ow close to the river. The upper Berg catchment has been suggested to be a groundwater recharge area, suggesting that groundwater is recharged during the high ow season (winter) with discharge taking place during the summer (Madlala, 2015). In addition, it has been suggested that the upper Berg catchment is characterised by relatively short residence times of aquifer water based on NaCl ion analysis (Madlala et al., 2018). This contrasts to the piston-ow hypothesis for these mountain systems proposed by Midgley & Scott, (1991).
The Table Mountain Group Aquifer is known for generally impeding layers of impermeable geological heterogeneity which result in mainly horizontal movements of water, primarily through the unsaturated zone of the aquifer, until reaching an area of connectivity, for example contact zones of springs or seeps (Ratcliffe, 2007). It is also known from sparse resistivity studies that areas of high conductivity appear to correspond to areas of deep fractures in the Berg catchment (Lasher, 2011). There is a paucity of data on where these fractures are exactly in the Berg and Breede catchments, and the exact conductivities of these fractures, as well as the seasonal and inter-annual variation in recharge/discharge dynamics. Despite this knowledge gap, the City of Cape Town and the national Department of Water and Sanitation in South Africa have embarked on a groundwater abstraction project to augment water supplies (City of Cape Town, 2019). The lack of available information and data on these aquifer systems is a major limitation to this study and for water resource management in the region. It is critical that South Africa, and other data scarce countries, invest in long-term empirical research to better understand critical water resources and the effect that ecosystem degradation as well as climate change may have on future water security (Blöschl et al., 2019).
Regardless of the conceptualisation of the saturated zone for this study, it has improved our understanding of the possible uncertainties around ecological infrastructure investments. In addition, since these models are validated against observed stream ow, it does provide some boundaries for understanding the possible contribution of aquifers to local water resources, and the results seem to suggest that groundwater may play less of a role in producing stream ow than previously thought. This could have critical implications for the uncontrolled groundwater abstraction in the region, and may raise red ags around the sustainability of abstractions. Further investigation is needed.

Recommendations
The ndings of this study have some implications for decision-making, therefore we make the following recommendations: Over the next decade, the UN Decade on Ecosystem Restoration, to build resilience in the face of a changing climate, investments in ecological infrastructure through clearing trees in grassland, savanna or shrubland ecosystems at scale (either exotic trees, or thinning bush encroachment) could yield signi cant water-related returns, especially in the dry season for climates with a strong seasonal signature, and dry (drought) years when water is needed the most. We recommend that countries prioritise investments into ecological infrastructure at scale (either through funding or incentivising private sector investment).
We advance the science of hydrological modelling of ecological infrastructure investments by providing more accurate (well-validated, detailed physicallybased model with eld and up-to-date remote-sensing inputs) and ner-scale results. There are many hydrological modelling studies, but a dearth of empirical hydrological monitoring studies of the bene ts of ecological infrastructure investments. Funders and governments should invest in long-term empirical studies that will improve our understanding of ecological infrastructure and the bene ts of investment. This includes improvements in basic monitoring, such as rainfall measurements in high altitude locations, and building an understanding of the dynamics of various aquifer systems.
Our research raises questions around what level of detail is needed in hydrological modelling efforts to capture the necessary uncertainty to reliably and transparently inform integrated water resource management and policy. Can integrated water resource management really take place where there is a limited understanding of rainfall or aquifer dynamics? Our hydrological modelling framework, although still with limitations and uncertainty, is a signi cant improvement on the modelling tools currently used to inform water resource management in data scarce regions. We recommend that in this age of big data and the fourth industrial revolution, water management in government keeps pace with these developments to ensure more accurate and judicious planning in the face of climate change.

Conclusion
In conclusion, this ne-scale, fully distributed hydrological modelling study on four catchments suggests that investments into ecological infrastructure yield signi cant returns in terms of water supply. These bene ts are most apparent in the dry season and in dry years, suggesting that such investments improve resilience as well as provide other co-bene ts such as improved water regulation. Furthermore, this study advances our understanding on one aspect of hydrological modelling uncertainty in a speci c use-case, generating understanding of the general applicability of the principles discovered. Due to increased resilience, bene ts arising from investments into ecological infrastructure are likely to accrue under anthropogenic climate change in the future. These ndings improve the evidence base of the bene ts of investing in ecological infrastructure, and may help to unlock private sector investments, needed to upscale interventions during the United Nations Decade on Ecosystem Restoration.
The location of the four study catchments in the upper reaches of the Berg and Breede Rivers, Western Cape, South Africa. The catchment boundaries are shown in orange, and the black lines on the inset map are provincial boundaries. Large rivers are overlain in dark blue. Service credit layer for base imagery: Drakenstein, Earthstar Geographics.

Figure 2
Map of simpli ed geology for all four catchments, with a map of South Africa (inset), and the conceptual framework of the hydrological ows of the two tributaries (Assegaaiboschkloof and Wolwekloof) forming the Upper Berg catchment, South Africa. The thick black lines on the geology map indicate the location of the cross sections forming the conceptual frameworks, looking upstream (west). The rough proportions of water using each of the hydrological pathways is indicated as percentages (dark blue as precipitation, green as evapotranspiration, light blue as overland ow, cream as inter ow).

Figure 3
Five scenarios of alien tree invasion modelled for the four case studies, Western Cape, South Africa.

Figure 4
Mean daily (left) and logged daily (right) simulated (red) and observed (black) stream ow (m 3 /s) for the baseline scenario (scenario 2) for three of the four case studies, the (a) Upper Berg, (b) Du Toits, and (c) Elands catchments, Western Cape, South Africa. Rainfall (blue) is displayed on s secondary axis. For the Dwars River comparison to two catchment gauges, please see Figure S7. For a two year demonstration period, please see Figure S5, and for monthly timeseries please see Figure S6, Supplementary Material.

Figure 5
Daily modelled stream ow output for the ve ecological infrastructure scenarios as predicted for the four case studies, the (a) Upper Berg, (b) Dwars, (c) Du Toits, and (d) Elands catchments, Western Cape, South Africa. The full time series is shown (left) as well as a two-year demonstration period (right).

Figure 6
The percentage change in stream ow (mean stream ow -m 3 .s -1 , mean annual volume -Mm 3 , mean annual runoff -mm) for the four ecological infrastructure scenarios relative to the baseline (pristine -Scenario 1) as predicted for the four case studies, Western Cape, South Africa. Catchment values are overlaid on the plot: blue = Berg, red = Dwars, green = Du Toits, and yellow = Elands.

Figure 10
Mean annual stream ow reduction from the four case studies overlain on the stream ow reduction curves from the South African paired catchment data. Mean predictions are represented by solid lines, and the curves tted by Scott and Smith (1997) for optimal and suboptimal growing conditions are represented by dashed lines. Dark shaded areas indicate 50% predictive intervals, while lightly shaded areas indicate the 95% predictive intervals (Adapted from Moncrieff et al., 2021). Data are from a 14-year period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) and are for the catchments fully invaded by pine (aged 7-15 years) for all current invadable area (i.e. excluding water impoundments, urban areas, roads and agriculture).