Study catchments and field data

The study catchments comprised three highly instrumented landscapes within the UK Defra Demonstration Test Catchments (DTC) Programme33: Newby Beck at Newby, Eden, Cumbria (54.59° N, 2.62° W; 12.5 km2); Blackwater at Park Farm, Wensum, Norfolk (52.78° N, 1.15° E; 19.7 km2); Wylye at Brixton Deverill, Avon, Hampshire (51.16° N, 2.19° W; 50.2 km2). Meteorological, flow and nutrient data are available at hourly resolution for parts of the period October 2011–September 2014. Further details of these catchments and the monitoring are available in Outram et al.34.

Hydrological modelling

We applied two hydrological models of different complexity to the three DTC catchments. Parameters for both models were identified on all, or part, of the period 1 October 2011–30 September 2013, and validated on the period 1 October 2013–30 September 2014, using observed field data.

The Hydrological Predictions for the Environment (HYPE) model8 is a hydrological model for simulating water flow and transport and turnover of nitrogen and phosphorus. The model is semi-distributed, dividing the landscape into classes according to soil type, land use and altitude. In agricultural lands the soil is divided into up to three layers, each with associated parameters. Soil, water and nutrient processes are simulated, with surface runoff, macropore flow, tile drainage and outflow calculated from individual soil layers. HYPE runs at a daily timestep, although the daily flows and phosphorus loads were calculated from the sub-daily observed data.

Nutrient inputs into the model in the form of fertiliser and manure were determined from the Defra British Survey of Fertiliser Practice. Nutrient inputs from point sources were added to each reach of the river or stream based on the measured discharges and pollutant concentrations from the Environment Agency national register of consented discharges for the period 2010–2012, as used in a national source apportionment tool35. No changes were made to point sources for the future.

Parameters in the HYPE model were chosen using the Generalised Likelihood Uncertainty Estimation (GLUE) methodology36 which samples the multi-dimensional parameter space to find sets of parameters which produce acceptable models which satisfy the evaluation criteria (termed behavioural models). While we recognise the sampling size is not extensive, given other experimental uncertainties we consider the chosen trajectories are indicative of the potential for change and in balance with all other computational and experimental design constraints. Models were evaluated using the Nash Sutcliffe Efficiency (NSE) coefficient37 for discharge (Q) and for P. The thresholds for behavioural parameter sets were chosen to give the best 10–15 behavioural parameter sets. For Eden these were: NSE ≥ 0.6 for Q and NSE ≥ 0.5 for P; for Wensum these were NSE ≥ 0.55 for Q and NSE ≥ 0.53 for P; for Avon these were NSE ≥ 0.6 for Q and NSE ≥ 0.6 for P. The ranges for parameter sets used in the projections are given in Supplementary Data 3. The HYPE model calibration and validation fit statistics are given in Supplementary Data 4. HYPE allows simple user-specified changes to land use or management by variation of the input parameters.

Data-Based Mechanistic (DBM) modelling38, using the CAPTAIN Toolbox for MATLAB39 identified transfer function models for rainfall-runoff and rainfall-phosphorus load directly from the high temporal resolution (hourly) data, requiring very few parameters. Either discrete-time or continuous-time transfer function models9, with the structures and parameters given in Supplementary Table 1, were identified directly from the hourly resolution observation data.

A second-order discrete-time linear transfer function with no noise model, denoted by [2, 2, δ] takes the form:

$$y\left( t \right) = \frac{{{b_1}{\rm{ + }}{b_2}{z^{ - 1}}}}{{1 + {a_1}{z^{ - 1}}{\rm{ + }}{a_2}{z^{ - 2}}}}u\left( {t - \delta } \right)$$ (1)

where y(t) is model output at time t, u(t) is model input, z −1 is the backwards step operator, i.e., z −1 y(t) = y(t−1). b 1 , b 2 , a 1 , a 2 are parameters determined during model identification and δ is the number of time steps of pure time delay. For a physical interpretation, models are only accepted if they can be decomposed by partial fraction expansion into two first-order transfer functions with structure [1, 1, δ] representing fast and slow pathways, with characteristic time constants and steady state gains, i.e.

$$y\left( t \right) = \frac{{{b_{\rm{f}}}}}{{1 - {a_{\rm{f}}}{z^{ - 1}}}}u\left( {t - \delta } \right) + \frac{{{b_{\rm{s}}}}}{{1 - {a_{\rm{s}}}{z^{ - 1}}}}u\left( {t - \delta } \right)$$ (2)

