Hydrologically driven ecosystem processes determine the distribution and persistence of ecosystem-specialist predators under climate change

Climate change has the capacity to alter physical and biological ecosystem processes, jeopardizing the survival of associated species. This is a particular concern in cool, wet northern peatlands that could experience warmer, drier conditions. Here we show that climate, ecosystem processes and food chains combine to influence the population performance of species in British blanket bogs. Our peatland process model accurately predicts water-table depth, which predicts abundance of craneflies (keystone invertebrates), which in turn predicts observed abundances and population persistence of three ecosystem-specialist bird species that feed on craneflies during the breeding season. Climate change projections suggest that falling water tables could cause 56–81% declines in cranefly abundance and, hence, 15–51% reductions in the abundances of these birds by 2051–2080. We conclude that physical (precipitation, temperature and topography), biophysical (evapotranspiration and desiccation of invertebrates) and ecological (food chains) processes combine to determine the distributions and survival of ecosystem-specialist predators.


Supplementary
. Results from GAMs describing observed bird abundance as a function of modelled cranefly abundance, including 2-dimensional tensor product smooth fitted to x and y coordinates to account for spatial structure in the data. For intercept and cranefly abundance term,

The dynamic water table model
Water-table variation is driven by input from rainfall, losses from runoff and evapotranspiration and a small, constant drainage out of the system to reflect loss into the bedrock. To drive these processes, each modelled peatland is assigned total monthly rainfall (mm) and mean monthly temperature (°C). Each column within the peatland is also assigned values for slope (°), aspect (°) and elevation (m).
Topography is used to modify weather inputs to provide locally-adjusted values for each peat column. Annual rainfall increases with elevation by 2.25 mm m -1 . Although this relationship varies from 0 mm m -1 to 4.5 mm m -1 around the UK 4 , this value represents a compromise to make the relationship widely-applicable. Equation 1 describes local rainfall adjustment: Here, R L = local rainfall, R O = observed rainfall, and E = relative elevation, i.e., the difference in metres between the modelled location and the elevation that weather data refer to. Therefore, a negative relative elevation decreases local rainfall, while a positive relative elevation increases local rainfall. For the monthly model, the 2.25 mm m -1 change is divided by 12, so that when summed over the year, the annual change is 2.25 mm m -1 .
Temperature decreases with elevation by 0.006°C m -1 , which should be applicable throughout the UK 5 . Temperature is also adjusted by slope and aspect to account for variation in incident radiation 6 . The equation used in the annual model 1 is retained for the monthly model: Water input occurs via precipitation; water is lost through evapotranspiration and runoff.
Evapotranspiration is calculated using the Thornthwaite equation 7 as this only requires precipitation and temperature inputs, and can be used for humid and wetland areas 8 . Potential evapotranspiration (PET) estimates are adjusted based on the WTD and vegetation composition to give actual evapotranspiration (AET); AET decreases as WTD becomes deeper 9,10 . The relationship between the AET:PET ratio and WTD changes between plant functional types (PFTs), and is estimated based on maximum rooting depths, root profile distribution and the shape of published relationships 11,12 : declines in AET occur sooner and faster in PFTs with shorter roots. For Sphagnum and other bryophytes, minimum AET is 50% of PET 12 ; for all other PFTs, minimum AET is 80% of PET 13 . AET:PET relationships are described in the table below. The relationship between PFT proportions and WTD is the same as in the annual model 1 .
Summary of relationships between AET:PET ratio and WTD for different plant functional types used in the model. Water inputs and outputs are summed to give a change value. As in the annual model, an exponential relationship is assumed between distance to the water table and available pore space, such that available space increases with distance from the water table. Total space is calculated by integrating over the available unsaturated peat cohorts. By combining the water entering the system with the available space, a new WTD is calculated.

Monthly runoff equations
The model does not use formal hydrological runoff process functions, but instead uses equations that aim to reproduce dominant water-table behaviour observed in blanket peat. Runoff varies from 90% to <10% of rainfall 14 , so equations were developed to reflect this. Runoff is strongly influenced by the existing water table, so runoff is a function of the previous time step's WTD; based on published data 14 , runoff was assumed to be related to WTD exponentially, but with runoff higher and more sensitive to WTD when the water table is within 5 cm of the surface. Runoff also increases with slope angle 15,16 , so a cos function was used to increase runoff as slope steepness increases. A condition was set such that runoff could never be higher than total precipitation. section Runoff equation parameterisation and sensitivity. As runoff data from situations with standing water were not available, the fitted model was set to produce a minimum of 95% runoff, but not to vary strongly with the depth of standing water.

