A coupled modelling framework to assess the hydroecological impact of climate change

Rivers are among the ecosystems most sensitive to climate change. Whilst methods quantifying the impact and uncertainty of climate change on flow regime are well-established, the impact on hydroecological response is not well understood. Typically, investigative methods are qualitative in nature or follow quantitative methods of limited scope, whilst the effect of uncertainty is frequently minimised. This paper proposes a coupled hydrological and hydroecological modelling framework to assess the impact of climate change on hydroecological response quantitatively. The characterisation and reduction of modelling uncertainties was critical to the development of the framework. The ability of the framework is illustrated through application to a case study river, the River Nar, Norfolk, England, using the UKCP09 probabilistic climate projections (high emissions scenario, SRES A1F1). The results show that, by the 2050s, a reduction in instream biodiversity is virtually certain if future emissions follow the assumptions of SRES A1F1. Disruption to the natural low flow processes, essential to ecosystem functioning, is also indicated. These findings highlight the importance of the framework in water resources adaptation, particularly with respect to future environmental flows management.


Introduction
The global climate system is changing, with changes to climatic behaviour (mean and variability) projected beyond the 21 st century (IPCC, 2014). Climate change is expected to amplify existing pressures on natural resources, as well as create new ones (IPCC, 2012). Amongst these, freshwater is considered the most essential (Vörösmarty et al., 2010); rivers and their ecosystems provide a diverse range of services upon which humans are dependent (Yeakley et al., 2016), these include: fresh water supply for human consumption; hydro-hazard regulation; and water purification (Gilvear et al., 2017). It is thus through freshwater resources, particularly rivers, that some of the most significant impacts of climate change will be felt (Ostfeld et al., 2012). Consequently, there are significant questions over the long-term sustainability of water resources (Gleick, 1998(Gleick, , 2016Klaar et al., 2014). It is clear that effective water management is central to successful climate change adaption (Ostfeld et al., 2012).
Climate is a major determinant of hydrological processes, where precipitation, temperature and evaporation represent the dominant drivers (IPCC, 2007). Consequently, a changing climate will inevitably lead to alterations of river flow regimes (Rahel and Olden, 2008;Arnell and Gosling, 2016). Attempts to model the impact of climate change on water resources have been ongoing since the mid-1980s (Arnell and Reynard, 1996;Christierson et al., 2012).
Climate projections are, however, subject to large unquantified uncertainties (Murphy et al., 2004), leading to concerns over their suitability for water resources adaptation and planning (Kundzewicz et al., 2008;Wilby, 2016). Examples of these uncertainties include (Clark et al., 2016;Wilby, 2016): (1) epistemic uncertainty, the inability to properly capture the underlying processes and feedbacks; and (2) accounting for variation due to natural climatic variability. In practice, uncertainty dictates the usefulness of the model. Inaccurate appreciation of this uncertainty precludes meaningful interpretation of the model, leading to sub-optimal decision-making (Warmink et al., 2010) when considering future projections. Clark et al. (2016) posit that, research which focuses on characterising, reducing and representing (quantifying) these uncertainties may allow for the provision of plausible flow projections under climate change. Difficulties with regards to the quantification of climate uncertainty may be addressed through the use of a perturbed physics ensemble: an ensemble of GCMs where variation of the model parameters allows quantification of uncertainty (Murphy et al., 2004;Clark et al., 2016). Such  T determinant of ecological health in riverine ecosystems (e.g. Power et al., 1995;Lytle and Poff, 2004;Arthington et al., 2006). Alteration of this natural flow regime threatens the ability to provide ecosystem services (Rahel and Olden, 2008;Vörösmarty et al., 2010;Arthington, 2012). Despite this, and perhaps surprisingly, the effects of climate change on river health are rarely considered, as observed in Durance and Ormerod (2007). Seven years later, Schlabing et al. (2014) describes little change, noting that even when accounted for, the methodology employed is often over simplified and rudimentary. A brief review of such work performed over the last two decades is thus indicated (see also Fig. 1). A large number of the studies investigating the impact of climate change on hydroecological response have been qualitative in nature ( Fig. 1, level 1 to 3 directly; not pictured); examples include Meyer et al. (1999); Ostfeld et al. (2012); Filipe et al. (2013) and Death et al. (2015). Whilst the number of quantitative studies has increased, their scope is often limited to the direct links between climate (temperature) and hydroecological response ( Fig. 1; for example, Durance and Ormerod, 2007;Kupisch et al., 2012 andJyväsjärvi et al., 2015). In these studies, the impact of the altered flow regime is not considered and rarely acknowledged. Döll and Zhang (2010) were the first to consider the impact of climate change on the flow regime at a global-scale. Assessment of the hydroecological impacts was qualitative in nature, with limited consideration of changes in the number of endemic fish species (counts are not considered meaningful bioindicators; for example see Li et al., 2010). The authors acknowledge that quantitative estimates of ecosystem response "have not yet been derived". Following this, studies of a similar nature have been undertaken at higher resolutions (catchment level); examples include Tang et al. (2015); Hassanzadeh et al. (2017) and O'Keeffe et al. (2018). Advances have also been made in the assessment of direct climate change impacts on the provision of freshwater ecosystem services (see conceptual framework in Pham et al., 2019); though again these are, at present, qualitative in nature. Merriam et al. (2017) perform a habitat assessment using a coupled stream temperature and hydrological model. Whilst the focus is on the availability of thermally suitable habitats for brook trout, and not flow alteration directly, the study represents an important advancement towards fully quantitative assessment of the instream ecological impacts of climate change. Nevertheless, significant questions arise as to the robustness of the applied methodology. The authors make assertions, using phraseology such as "high degree of certainty", when discussing results based upon R-squared values and RMSE, a statistical measure of average inaccuracy. Yet, the problems inherent to RMSE have been recognised for some three decades (Willmott et al., 1985), and more recently, Willmott and Matsuura (2005) conclude that, in the context of climate study, model-performance evaluations based primarily on RMSE are questionable and should be reconsidered. Further, issues arise in the calibration of the hydrologic model, where performance is assessed in terms of Nash-Sutcliffe, a statistic subject to long-standing, broad, and sustained criticism (Legates and McCabe, 1999;Seibert, 2001;Criss and Winston, 2008). Indeed, Clark et al. (2016) state that, when modelling the hydrological impacts of climate change, Nash-Sutcliffe (and similar efficiency criteria) introduces additional, unaccounted uncertainties. This disregard of uncertainty throughout the paper calls into question the validity of the results. It is clear that methods to quantify the impact, and associated uncertainties, of climate change on the flow regime are well-established ( Fig. 1, level 1 and 2). The hydroecological implications are less well understood and are rarely considered quantitatively; where attempts have been made, the effect of uncertainty is underplayed. Consequently, the impact of climate change on hydroecological response is unclear, and the fallout for ecosystem services poorly understood. This paper proposes a coupled hydrological and hydroecological modelling framework to assess the impact of climate change on hydroecological response quantitatively. The development of (each stage of) the framework has centred around the characterisation and reduction of uncertainty, in line with the recommendations in Clark et al. (2016). The outputs from this framework are quantitative hydroecological projections of climate change impacts. These outputs are intended to support water resources adaptation, for example in the equitable allocation of water for human use and the environment (known as environmental flows). In order to validate and demonstrate the ability of the framework, this paper features an application to a case study river, the River Nar in Norfolk, England.

Framework
An overview of the three main stages of the proposed framework is presented in Fig. 2. In stage 1, the hydroecological model is developed based on advances made in: (1) Visser et al. (2017), where lag in ecological response, an important component of flow variability (Monk et al., 2017), is accounted for through the consideration of multi-annual hydrological indicators; and (2) Visser et al. (2018a) present an information theory (IT) approach to minimise and quantify structural and parameter uncertainty. The second stage of the framework is the parameterisation of the hydroecological following a modified covariance approach (Visser et al., 2018b). The modified covariance approach focuses on the replication of specific hydrological characteristics (identified in stage 1), whilst also addressing a number of known limitations and uncertainties in hydrological modelling. In stage 3, climate projections serve as the input to the coupled model, providing the quantitative hydroecological projections of climate change impacts. Application of the framework to a case study river catchment is subsequently considered in 3. Case study application.
A holistic depiction of uncertainty was central to the development of the proposed framework. Additional commentary on the characterisation and reduction of the sources of uncertainty, following Clark et al. (2016), is provided in Appendix A.
In the development of this framework it is necessarily assumed that the hydroecological relationship remains stationary (as in A.G. Visser et al. Environmental Modelling and Software 114 (2019) 12-28 hydroclimatological modelling). The evolution of such relationships remains an unknown and is beyond the scope of this paper. A further important limiting factor of many hydroclimatological studies is the focus on extreme climatic events (e.g. Filipe et al., 2013;Thornton et al., 2014;Death et al., 2015). Climate models are well-known for their ineffective simulation of extreme climate, particularly with regards to precipitation (IPCC, 2014); knowledge of the impacts of extreme events is therefore limited. In an effort to address this, the UKCP09 climate projections distribution tails are clipped (5% and 95% probability levels; Murphy and Sexton, 2013). In addition, changes in climate mean and variability may lead to compound events or clustered multiple events; these events are not extreme in themselves but can lead to extreme events and/or impacts (IPCC, 2012). Essentially, severe impacts can occur from minor climatic events. The focus of the framework is therefore on these impacts rather than stochastic (individual extreme) events.

Stage 1 -hydroecological model
Statistical methods are well-established for the testing of hydroecological hypotheses, these include: multiple linear regression (for example, Clarke andDunbar, 2005 andMonk et al., 2007), and multilevel models (recent examples include Bradley et al., 2017 andChadd et al., 2017). Hydrological indicators (HIs) and ecological data serve as the basis for the development of these models. With their sensitivity to change, macroinvertebrates are ideal indicators of river health (Acreman et al., 2008;EA, 2013). This response is determined by considering macroinvertebrate flow velocity preferences as described by the Lotic-Invertebrate Index for Flow Evaluation (LIFE; Extence et al., 1999), a weighted index which takes into account the flow velocity preferences of the macroinvertebrate community; LIFE scores can range from one to twelve, indicating a preference for limited flow & standing water to rapid flows respectively.

Data
The structure of the benthic macroinvertebrate community is subject to change throughout the year. Typically, peaks of activity occur in spring (AMJ) and autumn (OND). Seasonal focus in hydroecological modelling is determined by factors such as the quantity and quality of the available data and the overall modelling objective. For example, fishing is vital to the communities along the case study river, the River Nar (Garbe et al., 2016). If the goal was to preserve future brown trout populations, then modelling efforts would seek to protect their primary food source, Ephemeroptera baetidae (mayfly), which hatch during the spring season. Macroinvertebrate data may be utilised at the species or family level; however, it should be noted that the use of family level data may mask species-specific information (Monk et al., 2012), leading to a reduction in accuracy (Extence et al., 1999).
A hydroecological dataset is created by pairing the ecological data with HIs. These indicators should be ecologically relevant, reflecting the five facets of the flow regime required to support the riverine ecosystem (Richter et al., 1996): magnitude, frequency, duration, timing and rate of change. To date, over 200 ecologically relevant hydrologic indices have been proposed (Olden and Poff, 2003;Monk et al., 2006;Thompson et al., 2013). Should seasonality be present in the hydrological regime, the time-series is split into relevant hydrological seasons; the HIs are then calculated for each. A number of studies (on groundwater-fed rivers) have observed a delay in macroinvertebrate response (for example, Boulton, 2003;Durance and Ormerod, 2007). Visser et al. (2017) and Visser et al. (2018a) propose the incorporation of time-offset HIs to account for this effect. The timeoffset may require fine-tuning if the number of indicators cannot be sufficiently reduced in the steps below; beyond this, no additional work is required.
With a large number of HIs, both variable redundancy and computational effort represent significant challenges. In response, Principal Component Analysis (PCA) is applied, allowing only those indices A.G. Visser et al. Environmental Modelling and Software 114 (2019) 12-28 which describe major aspects of the flow regime to be identified; following Olden and Poff (2003), the most relevant indices are selected proportionally from the five facets of the flow regime described above. The ecological data is then be paired with this set of ecologically relevant HIs.

Statistical modelling
As further aspects of hydroecological relationships are understood, such as ecological lag in response, the likelihood of modelling errors and uncertainty is increased. To account for this, the proposed framework makes use of an IT approach to determine the structure of the hydroecological model (after Visser et al., 2018a). The IT approach provides a robust measure of both structural and parameter uncertainty (see Appendix A.1) as well as a measure of the statistical importance of the model parameters (HIs; a central factor in the parameterisation of the hydrological model, stage 2).
The application of the IT approach consists of 4 steps; for a more detailed discussion see Appendix B.1 or Visser et al. (2018a). To summarise: (1) candidate models are evaluated with respect to the secondorder bias corrected Akaike Information Criterion (AICc, equation (B1); after Burnham and Anderson, 2002); (2) a best approximating model is inferred from a weighted combination of all the candidate models; (3) the parameters are ranked, such that the highest value (Akaike weight, equation (B3)) represents the most important in the model; (4) measures of uncertainty (structural and parameter) are made.
In the development and application of the framework, the IT approach was applied using the R package glmulti (Calcagno, 2013), developed and applied in a relevant discipline (see Isbell et al., 2011). In glmulti, a genetic algorithm (GA), a type of optimisation that mimics biological evolution, is used to select a subset of models (each assessed based on the above IT approach). The GA incorporates an immigration operator, allowing removed HIs to be reconsidered. Immigration sees the level of randomisation increase, and hence the likelihood of model convergence on the global optima rather than some local optima (Calcagno and de Mazancourt, 2010). Inference from a consensus of 5 replicate GA runs has been shown to greatly improve convergence (Calcagno and de Mazancourt, 2010). The multi-model average is subsequently derived from this subset of models. Parameters where the estimate and confidence intervals are zero (i.e. certainty that the index is not to be included), are then removed. In line with Anderson (2007), the set of model parameters is reduced to those accounting for 95% of the cumulative information (see Appendix B.2).
For validation of the hydroecological model see Appendix C.2.

Stage 2 -hydrological model
The HIs identified in stage 1 represent those characteristics of the flow regime which dominate ecological response. Driven by climate projections (Fig. 2), changes to these HIs may be determined from flow time-series simulated via hydrological model. Climate projections, input to the coupled hydrological and hydroecological model, allow the impacts of climate change on hydroecological response to be determined quantitatively (stage 3). This second stage of the proposed framework focusses on the parameterisation of the hydrological model. Clark et al. (2016) highlight model parameterisation as a major source of uncertainty. Typically, hydrological models are parameterised following a split-sample calibration-validation approach, with calibration focussing on the goodness-of-fit between observed and simulated flow. Limitations of the approach are widely acknowledged, these include (Westerberg et al., 2011;Clark et al., 2016): (1) bias in the model parametrisation as the result of disinformative data (Pelletier, 1988;Montanari et al., 2013); (2) the arbitrary nature of GOF statistics; and (3) equifinality (Beven, 2006).
In this proposed framework, the modified covariance approach (Visser et al., 2018b), based on Vogel and Sankarasubramanian (2003), is applied in an attempt to address these limitations (see also Appendix A.1). In Visser et al., 2018b, comparison relative to studies with similar modelling objectives (the simulation of ecologically relevant HIs) showed improvement in both model performance and consistency. A further major advantage of the approach lies in the focus on identifying the region of parameter space which best captures the characteristics of the HIs, providing a greater understanding of model suitability, limitations and uncertainty.

Hydrological model
To further minimise uncertainty, a parsimonious lumped hydrological model should be selected. In the development of the framework, the daily models from the GR (Génie Rural) suite of hydrological models were considered (GR4J, GR5J and GR6J; 4-6 free parameters). The GRJ models have been applied in a variety of hydrological con-  Coron et al. (2017). With observed moments lying outwith the simulated moments (see section 2.2.2), the five and six-parameter models GR5J and GR6J were rejected.
Continuous (daily) time-series of flow, precipitation and potential evapotranspiration serve as model input. The time-series should be of sufficient length for validation on the climate baseline in stage 3 (for example, the UKCP09 baseline is 1961-1990).

Modified covariance approach
The hydrological model is parameterised following the modified covariance approach, as set out in Visser et al. (2018b). In using this approach, the modelling objective is not the replication of a flow timeseries, rather, it is the identification of the region of parameter space which is best able to replicate the HIs. For a more extensive discussion of the modified covariance approach see Visser et al. (2018b).
In the application of this approach, the complete parameter space of the hydrological model is sampled. The number of parameter sets is dependent upon the number of free parameters and the level of uncertainty adjudged acceptable. To reduce bias, the parameter space should be sampled uniformly; for example, using Sobol quasi-random sequences (a Quasi-Monte Carlo method; Caflisch, 1998). The parameter sets thus established, the hydrological model is run in simulation mode using observed climate data. For each of the n time-series, the covariance (between observed climate and simulated flow) is calculated; this is repeated for the observed flow data. The HIs, identified in stage 1, are then determined from both the observed and simulated flows.
Prior to the selection of a parameter set, it is first necessary to validate the hydrological model structure. This is facilitated through plots of the observed and simulated relationships between the covariances and HIs. The model is validated if the moments agree, i.e. observed moments lie within the simulated moments (sampled parameter space). Error thresholds, in combination with index importance (determined in stage 1), are used to identify a suitable parameter set. A linear relationship between the minimum and maximum error thresholds and index importance is defined. Parameter sets which fall below this defined limit are rejected. For additional details see section 3, Case Study Application or Visser et al. (2018b).
The focus of the covariance approach is on the replication of specific hydrological characteristics in the catchment (the HIs), as opposed to flow. Consequently, the hydrological model should be assessed in terms of its ability to replicate these characteristics rather than the observed flow time-series. Indeed, the replication of the time-series is anticipated to be poor, consistent with similar work focussed on the replication of catchment characteristics (e.g. see Seibert, 2000).

Stage 3 -projections
The UKCP09 Weather Generator (WG) was selected due to its ability A.G. Visser et al. Environmental Modelling and Software 114 (2019) 12-28 to represent natural (climatic) variability (Murphy et al., 2009;Kay and Jones, 2012). This consideration of natural variability allows extraordinary (low probability) climatic events to be captured more effectively (Schlabing et al., 2014), which is particularly important for ecosystems (Wigley, 1985). The WG creates synthetic stochastic timeseries of climate variables based on observed climate statistics. The WG is perturbed to represent future climate through the application of change factors. Projections are at a 5 km resolution, allowing for representative simulation across smaller catchments (< 1000 km 2 , typical in the UK; Kilsby et al., 2007;Jones et al., 2010). Data requests are submitted using the UKCP09 web-based portal (http:// ukclimateprojections-ui.metoffice.gov.uk/ui/admin/login.php). The climate variables of interest are precipitation and potential evapotranspiration; note that, potential evapotranspiration may also be computed from the hourly time-series. The CMIP4 SRES scenarios upon which UKCP09 was based does not assign probabilities to specific emissions scenarios (Wigley and Raper, 2001;Meehl et al., 2007;Murphy et al., 2009); consequently, it is assumed that each emissions scenario is equally probable (Murphy et al., 2009).
For the validation of the WG output, UKCP09 recommend comparison of the observed and baseline climate data in the form of bimonthly and seasonal plots of the mean and 95% confidence intervals (for each climate variable; DEFRA, 2011). To this end, linear bias correction is applied bi-monthly (where necessary).

Baseline validation
The baseline climate data is used to validate the framework. The generated climate variables are input into the hydrological model, generating a range of possible flow time-series; for each time-series, the important HIs are calculated (per hydrological year/season). These indices can be assessed relative to the observed indices (determined in Stage 2). Validation is through cumulative distribution and probability density functions (CDF and PDF respectively), comparison of the mean and 95% confidence intervals.
With these indices and the hydroecological model, a range of possible LIFE scores for the baseline period may be determined; validation is as above. If the length of the ecological time-series is insufficient, an alternative approach may be applied; this is further considered in the application of the framework.

Future hydroecological projections
Simulation of future hydroecological projections (LIFE) is analogous to the validation, with the exception that the future climate projections serve as the input data. Each emissions scenario/time-period should be considered distinct.

Case study application
The ability of the framework is both validated and demonstrated through application to a UK case study river. Descriptions focus on the case study specific data acquisition and preparation, the subsequent analysis being as per the described framework.

Study area
The Nar represents a vulnerable and important river type (groundwater-fed chalk stream), already subject to significant stress (NRT, 2012). The additional threat of climate change to its ecological potential cannot be understated. It is intended that the power of this new proposed framework be illustrated in its application to this case study river.
The spring-fed River Nar rises in the Norfolk chalk hills (52.749°N, 0.812°E), 60 m above sea level, flowing west for 42 km before joining the River Great Ouse (52.748°N, 0.394°E). The formation of the fen basin, and resultant dissection of the chalk, created two distinct river units, delineated by a significant gradient change at Narborough (Fig. 3; Sear et al. (2005)). With a greater abundance and quality of ecological data (Visser, 2015), the focus is on the 153.3 km 2 chalk sub-catchment.
The hydrology of the River Nar is characteristic of pure chalk streams (Sear et al., 1999), with a high Base Flow Index (0.91; Sear et al., 2005) and relatively low flows: mean 1.12 m 3 /s, Q10 2.03 m 3 /s and Q90 0.49 m 3 /s over the available record, where Q10 and Q90 represent the 10% and 90% flow exceedance respectively (equivalent to 90 th and 10th percentiles). A reliance on groundwater results in a highly seasonal flow regime, where aquifer recharge primarily occurs in winter months, leading to a progressive rise in river flow until March/ April.
Hydrological data was extracted from the National River Flow Archive (1990-2012CEH, 2018) at the Marham gauge (52.678°N, 0.548°E; Fig. 3). The hydrological data was subdivided into six subsets: two hydrological seasons, winter (ONDJFM) and summer (AMJJAS) and three time-offsets (0-2 years). A total of 63 6 × ecologically relevant HIs were considered; this was reduced to 29 through PCA. This reduced set of HIs were then paired with the LIFE scores to create the hydroecological dataset.

Stage 2 -hydrological model
For parameterisation of the hydrological model, 54 years of daily average flow recorded at the Marham gauge were extracted (September 1961(September to 2015. The corresponding climate variables were computed from daily average rainfall and hourly temperature data at 5 MIDAS stations in and around the catchment ( Fig. 3; Met Office, 2016). The parameters of interest are the average daily precipitation (P) and potential evapotranspiration (PE); P is determined via the computation of the daily catchment average rainfall, whilst PE is estimated from hourly temperature data using the temperature-based PE model from Oudin et al. (2005).
In order to verify the method of investigation, n = 100,000 parameter sets were generated using Sobol sequencing. The HIs used to parameterise the model are those indicated by the hydroecological model in stage 1.
In the parameterisation of the hydrological model, the minimum and maximum error were specified as 17.5% and 35% ( error 2 min × ) respectively, from which the linear threshold was determined (the relationship between the minimum and maximum allowable error and the relative importance of the variables; covariances were assigned a relative importance of 1).

Stage 3 -projections
To address a number of the uncertainties indicated in the introduction (see also Appendix A.2), the UKCP09 probabilistic climate projections are used. The 2050s (2040-2069) high emissions scenario (A1F1 SRES) is considered. This emissions scenario is approximately equivalent to a change in temperature of 4.3°C by 2081-2100 (relative to the pre-industrial period 1850-1900; Riahi et al., 2011;Met Office, 2018b). For demonstrative purposes, the UKCP09 WG was run for the full range of 10,000 variants.

Stage 1 -hydroecological model
The hydroecological model, a linear multi-model average, is depicted in Equation (1) (overleaf). Summaries of the HIs are provided in Table 1; importance represents the relative weight of evidence in support of each index in the model (according to IT), whilst the relative parameter uncertainty is the 95% confidence interval relative to the parameter estimate. The underlying hydroecological processes are first considered, followed by a review of the predictive ability and uncertainty associated with the hydroecological model. (1)

Underlying hydroecological processes
The winter hydrological season, when the chalk aquifer recharges, features both the most and least important HIs, 10R90Log and riseMn respectively. The indicator 10R90Log, ratio of low flows to high flows (10th to 90 th percentiles), is described as an indicator of responsiveness (Richards, 1990). The log-scale of the index, coupled with its importance, means that there is scope for 10R90Log to dominate the hydroecological model, both positively and negatively. Fig. 4 clearly illustrates that large values of 10R90Log correspond with the highest LIFE scores, and vice versa. It is only when 10R90Log is small (∼0) that the other six indices contribute to LIFE score. Varying high and low flows shows that the highest values of 10R90Log, and hence LIFE, are achieved when high flows are medium-high ( P Q 0 exp( 90(log( ))) 1 < < ). Given the log-space, the scope for a negative impact ( P Q exp( 90(log( ))) 1 > ) is large. Surprisingly then, whilst magnitude of flow is of importance for the recharge of the aquifer, higher winter flows may actually negatively impact the macroinvertebrate community. The negative sign of the HI riseMn indicates a preference for a low mean rise rate in winter flows. However, the low importance of the index (Table 1) sees it consistently contribute less than 2.5% to the LIFE score (Fig. 4).
In terms of hydrological season, the summer months (AMJJAS) dominate the hydroecological model. There is an indication that, in the summer months, consistency in flow (low range/variation) is preferred: (1) a sustained increase in flow (RevPos) sees a large negative impact on LIFE (importance = 0.80); (2) though not as important, logQVar similarly implicates large variation in flow; and (3) minimising the range between low and median flows (Q70Q50, Q80Q50 and Q90Q50) has an increasing effect on LIFE.
Looking to Fig. 5, four out of the five summer HIs are lagged (S-1). Essentially, these indicators are influencing the health of the river two years in advance; should there be a bad summer, with lots of variation, the consequences could be severe. However, the presence of the Q80Q50 indicator in the immediately preceding season gives some scope for improvement. However, it is also worth noting that the negative impacts of a 'bad summer' would only be felt if the value of the index 10R90Log was low, whilst if 10R90Log is largely positive or negative, the preceding summers are of limited importance.
In terms of management, it is clear that summer floods in particular could be detrimental to the health of the river; perhaps representing an argument for improving connections to flood plains. Similarly, extremely high winter flows may be harmful, indicating there may be scope to abstract and store waters during the winter months for use in summer. However, it is worth noting that negative impacts are also a necessary component of the proper ecosystem function; for example, they might act as a 'natural reset button' (Everard, 1996;Lake, 2003). Interestingly, the majority of the indices are dimensionless (with the exception of logQVar and riseMn), this allows for some scope for variation in flow without causing excessive damage; for instance, in summer, a need for increased abstraction need not necessarily be a detriment to river health (though this assumption ignores the other effects of decreased flow).

Predictive ability and uncertainty
The predictive ability of the model is first indicated by the relative parameter uncertainty (unconditional variance, or 95% CI, relative to Fig. 3. River Nar catchment. Data sampling and recording sites/stations are marked; climate data is recorded at the MIDAS stations. Inset: Location of the River Nar in the UK.

Table 1
The HIs indicated in the hydroecological model in descending order of importance. The facets of the flow regime are denoted as M (magnitude) and R (rate of change); the hydrological seasons are indicated by W (winter) and S (summer). the parameter estimate; Table 1). Generally, as relative parameter uncertainty increases, the importance of the index decreases; this is one of the advantages of the weighting of the HIs in the IT approach (Visser et al., 2018a). The fact that the most important index, 10R90Log, has greater uncertainty than the second most important, RevPos, suggests that this may be the best parameterisation possible in the model. With regards to the implications of parameter uncertainty, further inference may be made through the consideration of the 10,000 Monte Carlo simulations (Fig. 6). The plot shows that the hydroecological model performs well (low interquartile range of 0.44 and relative error centred around one; perfect agreement). This level of uncertainty is considered satisfactory. Fig. 7 (a) depicts the observed and simulated relationship between the covariance of precipitation and flow, P Q ( , ), and the HI Q80Q50; Fig. 7 (b) depicts the same relationship for the climate variable potential evapotranspiration, PE Q ( , ). For all seven HIs, the observed moments lie within the simulated moments, validating the use of the hydrological model.

Stage 2 -hydrological model
The capacities of the production (x1) and routing (x3) stores were estimated as 511 and 311 mm respectively; the time elapsed for flow routing is approximately 1.17 days (x2). Inflow from the chalk aquifer is represented by a positive groundwater exchange coefficient (x4) of 2.84 mm per day. The level of agreement for all seven HIs is summarised in Table 2. With a value of 0.8, the largest covariance relative error is for potential evapotranspiration; this is considered acceptable   as precipitation is considered the principal determinant of flows in the East Anglia region (Kay et al., 2013). The HI relative errors are below ± 11%, with the exception of the least important index, riseMn (relative error = 34%). Overall, the level of relative error in the hydrological model is considered satisfactory; the impacts of the error in the index riseMn are likely negligible (based on the findings from the hydroecological model). For standard model validation, see Appendix C.3.

Stage 3 -projections 4.3.1. Baseline validation
The ability of the framework to reproduce the observed data, hydrological and hydroecological, is assessed via CDFs and PDFs. For the CDF plots, the observed function should situate within the boundaries of the baseline projections; ideally, centrally. The PDF plots focus on relative error, where a value of 1.0 indicates perfect agreement; here the objective is on a low interquartile range (IQR). In the interests of concision, validation of the hydrological projections centres on the index Q80Q50, selected both due to its high importance (0.80) and ease of interpretation (ratio of moderate low flows to median flows); summary tables for all seven HIs are available in Appendix C.4.

Hydrological model validation.
The validation plots for the HI Q80Q50 are presented in Fig. 8. Fig. 8a/d are based on the UKCP09 baseline ; both satisfy the objectives outlined above: the CDF of the observed values lies within the projections and the PDF shows a low IQR. Comparatively, the 95% confidence interval (CI) appears large, however, given the probabilistic nature of the projections this is not unexpected. The baseline interval for which ecological data is available (1986)(1987)(1988)(1989)(1990) is summarised in Fig. 8b/e. The IQR is similar to the 30-year baseline, with a minor improvement in the 95% CI; given the limited time-period, the right-skew of the PDF (Fig. 8e) cannot be ascribed significance. On the alternative baseline (2010-2017), the CDF is notably stepped (Fig. 8c); this is reflected in the PDF (Fig. 8f) with a local maximum, and a median not equivalent to one (perfect agreement). Despite this, the IQR is the lowest of the three validation plots.
Overall, for the CDF's, the observations lie centrally within the probabilistic projections, whilst the PDF's reveal low IQR's. The plots satisfactorily validate of the use of the UKCP09 projections through the hydrological model.

Coupled hydrological-hydroecological model validation.
There is no ecological data (species LIFE) available for the period 1961-1990 (baseline validation period). However, sampling and identification of macroinvertebrates to the family level was carried out during the period 1989-1990, allowing for some comparison. To further address this, an alternative baseline was established through consideration of the earliest time-period considered by the UKCPO9 WG (2010-2039; reduced to 2010-2015), run for the medium emissions scenario for 1000 randomly sampled variants. Subsequently, the climate variables are bias corrected relative to the observed data in this period.
Validation plots for the hydroecological response, LIFE, as predicted by the coupled model, is presented in Fig. 9a/c for the baseline (1986)(1987)(1988)(1989)(1990) and Fig. 9b/d for the alternative baseline (2010-2017); recall that for the period 1986-1990 only family LIFE data is available. On this baseline period, the CDF's ( Fig. 9a) are in agreement, with a small IQR for the relative error of approximately 0.1 (Fig. 9c). The somewhat swollen 95% CI may have a threefold explanation: 1) family level application of the LIFE methodology tends to underestimate hydroecological response (Extence et al., 1999;Monk et al., 2012); 2) the limited number of years/data points; and 3) the probabilistic nature of the projections. For the alternative baseline (2010-2017), the CDF (Fig. 9b) is in agreement. The PDF (Fig. 9d) also reveals a lower IQR (relative to Fig. 9c) as well as an improved 95% CI.
Although the temporal range of the validation is limited, both time periods are able to achieve a satisfactory level of performance, thereby validating the use of the UKCP09 projections and the coupled hydrological-hydroecological model. The use of the coupled model is thus considered fit for purpose in application to future projections.

Future projections
The UKCP09 climate projections for the 2050s time slice of the high emissions scenario were inputted to the coupled hydrological and   (1961-1990; left); the baseline period for which there is corresponding ecological data (1986-1990; middle); and the alternative baseline (2010-2017; right). A.G. Visser et al. Environmental Modelling and Software 114 (2019) 12-28 hydroecological model. The focus herein is on this hydroecological response. The projections of hydroecological response are first considered relative to the baseline through the CDFs and PDFs in Fig. 10. Looking to the means first (dashed lines), the projected change is relatively small, however, there is a consistent increase in the range of possible LIFE scores across the distribution. The increase in maximum LIFE scores appears responsible for the majority of this change, though, some may also be attributes to the minimum values, specifically, the tails of the distribution (percentile < 0.375). In Fig. 4 it was shown that the index 10R90Log was the main determinant of higher LIFE scores; it can thus be presumed that from percentile > 0.75, the increase in LIFE scores is the result of an increased stability in the ratio of high to low flows in the winter season. Where percentile < 0.75, the five summer HIs are likely to dominate, a further indication of increased stability of flows in the river. Fig. 11 gives an indication of the probability of these hydroecological projections. These probabilities are in line with calibrated language used by the IPCC since AR4 (Treut et al., 2007;Mastrandrea et al., 2010, Table 3). Widening the confidence interval, from about as likely as not to virtually certain increases likelihood but results in a wider estimate. Overall, the bounds of uncertainty are relatively narrow, however, it is clear that the greatest confidence lies in the interquartile range, rather than the tails of the distribution. It should also be noted that, whilst the change to maximum LIFE scores is still in evidence, the decrease in LIFE scores at the lower distribution has disappeared, indicating a lack of certainty in those projections.
Whilst Schlabing et al. (2014) also observed limited changes in the central tendencies, they also note that it is important to consider the tails of the distribution. Although these events lie outwith the probabilities indicated in Fig. 11, this may be justified due to the ability of the WG to capture low-frequency events (Dubrovský et al., 2004;Mehrotra and Sharma, 2007). As in Schlabing et al. (2014), Fig. 12 looks to the hydroecological response at the 5th and 95th percentiles. It is important to note that LIFE scores at these percentiles will be primarily determined through the winter HI 10R90Log.
At the 5th percentile, LIFE scores < 4.5 account for only 0.016% of observations and are therefore omitted. Broadly speaking, the frequency of lower LIFE scores appears to decrease under the future projections. Consequently, there is almost a 10% increase in the number of LIFE scores equal to 7. For the 95th percentile, LIFE scores > 9.5 account for 0.0244% of the total observations and are therefore omitted. With a marked increase in the frequency of LIFE ≥ 8, the positive change in hydroecological response previously observed (Figs. 10 and 11) is clearly in evidence.

Implications for the case study river
The proposed framework has indicated a clear hydroecological response to the projected changed climate under the A1F1 high emissions scenario in the 2050s. However, the magnitude and direction of change is projected to be both small and positive. The scale of this change is in line with the UKCP09 projections for the East Anglia region under this scenario. In this region, the projected change in mean annual precipitation is small, ranging from ± 5% across the 10th to 90 th percentiles (Met Office, 2018); note that Kay et al. (2013) observe that hydrological response in East Anglia is principally determined by precipitation. Further, the range of LIFE values is known to be small,  A.G. Visser et al. Environmental Modelling and Software 114 (2019) 12-28 particularly in the IQR (based on the BIOSYS database of ecological data across 548 catchments in England (Environment Agency, 2018), the IQR for 546 catchments is 1). Therefore, the observed change signal may be presumed to represent a true change in community structure. Figs. 10 and 11 indicated a clustering of LIFE scores; visually this is most clear with the 1.375 increase in the mode (Fig. 10b). This is reflected in the summary statistics, for example, the kurtosis of the LIFE score distribution increases from 16.6 to 18.4. This flattening of the hydroecological response is a possible indication of a reduction in biodiversity. If this were to be the case, this would increase the vulnerability of the river overall, monocultures being far more susceptible to local extinction. Further, the reduced frequency of events where LIFE scores fall very low, could impact negatively upon the river, robbing the ecosystem of vital, natural reset events (Everard, 1996;Lake, 2003).

Framework limitations
Limitations of the framework centre around the assumptions of stationarity and data availability. In climate change modelling, projections are often based on historic climate, with the assumption that the statistical properties of the climate remain stationary. This assumption is inherited under both hydroclimatic and coupled modelling. The corollary is an enforced assumption that ecological response remains the same as it is now; consequently, at this time, it is not possible to account for the adaptive response of the riverine community.
A barrier to hydroecological studies has been the lack of paired long-term hydrological and ecological time-series (Monk et al., 2007(Monk et al., , 2012. This problem persisted in the development of the hydroecological model. For example, in the UK, routine macroinvertebrate sampling began circa 1990 (Murray-Bligh, 1992). Therefore, given the baseline of 1960-1990, validation is limited. To address this, an alternative baseline was derived. The use of climate projections with a more up-to-date baseline, for example, the soon to be released UKCP18 projections or a WG trained using CMIP5 or CMIP6 climate data would also address this.

Conclusions
The implications for flow regime make rivers among those ecosystems most sensitive to climate change (Death et al., 2015;Watts et al., 2015). Whilst studies have attempted to assess the impact of climate change on hydroecological response, methods are often qualitative or follow quantitative methods of limited scope. The resulting lack of clarity renders the fallout for ecosystem services effectively an unknown. In answer, the proposed framework provides a quantitative approach, developed using methods which minimise uncertainty (in line with recommendations in Clark et al., 2016). Fig. 11. Projections of hydroecological response bounded by the IPCC likelihood confidence intervals set out in Table 3.

Table 3
Uncertainty terminology as used by the IPCC (Treut et al., 2007;Mastrandrea et al., 2010); the third column indicates the confidence interval specified in this study.
A.G. Visser et al. Environmental Modelling and Software 114 (2019) 12-28 The ability of the framework has been illustrated through application to a case study river. The projected hydroecological response in April-June, the period of peak MI activity in the river, is considered under the A1F1 high emissions scenario in the 2050s. The hydroecological response is in line with climate projections for the East Anglia region. The projections indicate that a reduction in biodiversity is virtually certain; a possible disruption to low flow processes essential to ecosystem functioning is less strongly indicated. It should be noted that, whilst the projected hydroecological change may be limited, the River Nar is strongly influenced by groundwater (BFI = 0.91). Consequently, the impact of changes in precipitation may be reduced; thus, greater change in response might be expected in catchments where surface runoff dominates.
In summary, the proposed framework serves as a new and dynamic tool with the potential to provide valuable information in the pursuit of more accurate assessments of the impact of climate change on river ecosystems. Critically, and possibly uniquely in the field (Bennett et al., 2013), the end user will also be provided with a quantifiable measure of uncertainty in the hydroecological projections. Further applications of the framework include water resources planning and future environmental flow management. In recent years, hydroecological modelling is often undertaken using a regime-based spatial framework (for example Monk et al., 2011;Zhang et al., 2012). In a similar manner, the proposed framework may be extended to cover multiple rivers of similar flow regime classification. Such generic projections of the impact of climate change on hydroecological response might thus be used to plan wider adaption measures, including for ungauged rivers, where appropriate. The projections may also be used to assess the implications of climate change on the provision of instream ecosystem services (e.g. through the framework set out in Ncube et al., 2018).

Data availability
Consent has not been given to share the data used in this study; however, these data are freely available from the original sources: the Environment Agency (EA, 2018; macroinvertebrate sampling records); Met Office (2016; climate, precipitation and temperature); National River Flow Archive (CEH, 2018; gauged flow at Marham); and data requests for the climate projections may be submitted using the .
Supplementary data to this article can be found online at https://doi.org/10.1016/j.envsoft.2019.01.004. Table A1 Types and sources of uncertainty that are addressed in stages 1 and 2 (hydroecological and hydrological modelling) of the proposed framework. The modified covariance approach focusses on replicating the essential characteristics of the catchment explicitly. The relative importance of the HIs is known; error thresholds for each HI may be specified accordingly. Equifinality

A.1 Stages 1 and 2
The modified covariance approach considers the full parameter space which is narrowed down to a small region which is able to best replicate the HIs of interest.

A.2 Stage 3 -Climate projections
Climate change projections are central to the application of the proposed framework. It is recommended that probabilistic climate projections, which consider a range of impacts, be used when applying the proposed framework. In the case study application, the UKCP09 probabilistic climate change projections are used, specifically, the weather generator product. The application of the framework is not limited to UKCP09, other sources of probabilistic climate change projections include: the COMEPRO project in the Mediterranean region (Kaspar-Ott et al., 2016;Ring et al., 2018); the MIT IGSM-CAM framework (Monier et al., 2013a), applied over Northern Eurasia (Monier et al., 2013b); and UKCP18, the next iteration of UK Climate Projections, based on the Research Concentration Pathways from AR5 (Met Office, 2018a). Equally, projections may be produced via weather generator may be used directly; for example, the Vector-Autoregressive Weather Generator (Schlabing et al., 2014).
The UKCP09 identifies three major sources of uncertainties in their climate projections (Murphy et al., 2010): epistemic uncertainty (incomplete understanding of climate system processes), internal climate variability, and scenario uncertainty. A summary of the controls introduced in UKCP09 to minimise this uncertainty is detailed in Table A2. A.G. Visser et al. Environmental Modelling and Software 114 (2019) 12-28  Sources of uncertainty in the UKCP09 weather generator climate projections and the controls in place to minimise this; based on Murphy et al. (2010).

Source Controls
Epistemic uncertainty A perturbed physics ensemble of the variance (Clark et al., 2016), e.g. different mathematical representation of the processes, interactions and feedbacks.

Variability
Multiple runs with the same initial conditions for each ensemble. Weather generator simulations based on statistical characteristics in the observed data. Simulations pick up more extreme climatic events (Schlabing et al., 2014). Scenario uncertainty There is a lack of agreement in how relative probability should be assigned to emissions scenario. To address this, UKCP09 presents three emissions scenarios: low, medium and high.

B.1 Model evaluation
Although the overall concept of information theoretics may be unknown to the reader, certain aspects should be familiar. Based on "deep theoretical foundations" (Burnham and Anderson, 2001, p. 244), the concept and application are conspicuously simple. Candidate models are evaluated over three steps: 1) measuring the information lost in each approximating model; 2) determination of the evidence in support of each model; and 3) multi-model inference of a final model structure from the candidate set.
Step 1 Loss of information from model f Kullback-Leibler (K-L) gives a measure of the amount of information that is lost when model g is used to approximate reality, f . A model which loses the least information, i.e. has the most supporting evidence out of the candidates, can be considered the best approximation of reality.
The information loss I f g ( , ) is determined through computation of information criteria (IC). A multitude of IC exist, some of which with the reader is undoubtedly familiar. The Akaike Information Criterion (AIC) represents the standard estimate of Kullback-Leibler information (Burnham and Anderson, 2002). In hydroecological modelling, the sample size, n, is often small relative to the number of variables, K . A second order bias corrected version of AIC, AICc, can account for this through the addition of a second penalty (Burnham and Anderson, 2002): Step 2 Evidence in support of model g i The value of AIC is dependent on the scale of the data, the goal is to achieve the smallest loss of information given the data. This difference is rescaled and ranked relative to AIC min : The value of i may be interpreted through a rule of thumb (based on likelihood intervals): 2 i < , there is substantial supporting evidence for model g i ; 4 7 i , the models are not as competitive; and if 10 i > it can be assumed that there is strong evidence against model g i (Burnham and Anderson, 2002). From this measure of evidence, the likelihood that model g i is the best approximating model can be determined. This is known as the Akaike weight, w, ranging from 1 to 0, for the most and least likely models respectively: The use of the Akaike weight allows for clearer inference when considering the candidate models.
Step 3 Multi-model inference When using information theory model selection, the best approximating model is inferred across the entire candidate set. This is achieved through consideration of a weighted combination of all candidates. Parameter averages,ˆ, are simply the sum of the Akaike weights for each model that contains the predictor,ˆ: As a result, the parameter averages are ranked, such that the highest value represents the most important in the model. This eliminates the problem of multiple equally plausible models with different parameter structures (equifinality).

B.2 Application using glmulti
The package glmulti streamlines the above steps into a single function (Calcagno, 2013). The fundamentals of the algorithm and approach are available in Calcagno and de Mazancourt (2010). A.G. Visser et al. Environmental Modelling and Software 114 (2019) 12-28 A subset of models was determined using the function glmulti. The function was applied five times to ensure convergence to a consensus of model subsets, with the function coef applied to determine the IT multi-model average. The number of indices is reduced by removing those indices where both coefficient and standard error are zero and to within the 95% confidence interval, by ordering by descending importance: (importance i / cumsum (importance)) < 0.95; this is illustrated in Table B1 overleaf. Table B1 The structure of the hydroecological model prior to final removals (detailed above). Hydrological seasons are indicated by W (winter) and S (summer); the facets of the flow regime are denoted as M (magnitude) and R (rate of change). The removed indices occupy the final five rows, reasons indicated in bold.

C.1 Hydroecological model -Index contribution
The contribution of each HI to hydroecological response was determined using the baseline data . Each of the 10,000 WG variants and hydrological year were considered independently. For each of these iterations, the relative contribution of the HI was determined: index value index coefficient index value index coefficient (C1)

C.2 Hydroecological model -Validation with observed data
Data represents a key limiting factor to hydroecological modelling, with long-term (> 15-20 years) macroinvertebrate community data, at the species level, uncommon (Monk et al., 2012); Consequently, the length of the time-series in hydroecological modelling is insufficient for splitsampling (calibration and validation); this is commonplace in hydroecological modelling (Monk et al., 2012;Environment Agency, 2018).
The exploration of the model uncertainty serves as one approach to address the validation. Further validation is considered through comparison of simulated species LIFE scores  to three data sources summarised in Table C1; see Figure C1 for validation. The following should be considered when interpreting Figure C1: • As discussed previously, LIFE score differences across taxonomic level are inevitable; • Differences in LIFE scores of the same taxonomic level are due to known errors within the BIOSYS records; BIOSYS stat that corrections to address these inconsistencies are in progress (Environment Agency, 2018); • April 1995-September 1997 saw extremely low rainfall, leading to errors in the recording of low flows (NRFA, 2014). This discrepancy may be the reason for the differences in observed and simulated values. It should be noted that the hydroecological model was not trained on data marked as suspect; this data is included in Figure C1 to allow for the fitting of trendlines.
The two dashed lines represent the comparison of the training data (light green) and model simulations (dark green); the similarity in the slope of the lines indicates a high level of agreement in LIFE scores. Only two additional species LIFE scores are available (2013 and 2014; solid blue circles); it can be seen that these values are consistent with the observed training data trendline (light green dashed line).
Two sources of observed family LIFE scores are available; see above for notes on differences between the data sources. The slope of the trendline for family LIFE scores, determined as part of this study (solid blue line), is similar to the observed training data; the underestimation of LIFE scores may be attributed to the difference in taxonomic level. For the EA BIOSYS data, a good level of agreement is again indicated; though it can be seen that the validity of the model improves over time. Fig. C2. Validation of the hydrological model using observed data. (a) Comparison of the mean and 95% confidence interval (2 standard deviations). (b) Probability density function of the relative error; a value of one indicates perfect agreement between model∼observations.

C.3 Hydrological model -Validation with observed data
A.G. Visser et al. Environmental Modelling and Software 114 (2019) 12-28 C.4 Hydrological model -Validation with projections Differences in the 95% CI for the index 10R90Log are the result of a single outlier observation during the 30-year time period (note that the other two time periods are 4 and 7 years in length). It is also observed that the impact on hydroecological response reduces as 10R90Log increases to more extreme values ( ± 10 specifically).