Hydrological impact of widespread afforestation in Great Britain using a large ensemble of modelled scenarios

Catchment locations and input data

To determine the impact of afforestation on catchment hydrology we select twelve varied catchments from across the British Isles (Supplementary Material-S1 and Fig. 1). These catchments capture a range of hydrological regimes, drainage patterns and catchment soil and land-cover properties to determine how such factors may influence catchment response to afforestation. Being predominantly >1000 km² in area (ranging from 511 to 9931 km² in size), they are adequately represented in a hydrological model to integrate processes at a 1 km² spatial resolution^54,55. Two catchments are nested within larger ones, the Ure within the Ouse, and the Severn at Bewdley (Severn-B) within the Severn at Haw Bridge (Severn-HB) (Fig. 1).

The period 2000–2010, a flood-rich period for the British Isles^36,37, is chosen to assess afforestation influence on streamflow as it allows us to avoid the uncertainty that would be associated with land-cover changes over a longer period when comparing to baseline results. This length of the simulation period also reduces the computational demand with a large ensemble of land-cover scenarios. Accordingly, the CEH land-cover map for the year 2000⁵⁶, in the form of the CHESS-land dataset⁵⁷, is used to provide configurational datasets specifying soil hydraulic and thermal properties, vegetation characteristics, and orography, for the model at a 1 km² spatial resolution for the unaltered land-cover scenarios. This dataset has successfully been used in other studies^55,58. The 25 m rasterised land-cover map is reclassified into eight different land-cover types (Supplementary Tables S4 and S5) and used to derive afforestation scenarios related to land cover before being converted to a percentage land-cover fraction at a 1 km² spatial resolution. To provide the required meteorological driving data, we use the CHESS-met dataset⁵⁹ which includes long-wave and short-wave radiation, air temperature, specific humidity and pressure. The 50 m CEH Integrated Hydrological Digital Terrain Model elevation data is used to derive topographical and catchment attributes as well as catchment boundaries and river networks⁶⁰. Soil hydraulic information comes from the Harmonised World Soil Database and was made uniform across each grid cell⁶¹.

The model

The Joint UK Land Environment Simulator (JULES) is a physically based land-surface model that simulates the fluxes of carbon, water and energy at the land surface when driven by a time series of the atmospheric data^23,24. Multiple studies have used JULES before including the investigation into evapotranspiration drivers across Great Britain^58,62, atmospheric river formation over Europe⁶³, the impact of solar dimming and carbon dioxide on runoff⁶⁴ and developing river routing algorithms with a Regional Climate Model⁶⁵. JULES is routinely used at the Met Office, where it is coupled with several other models to understand future changes globally and across the UK, by bridging the atmosphere, land surface and ocean⁶⁶. This study is predominantly a theoretical, scenario-based modelling study designed to draw out general principles and to quantify the relation between afforestation and hydrological response, and as such the results are not intended to provide detailed guidance for specific practical actions.

The use of a process-based model enables us to investigate physical explanations for the hydrological impacts of changes in land cover and the explicit representation of vegetation that will influence fluxes, partitioning and storages within the realm of epistemic uncertainty for other conceptual and hydrological models where vegetation is not included. JULES models both plant phenology and canopy storages^23,24. When changing the plant functional type in JULES, both the properties of the above-ground vegetation (such as canopy height and leaf area index) and the soil infiltration factor and the root depth are altered²³. However, there are several caveats that must be considered with this approach. First, the model configuration used in the present study is uncoupled from the atmosphere and so large-scale land-cover changes cannot alter nearby weather⁶⁷. Second, each grid cell is hydrologically separated from adjacent cells, with streamflow and runoff hydrologically uncoupled from the rest of the system. Soil water also does not flow between grid cells. Third, soil thermal and hydraulic properties are uniform across a grid cell. This reduces the impact of hydrological pathways within a cell and the interaction of vegetation with these varying soil types that could have ramifications at multiple temporal and spatial scales. For example, within the cell there may be vegetation that is water-stressed (e.g. valley sides) compared with vegetation where water is not limited (e.g. floodplain) which would change how much transpiration is possible and thus runoff⁶⁸.