where b f and b s are gains on the fast and slow pathways, respectively, and a f and a s are parameters characterising the time constants of the fast and slow pathways respectively. a f and a s are roots of the denominator polynomial in the second-order transfer functions above (Eq. (1)).

A second-order continuous-time linear transfer function with no noise model takes the form:

$$Y\left( s \right) = \frac{{{b_1}s{\rm{ + }}{b_2}}}{{{s^2}{\rm{ + }}{a_1}s{\rm{ + }}{a_2}}}{{\rm e}^{ - s\tau }}U\left( s \right)$$ (3)

where, Y(s) and U(s) represent the Laplace transforms of the output and input, respectively. b 1 , b 2 , a 1 , a 2 are parameters in the denominator and numerator polynomials in the derivative operator \(s = \frac{\rm d}{{{\rm d}t}}\) that define the relationship between the input and the output, and τ represents the delay. Models are only accepted if they can be decomposed by partial fraction expansion into two parallel, first-order transfer functions, i.e.

$$Y = \frac{{{b_{\rm{f}}}}}{{s + {a_{\rm{f}}}}}{{\rm e}^{ - s\tau }}U{\rm{ + }}\frac{{{b_{\rm{s}}}}}{{s + {a_{\rm{s}}}}}{{\rm e}^{ - s\tau }}U$$ (4)

where a f and a s are direct reciprocals of the fast and slow time constants respectively, which define the fast and slow components of the response. b f and b s are parameters which determine the gain of the fast and slow components, respectively.

Note that parameters b 1 , b 2 , a 1 , a 2 (and parameters b f , b s , a f , a s ) have different interpretation, and therefore different values between discrete-time and continuous-time models. The relationship between the parameters (see most Control Engineering textbooks40) between discrete model denoted by superscript d and continuous time model denoted by superscript c is as follows:

for instance, for denominator parameter a f

$$a_{\rm{f}}^{\rm{d}} = {{\rm e}^{ - a_{\rm{f}}^{\rm{c}}{\rm{\Delta }}t}}$$ (5)

while for b f we have:

$$b_{\rm{f}}^{\rm{d}} = \frac{{b_{\rm{f}}^{\rm{c}}}}{{a_{\rm{f}}^{\rm{c}}}}\left( {1 - {{\rm e}^{ - a_{\rm{f}}^{\rm{c}}\Delta t}}} \right)$$ (6)

Models are evaluated according to the Nash Sutcliffe Efficiency (NSE)37 and the Young Information Criterion (YIC)38, an objective statistical measure which combines how well the model fits the data together with a measure of over-parameterisation.

Linear models, using observed rainfall, were identified first, and then improved where possible by use of the ‘effective’ rainfall. The rainfall-runoff non-linearity41 was based on the storage state of the catchment, for which the discharge was used as a proxy, i.e.

$$R{\rm{e}}\left( t \right) = R\left( t \right){\left( {Q\left( {t - 1} \right)} \right)^\beta }$$ (7)

where Re(t) is the effective rainfall at time t, R is the observed rainfall, Q is the observed discharge and β is a constant exponent which is optimised from the observed data at the same time as model identification. For rainfall-phosphorus load models, linear models using observed rainfall were identified first, and then improved using the same effective rainfall relationship as identified with the rainfall-runoff model, with the effective rainfall generated one step ahead, using Eq. (7) and the simulated discharge at time (t−1).

The DBM models are able to make best use of the high-frequency data to capture P dynamics which typically occur at time scales of hours rather than days. High temporal resolution measurements of nutrient dynamics have previously demonstrated that a daily time step is insufficient to capture sediment and P dynamics42, resulting in an underestimation of export loads. Definitions, structure, and parameters for the models identified are provided in Supplementary Table 1, along with model fit statistics for the identification and validation periods (Supplementary Data 5).

Both the water quality model HYPE and the simple DBM model assume that relationships and processes identified on the basis of present day data will still hold in the future. However, the HYPE model, despite the uncertainty associated with the large parameter set, allows variation of the parameters to simulate changing environmental and management conditions, whereas the simple DBM model reduces parameter uncertainty but has no explicit way to change the identified parameters. We have used the DBM models to investigate the P response to climate change only.

Future climate data

We drive the models with future climatic rainfall and meteorological data from a convection-permitting (1.5 km grid) regional climate model (RCM-1.5 km), which is a configuration of the Met Office Unified Model43, and with data from the UKCP09 Weather Generator11, at hourly resolution for DBM and daily resolution for HYPE. The differences between UKCP09-WG and RCM-1.5 km are given in Supplementary Table 2