Stabilising modelled WTD behaviour
After runoff is calculated, total water input for the time step is calculated as precipitation -(evapotranspiration + runoff). However, individual months with relatively high or low precipitation could cause unstable behaviour, with rapid flooding and droughts occurring more easily and frequently than observed in real blanket bogs 14,17 . Therefore, after water input was calculated, functions were applied to stabilise model behaviour.
In months of relatively high precipitation, runoff and evapotranspiration were not always sufficient to account for incoming precipitation, causing available pore space to be rapidly used up, in turn causing rapid and long-lasting floods. This probably reflected the inability of the model to represent processes such as interception of rainfall by vegetation, increased surface runoff during high intensity rain events, or snowfall during winter months. Hence, for relatively wet sites, in months where precipitation was greater than average, a correction to water input was applied, shown in Here, W C is the corrected water input, W U is the uncorrected water input (i.e., precipitation -(evapotranspiration + runoff)), R ratio is the ratio of observed monthly rainfall to that year's average monthly rainfall (i.e., annual rainfall/12), and μ is the power-law correction factor. The correction was only used when R ratio > 1 and for sites with average annual rainfall (R av_ann ) > 1500 mm.
Initially derived with μ = 1, it was found that wetter sites were more susceptible to flooding, so μ varied with R av_ann : 1500 ≤ R av_ann < 1625, μ = 1; 1625 ≤ R av_ann < 1750, μ = 2; 1750 ≤ R av_ann < 2000, μ = 3; 2000 ≤ R av_ann , μ = 4. The correction operates such that as the disparity between that month's precipitation and the mean monthly precipitation increases, water input becomes smaller, representing processes such as changes in surface runoff with rainfall intensity. The correction is applied to water input rather than runoff, because the model erosion term, which is calculated as a function of runoff, becomes too large if runoff is dramatically increased. Hence, the correction ensures that in a heavy rainfall event, areas with moderately high WTDs do not experience rapidlyoccurring, long-lasting floods, but areas with lower WTDs still experience re-wetting.
Further, whilst rapid, short-term droughts under very low rainfall conditions were well-replicated, under moderately low rainfall, relatively wet sites could switch rapidly to drought conditions. This was likely either because evapotranspiration or runoff was too high under these conditions. Therefore, a further stabilisation factor was added, shown in Equation 7: 0.5 0.75; * 0.25 0.5; * 1 Variable definitions are as for Equation 6. The correction reached a maximum when R ratio = 0.5, and then decreased again as R ratio approached 0.25 or 0.75; a graduated correction was used to avoid sharp transitions from uncorrected to corrected water input. Hence, the correction allowed the model to continue representing droughts in very dry conditions well, but in months with only moderately low rainfall, water tables fell more slowly. Both corrections were derived by trialling the model for the evaluation datasets (see below) and observing resulting behaviour of model predictions; correction factors were accepted when flood and drought behaviour stabilised across all evaluation datasets. A simple sensitivity analysis was carried out, where each value was raised and lowered by 10%.

Runoff equation parameterisation and sensitivity
The range in R 2 was 0.15, but this was primarily caused by ε being lowered; excluding ε-10%, the range was only 0.06. The range in mean WTD was 2.6 cm; this variation appeared to be caused by ε-10% and κ+10%, and when these were excluded, the range was only 1.5 cm. The most sensitive value was minimum WTD, which had a range of 5.2 cm, but this was again caused by ε-10%; when this was excluded, the range was only 2.3 cm. Maximum WTD was relatively insensitive, with a range of only 1.6 cm. The mean absolute error relative to observed data was 0.6 cm for mean WTD, 0.5 cm for minimum WTD and 1.9 cm for maximum WTD, indicating relatively robust performance. Qualitatively, there was little effect on model predictions from parameter value variation; only ε-10% caused substantially different behaviour, allowing floods to occur more frequently.
The strength of the slope effect (parameter λ could not be parameterised using the Moor House data, as data were not available over a range of slopes. To set the value, modelled peat depth (which, via erosion, is ultimately driven by runoff) was compared to observed peat depth data across shallow, moderate and steep slopes at Lake Vyrnwy (A. Heinemeyer, unpublished data). The slope parameterisation was accepted when it broadly reproduced observed patterns. The use of cos(4.5*slope) meant that no slopes over 40° could be modelled with this parameterisation.
Parameters used in monthly runoff equations and fitted values.

Parameter Equation Description
Value Evaluating model performance The model was run for three British blanket bogs for which observed WTD data were available (see Fig. 1 in main text). The sites were Moor House (northern England), Lake Vyrnwy (mid Wales) and Oughtershaw Moss (northern England); summaries of datasets are presented in the table below.
Observed water tables were converted to monthly means for use in evaluation.
Sites differed in condition and monitoring method so datasets may vary in suitability for use in evaluation. As the model is parameterised for an intact peatland, it may not represent the hydrology at drained sites well. At Lake Vyrnwy, manual dipwell data were not frequent enough to represent true monthly means, so may not be well-predicted. Oughtershaw Moss data allow effects of peatland drainage to be explicitly examined, but as data span only 18 months, longer-term behaviours cannot be evaluated. Consequently, results should be interpreted in the context of site condition and monitoring method.
To generate model predictions, the model was driven with monthly climate data. Monthly predictions started in 1914, at the start of the UK Met Office 5 km gridded climate data 19 . All models were driven by these data until 2010, apart from the Moor House run, which was driven by data from a local ECN weather station from 1999 onwards; these data were gap-filled using the UK Met Office gridded data, which were first locally calibrated via regression against the ECN weather station data.

Evaluation results
Modelled WTDs were compared to observed WTDs for all sites (see main text Fig. 2 and figure   below). Model performance was best at sites where dipwells monitored intact peat, producing the best results in both WTD position and fluctuations. Fluctuations were reasonably well predicted even in drained sites, but WTD position became less well predicted. This is to be expected, as after drain RMSE was 2.9 cm for both of these sites. RMSE at Lake Vyrnwy was only 3.9 cm, and minimum WTD was predicted to within 2 cm. However, for Oughtershaw Moss blocked and drained sites, minimum WTD was only predicted to within 3 -5 cm and RMSE was 4.5 cm and 6.6 cm respectively. Maximum WTD was predicted to within 1.8 cm at Moor House, and 0.7 cm at the Oughtershaw Moss drained site, but within only 2 -5 cm elsewhere. Therefore, values reported here are of a similar magnitude to those reported for other models.