Precipitation in the model is partitioned by vegetation and when it reaches the soil surface it is portioned into either infiltration excess overland flow, at a rate controlled by the hydraulic conductivity of the soil, or saturation-excess overland flow as determined by the Probability Distributed Model (PDM)^69,70. Throughfall (T_F) through the canopy is dependent on the rainfall and the existing water in the canopy:

where P is the rainfall rate (kg m⁻² s⁻¹), C is the amount of water in the canopy (mm), C_m is the maximum water storage of the canopy (mm) and ε_r is the fraction of the grid cell occupied by convective precipitation. The maximum amount of canopy water storage is a function of the leaf area index (L):

where A_m is the ponding of water on the soil surface and interception by leafless vegetation (mm) and B_m is the rate of change of water holding capacity with leaf area index. At each timestep (n) the canopy storage is updated thus:

Based on the surface energy balance, the fraction of the proportion of water stored in the canopy compared with the maximum canopy capacity of that plant type is used to calculate the effective surface resistance to determine tile evapotranspiration.

Surface runoff is generated by two processes in JULES: infiltration excess, where the water flux at the surface is greater than the infiltration rate of the soil, and saturated excess overland flow where the water flux at the surface is converted to runoff when the soil is completely saturated. To calculate the saturation-excess overland flow, the PDM⁶⁹ is used to determine the fraction of the model grid cell that will be saturated (f_sat) which is used as a multiplier to convert any excess water reaching the surface to runoff:

where S is the fraction of the grid cell soil water storage, S₀ is the minimum storage at and below which there is no surface saturation (mm), S_max is the maximum grid cell storage (mm) and b is the Clapp and Hornberger⁷¹ soil exponent. We use the topography-derived parameterisation for the b and S₀/S_max parameters to reduce individual calibration with the following relationship⁵⁵:

where s is the grid cell slope (°) and s_max is the maximum grid cell storage (mm). Once interception and surface runoff have been calculated, the remaining water enters the soil. This water is allocated to the different soil layers within the soil column by using the Darcy–Richards equation:

where W is the vertical flux of water through the soil (kg m⁻² s⁻¹), k is soil conductivity (kg m⁻² s⁻¹), φ is suction (m) and z is the vertical flux of water through the soil (m). To calculate suction and soil conductivity, we use the van Genuchten⁷² scheme:

where θ is the average volumetric soil moisture (m³ m⁻³), θ_s is the soil moisture at saturation (m³ m⁻³), α and m are van Genuchten parameters dependent on soil type. The hydraulic conductivity is calculated thus:

where K_h is the hydraulic conductivity (m s⁻¹) and K_hs is the hydraulic conductivity for saturated soil (m s⁻¹). ε is an empirical value set at 0.5 in JULES and S is found by:

where θ_r is the residual soil moisture (m³ m⁻³). Vegetation can access water from each level in the soil column as a function of the root density where the fraction of roots (r) in each soil layer (l) from depth z_l-1 to z_l is:

where z_l is the depth of the lth soil layer, d_r is the root depth (m) and z_t is the total depth of the soil column (m). The water flux extracted from a soil layer is e_lE where E is transpiration (kg m⁻² s⁻¹) and e_l can be found by:

and β_l is defined by:

where θ_c and θ_w are the volumetric soil moisture critical and wilting points respectively (m³ m⁻³) and θ_l is the unfrozen soil moisture at that soil layer (m³ m⁻³). In this configuration of JULES, when a soil layer becomes saturated, the excess water is routed to lower layers. When the bottom layer becomes fully saturated any excess water is added to the subsurface runoff. Both the surface and subsurface runoff are then passed to the River Flow Model^65,73 which routes the flows according to a flow direction grid⁷⁴.