The UKCP09 Weather Generator11 creates synthetic time series of weather variables at daily and hourly frequency, at 5 by 5 km grid square resolution, which are consistent with the underlying UKCP09 Climate Projections44 (from an ensemble of 11 RCMs, 25 km grid). The Weather Generator applies stochastic models to generate many different, but statistically equivalent, time series which are stationary for a given time slice (30-year baseline or future time period) and emissions scenario.

The Weather Generator uses a stochastic rainfall model, with other variables generated according to the rainfall state, based on empirical relationships between climate variables in a baseline dataset of observations (1961–1995, 5 km grid). Future weather generator time series are generated by perturbing with change factors taken from the probabilistic projections of UKCP09. This ensures that the overall statistics (means and standard deviations) of the weather generator distributions are the same as those projected by UKCP09. The generated climate distributions include both the natural climate variability (including the significant spatial signature as observed in the 5 km grid baseline observations) and some uncertainty from the UKCP09 projections (as change factors are derived from the 11 member ensemble), giving a statistically-based distribution of climate at each chosen location, for each emissions scenario.

This study employed random number seeds to generate 100 plausible time series of 30 years (30 year baseline period and 30 year future time period), for each location, for each specified emissions scenario, as detailed below. The parameters used for generation of the UKCP09-WG time series data are given in Supplementary Table 3 (Newby Beck, Eden), Supplementary Table 4 (Blackwater, Wensum) and Supplementary Table 5 (Wylye, Avon).

RCM-1.5 km is much more computationally demanding and, therefore, we use results from two 12-year simulations: one for the period 1996−2009 and a 12-year simulation representative of the year 2100, under the RCP 8.5 scenario. The high resolution, convective permitting model allows more realistic representation of rainfall at a local scale, particularly for extreme events13.

Expert elicitation to guide future agricultural changes

We determined likely future agricultural changes through expert elicitation with stakeholders at workshops held in each catchment. Discussion on identical topics in each catchment was followed by completion of a questionnaire (Supplementary Tables 6 and 7). The responses to agricultural change options (selected responses in Supplementary Fig. 1) guided the modelling of future changes. Much of the discussion of likely agricultural changes was indicative of intensification of agriculture, either through an increase in livestock numbers or the need to grow more crops on the same land. The results of the discussions and the questionnaires were used to define simple land management change scenarios that could be easily incorporated into the HYPE model. Due to the semi-distributed nature of the model, it was not possible to make spatially explicit land use changes. Using lumped land use changes (not related to location, but related to the hydrological response unit) increased uncertainty in the result, as there was large variability in the result from different spatial representations. As the discussions and questionnaires pointed to intensification in farming, it was decided to include this in as simple and transparent a way as possible. The percentage changes for agricultural P inputs (±20, ±50, ±80% relative to baseline) are arbitrary, but cover a range of different scenarios for intensification. This simplified representation of likely agricultural changes could be applied transparently to all three catchments without the spatial uncertainty associated with more specific measures.

Projections under climate change only

We make future projections of discharge and P load under climate change only using both HYPE and DBM with the same parameter sets as identified with the observed field data. Using climate data from UKCP09-WG, annual loads (for baseline and scenario conditions) were calculated by taking the average over the last 26 years of each 30 year run, with the first 4 years being used as a spin-up period. For each 30 year run with HYPE, the averages were also taken over all behavioural runs. Winter loads were calculated by summing over the months December, January and February for each run. Summer loads were summed over June, July and August. 5th and 95th percentiles were calculated from the spread over 100 runs. Using climate data from RCM-1.5 km (12 year run), annual and seasonal loads were calculated by taking the mean over years 3–12.

Projections under combined climate and agricultural change

We made future projections of P load under combined climate change and agricultural change scenarios using HYPE. The agricultural scenarios, indicating a change in intensification of agriculture, were represented by increasing or decreasing the fertilizer and manure inputs (−80, −50, −20%, no change, +20, +50, +80%). No changes were made to point source inputs. Annual and seasonal loads were calculated as above.

Data Availability

The DTC data are available at http://www.environmentdata.org/dtc-archive-project/dtc-archive-project (Browse by Collection>Defra Collections>Eden/Hampshire Avon/Wensum Demonstration Test Catchment Data). UKCP09 data and access to the Weather Generator is available at http://ukclimateprojections.metoffice.gov.uk/.

The data underlying the figures in this manuscript are openly available from Lancaster University data archive at https://dx.doi.org/10.17635/lancaster/researchdata/110.