This study uses a combination of calibrated model parameters from the previous work of Robinson et al.⁵⁹ and Martinez-de la Torre et al.⁵⁵ (Rose suites u-bi090 and u-au394, respectively, which can be found using the Rose/Cylc suite control system: https://metomi.github.io/rose/doc/html/index.html). We compare observed streamflow from the NRFA database⁷⁵ with model output for the years 2000–2010 using the base land and CHESS-met datasets. The model is spun-up for the years 1990–2000 to ensure soil moisture content has been equilibrised. To quantify the accuracy of the model, we use a range of standard error metrics. These include the Nash–Sutcliffe efficiency⁷⁶ measure:

Kling–Gupta efficiency⁷⁷:

Root-mean-squared error:

Mean absolute error:

where Q_sim is the simulated discharge, Q_obs is the observed discharge, r is the linear correlation between observation and simulations, σ_sim|obs is the standard deviation of discharge, μ_sim|obs is the mean of discharge and n is the number of observations. We also use NSE(log(Q)) and KGE(1/Q) to understand how well the model can reproduce low flows. Using these measures, we find that JULES performs satisfactorily apart from the Avon Catchment which may be due to fast subsurface flows generated by its geology⁵⁵ (Supplementary Table S7). With process-based models, it is difficult to both accurately reproduce physical processes and make the output faithful to reality due to epistemic uncertainties⁷⁸. Even though model performance is not the same as achieved with calibrated conceptual or empirical models, it allows us to determine the effects of vegetation changes on the hydrological cycle.

JULES’ ability to faithfully represent hydrological land-surface processes in Great Britain has been evaluated in several studies^58,79,80 and the plant functional type parameters it uses at global scales^81,82. To validate the ability of our configuration of JULES to represent soil moisture and potential evapotranspiration rates, we compare the model output with observations from twelve COSMOS-UK sites within our catchments covering grasslands, croplands, coniferous and broadleaf woodland⁸³ (Supplementary Fig. S8). We evaluate model performance from the start of the COSMOS-UK station records until January 1, 2018 so that we use the same forcing data as our experiments. Station start dates vary from October 2013 to August 2017. We compare COSMOS-UK observed soil moisture to the first 0.1 m of the soil column in JULES and evaporation to the sum of the soil evaporation and plant transpiration. We find a median KGE score of 0.44 for the topsoil moisture and 0.53 for potential evaporation (Supplementary Tables S9 and S10). Low error metrics observed for topsoil moisture are due to systematic undercalculation by JULES⁸⁰. At our broadleaf sites, Alice Holt and Wytham Woods, we find both systematic over and underprediction of the topsoil moisture respectively. In line with other studies, we find that there is a slight overestimation of evaporation in JULES^58,84. As illustrated by the median coefficients of determination between the COSMOS-UK and JULES data of 0.62 and 0.60 for the topsoil moisture and evaporation respectively, JULES broadly represents changes in these variables over time.

Land-cover scenarios

Modelling the influence of afforestation on catchment hydrology has been attempted before but usually only at the scale of a single catchment for a limited range of scenarios. In this study, we focus on the theoretical effect of widespread planting of broadleaf trees to examine whether planting location is a stronger control on hydrological response than afforestation extent by using a large ensemble of up to 288 land-cover change scenarios. We choose to focus just on broadleaf woodland for several reasons. First, we are trying to replicate a landscape that could be considered a natural climatic climax community that might occur if it had not been for human intervention during the Holocene. Second, broadleaf woodland has the potential to absorb and store carbon in soils for longer time periods. Finally, to reduce computational cost and the issue of potentially expanding the errors induced by potentially spurious parameters of needleleaf woodland in this version of JULES⁸⁵. Although potential woodland planting locations have been suggested by the Environment Agency and authorities in Scotland and Wales^86,87,88, the differences in planting criteria means it is not possible to systematically compare hydrological changes across our chosen catchments. Here we attempt to create afforestation scenarios related to both catchment river network structure and land use that are directly comparable across a range of catchments. Afforestation was in grassland areas to reduce the complexity of the decisions made and enable an understanding behind catchment sensitivity to land-cover changes related to soil and catchment structure.

Three metrics were selected to discretise the catchment into distinct areas for afforestation: the Topographic Wetness Index (TWI)²⁸, Strahler²⁷ and Shreve orders²⁶. These metrics capture different parts of the catchment such as propensity for saturation, drainage network location and relative contributing areas. TWI is calculated by:

where a is the upslope area draining through a point, per unit contour length, and γ is the local surface topographic slope in radians. All three metrics were calculated using the 50 m IHDTM⁶⁰, thresholding stream formation at an accumulation of ten pixels using the D8 flow direction algorithm within ArcGIS 10.6.1⁸⁹. Strahler order ranges from one (headwaters) to seven (lowlands). Due to the continuous nature of TWI (0.05–31.49) and the large ordinal range of Shreve order (1–9523) calculated for the entire British Isles, we group TWI orders into five quantiles and seven quantiles for Shreve. Increasing TWI order in this case indicates increasing propensity for saturation, or potential maximum saturation level, and increasing Shreve order indicates increasing contributing area. Catchments were broken down to watersheds from the downstream point of the Shreve and Strahler orders. Due to the nature of the data, this led to some first order Strahler catchments being incorrectly generated for some catchments (Supplementary Table S7). Using these generated catchment areas, we plant both inside and outside of these watersheds to understand the hydrological difference between opposing planting locations. In each of the catchment areas, two different levels of afforestation were tested of ~25 and 50% of the possible planting area. Planted area was assigned at random in the catchment and was produced by calculating the area available for afforestation and randomly producing points that covered the area required using the Create Random Points tool in ArcGIS 10.6.1.

Discussions exist about where to plant woodland in relation to existing land cover, to provide ecosystem services, including around watercourses^29,30, urban areas^31,32 and woodland^4,33. Therefore, in this study we try to understand how these potential planting scenarios will affect hydrology in general. Using the CEH 2000 land-cover map⁵⁶ buffers of broadleaf land cover were created at 25 and 50 m around these three land uses (Supplementary Fig. S7). These were then discretised according to the catchment areas. As an example, one scenario would be afforesting up to 50 m around existing broadleaf woodland inside the Shreve order one catchment area, whilst another would be randomly afforesting within 25% of the available area outside of TWI order five areas.

Afforestation according to different catchment areas and land-cover uses between 234 and 288 scenarios for each catchment and between 0 and c. 40 percentage point increase in broadleaf woodland (Supplementary Fig. S2 and Table S8). Due to the structure and size of the different catchments, and thus differences in Strahler and Shreve orders, not all catchments had a comparable number of higher orders. Produced scenarios were converted to the 1-km² grid scale by altering the fraction of land-cover types within each grid cell. It should be noted that this work only considers the impact of mature broadleaf woodland and neglects the influence of the initial planting and growing of the woodland that would likely have its own impact on catchment hydrology as frequently reported^13,49. Furthermore, it does not include the period when there would be the highest amount of carbon sequestration. This study seeks to understand the theoretical impact of woodland on catchment hydrology when fully developed to understand the long-term implications of management decisions.

Hydrological signatures and analysis

Several hydrologic indices can be used to characterise the influence of afforestation on streamflow regime^34,35. To analyse average streamflow and extremes, we look at the top 1% (very high flow), 5% (high flow), 50% (median flow), 90% (low flow) and 95% (very low flow) quantiles of daily streamflow. To quantify flow variability, we use the slope of the flow duration curve^38,40 calculated thus:

where Q_33% is the 33rd flow exceedance quantile and Q_66% is the 66th flow exceedance quantile. To ascertain catchment responsiveness to climatic forcing, we use median streamflow elasticity^40,41:

where dQ and dP are the annual changes in yearly discharge and precipitation, respectively. Finally, we use the runoff ratio to quantify water balance changes related to streamflow and evapotranspiration⁴²:

where µ_Q and µ_P are the average yearly discharge and precipitation using daily values, respectively. We also qualitatively assess the largest peak flow daily event in the 10-year record used in this study to determine the impact of afforestation on the highest possible flows in each catchment.

To determine how afforestation influences streamflow metrics, percentage changes in flow metrics are plotted as a function of percentage point increases in afforestation (calculated using the difference between original and afforested scenario). Quantile regression is applied to determine the median regression slope of the trend for the entire period⁴³. The benefit of using quantile regression is that it identifies the median response of the input variable (in this case the level of afforestation in both percentage and absolute terms) without being influenced by extreme outliers. In this way, we can estimate the proportional streamflow response to afforestation over the period. We use the regression slope coefficient as a proxy of catchment sensitivity to afforestation for each streamflow metric. The slope coefficient is then correlated to catchment attributes, as stated in the CAMELS-GB dataset⁴⁴, using Spearman’s rank correlation. This allows us to determine the direction and significance of the catchment property influences on the sensitivity of catchments to afforestation for the different hydrologic signatures. To determine the impact of different planting locations according to catchment and land-cover location a one-way analysis of variance (ANOVA) test is undertaken using R⁹⁰.

Source: Resources - nature.com