Modelling indicators of water security, water pollution and aquatic biodiversity in Europe

The GWAVA (Global Water AVailability Assessment) model for indicating human water security has been extended with a newly developed module for calculating pollutant concentrations. This module is first described and then illustrated by being used to model nitrogen, phosphorus and organic matter concentrations. The module uses solely input variables that are likely to be available for future scenarios, making it possible to apply the module to such scenarios. The module first calculates pollutant loading from land to rivers, lakes and wetlands by considering drivers such as agriculture, industry and sewage treatment. Calculated loadings are subsequently converted to concentrations by considering aquatic processes, such as dilution, downstream transport, evaporation, human water abstraction and biophysical loss processes. Aquatic biodiversity is indicated to be at risk if modelled pollutant concentrations exceed certain water quality standards. This is indicated to be the case in about 35% of the European area, especially where lakes and wetlands are abundant. Human water security is indicated to be at risk where human water demands cannot be fulfilled during drought events. This is found to be the case in about 10% of the European area, especially in Mediterranean, arid and densely-populated areas. Modelled spatial variation in concentrations matches well with existing knowledge, and the temporal variability of concentrations is modelled reasonably well in some river basins. Therefore, we conclude that the updated GWAVA model can be used for indicating changes in human water security and aquatic biodiversity across Europe. Editor Z.W. Kundzewicz Citation Dumont, E., Williams, R., Keller, V., Voss, A., and Tattari, S., 2012. Modelling indicators of water security, water pollution and aquatic biodiversity in Europe. Hydrological Sciences Journal, 57 (7), 1378–1403.

Le modèle GWAVA (Global Water AVailability Assessment), caractérisant la sécurité de l'eau à usage humain, a été complété par un module développé récemment pour calculer les concentrations en polluants. Dans cet article, ce module est d'abord décrit puis appliqué à la modélisation des concentrations en azote, phosphore et matière organique. Ce module utilise uniquement des variables d'entrée susceptibles d'être disponibles dans les scénarios futurs, afin que le module soit applicable à de tels scénarios. En premier lieu, le module calcule la charge polluante provenant de la surface du sol vers les rivières, les lacs et les zones humides en considérant des facteurs tels que l'agriculture, l'industrie et le traitement des eaux usées. Les charges polluantes calculées sont par la suite converties en concentrations en prenant en compte un certain nombre de processus aquatiques tels que la dilution, le transport vers l'aval, l'évaporation, le captage d'eau, et les pertes dues aux processus biophysiques. Un risque est estimé sur la biodiversité aquatique si la concentration en polluants modélisée excède certaines normes de qualité de l'eau. Selon les résultats obtenus, ce serait le cas pour environ 35% de la zone européenne, en particulier là où les lacs et zones humides sont abondants. Un risque est estimé sur la sécurité de l'eau à usage humain dans les zones où les demandes en eau ne peuvent être satisfaites en période de sécheresse. Ce serait le cas pour environ 10% de la zone européenne, en particulier dans les zones méditerranéennes, les zones arides et les zones ayant une forte densité de population. La variabilité spatiale des concentrations estimées est en accord avec les connaissances actuelles, et la variabilité temporelle est modélisée de manière satisfaisante dans certains bassins. Ainsi nous concluons que cette mise à jour de GWAVA peut être utilisée afin d'estimer les changements au niveau de la sécurité de l'eau à usage humain et de la biodiversité aquatique à travers l'Europe.

Keywords: water resources; eutrophication; modelling; BOD; nitrogen; phosphorus; Europe; GWAVA; water scarcity; water quality; EU Water Framework Directive; ressources en eau; eutrophisation; modélisation; DBO; azote, phosphore, Europe, GWAVA, rareté de l'eau, qualité de l'eau, directive cadre européenne sur l'eau

The requirement to manage water resources in an integrated and sustainable manner has become a driving force behind the use of large-scale gridded models (Xu and Singh [89]). Such models have traditionally focused solely on availability of water resources for direct human use. However, the quality of freshwater resources for supporting the diversity of native species is an important aspect (MEA [59]), especially because native species can sustain important ecosystem services such as food and recreation (Loomis et al.[54]).

Modelling impacts of changing water resources on biodiversity requires consideration of water quality and its variability over large spatial and temporal scales. This large-scale variability is important because global drivers such as population growth and climate change are likely to have widespread and long-term impacts on water quality. Population growth has widespread impacts on water quality due its effects on agricultural and industrial production and waste disposal (Billen et al.[10], Hoekstra [44]). Global climate change impacts may decrease water quality through reduced dilution capacity of some rivers because of more frequent droughts, or increased pollutant loadings to other rivers due to changed rainfall patterns (Bates et al.[8]).

Five-day biochemical oxygen demand (BOD5), total nitrogen (TN) concentration and total phosphorus (TP) concentration are water quality parameters known to threaten biodiversity if they exceed certain levels (WFD UK TAG [83], [84]). The BOD5 affects oxygen availability in rivers. Native macro-invertebrate communities are often sensitive to changes in oxygen availability, causing them to disappear as oxygen availability deteriorates (Carvalho et al.[15]). BOD5 levels are increased by loading of organic matter, such as discharges from sewage treatment works and storm overflows, and agricultural loadings from slurry and silage liquor. Elevated TN and TP concentrations can lead to increased phytoplankton growth, which can cause eutrophication and reduced biodiversity. Diatoms are sensitive to altered TN and TP concentrations and, during eutrophication, they are often quickly replaced by other (often undesirable) blue-green algae (Hering et al.[43]). Elevated TP concentrations in surface waters mainly result from sewage discharges into rivers, whereas elevated TN concentrations mainly result from agricultural practices, such as manure and fertilizer application and cultivation of N2 fixing crops (Seitzinger et al.[71]).

Until recently, pollutant concentrations in surface waters were not modelled on a continental to global scale. At such scales, only pollutant loading models existed that often treated entire river basins (e.g. Dumont et al.[18]), or sometimes large sub-basins as the basic unit (e.g. Grizzetti et al.[38]). Such models are useful for assessing the impacts of accumulation of pollutants in large water bodies at the outlets of large (sub) basins, such as coastal shelves, lagoons or major lakes (e.g. Garnier et al.[35]). However, we expect that concentrations are more useful than loadings for testing of compliance with water quality standards such as those developed for the European Water Framework Directive (WFD). Many distributed (high-resolution) surface water quality models that are suitable for modelling concentrations on the catchment to country scale have been reported and reviewed in the scientific literature, for example by the EUROHARP project (Andersen et al.[4]). Such models usually require too detailed data (for input and calibration) to make application on a continental scale feasible. Nevertheless, three studies have been published recently of distributed models of pollutant concentrations on the continental scale:

• a. The HYPE model (Lindström et al.[53]) has been used to model daily N and P concentrations across Europe at a 120-km2 resolution. Donnelly et al. ([17]) validated these results by evaluating the bias of modelled long-term average concentrations (referred to as "relative volume errors") at the outlets of large European river basins.
• b. The model of He et al. ([42]), like the HYPE model, explicitly describes the terrestrial N cycle allowing the calculation of daily nitrate concentrations. Their modelled nitrate concentrations cover the whole globe on a 0.5-degree resolution, but validation only took place using multi-year average dissolved inorganic N loads (t N year-1) and yields (t N km-2 year-1) at major river basin outlets.
• c. Vörösmarty et al. ([82]) modelled long-term average BOD, N and P concentrations on a 0.5-degree resolution in combination with 20 other geospatial drivers in order to carry out global analysis of current threats to human water security and biodiversity. In their analyses they thematically grouped BOD, N and P under "pollution" together with soil salinization, mercury deposition, pesticide loading, sediment loading, potential acidification and thermal alteration. They did not validate their modelled long-term average pollutant concentrations, but instead compared their integrated analyses of threat to human water security and river biodiversity to other threat studies.

In this paper, we describe, illustrate and validate a distributed continental-scale model of TN, TP and BOD concentrations that should overcome some of the previously-mentioned shortcomings of existing approaches by modelling on a high spatial and temporal resolution and by validating the spatial and temporal variability of modelled concentrations (not loads). The model described in this paper is based on the GWAVA (Global Water AVailability Assessment, Meigh et al.[57]) model. The GWAVA model is a gridded model for prediction of water resources scarcity at continental and global scales; it has been further developed to include a water quality module. This module enables GWAVA to model concentrations of TN, TP, BOD5 and other pollutants in a spatially-distributed manner at scales ranging from the river basin scale to the global scale. In addition, the new module can be a basis for improving existing continental- and global-scale gridded models of water resources that are currently largely based on water quantity (e.g. Arnell [6], Cherckauer et al.[16], Takata et al.[74], Hanasaki et al.[41]), or that do not include seasonal variability (Vörösmarty et al.[82]).

We validate and illustrate the current risk to human water security and biodiversity predicted by the updated GWAVA model. We developed this model such that it can be run solely using input data available for future projections (e.g. Kämäri et al.[48]). We also developed the model such that it is feasible to parameterize the model for pollutants whose behaviour in the environment is less well known than that of TN, TP and BOD5.

We first describe the pre-existing GWAVA model; then how the updated GWAVA models the levels and pathways of TN, TP and BOD5 from their sources (households, paved surfaces, industry, agriculture) to rivers, lakes and wetlands. We compare the resulting monthly gridded maps of 5-arc-minute (5′) resolution of TN, TP and BOD5 levels with European water quality standards to indicate current aquatic biodiversity risk across Europe. Finally, we combine this with modelled human water security to provide an integrated view.

Water flows and pollutant fluxes were modelled with a monthly time step and on a 5′ grid resolution. Based on this, pollutant levels and indicators of human water security and aquatic biodiversity were mapped across Europe. The methods used are described below.

The GWAVA model was developed by Meigh et al. ([57]) and later improved and extended in different regional and global research projects (e.g. Tate et al.[76], [75], Tate and Meigh [77], Meigh et al.[56], Folwell and Farquharson [33], Farquharson et al.[30], IVL et al.[47]). GWAVA estimates water scarcity on a cell-by-cell basis by comparing modelled river flows with modelled human demand for water (Fig. 1). First, snowmelt, icemelt and rainfall interception by vegetation are modelled according to Bell and Moore ([9]), Rees and Collins ([68]) and Calder ([14]), respectively. Then runoff from each cell is modelled with a daily time step, considering vegetation, soil types and climate. The module responsible for this is largely based on the Probability Distributed Model (PDM, Moore [63]) and has been shown to be capable of producing realistic estimates of flows across Europe (Arnell [5]). Modelled cell runoff is routed through the river network including lakes, reservoirs, wetlands and artificial water transfers such as canals. For river channels this is done using a simple Muskingum method, whereas the routing formulation for other water bodies depends on their management and use. The amount of irrigation water required by crops is modelled on a monthly basis considering 24 different crop types across Europe and up to eight crop types per cell. This irrigation water requirement is modelled as the total crop-water requirement minus the part of rainfall taken up by the crop (effective precipitation). Total water requirement of irrigated crops is modelled according to FAO guidelines (Allen et al.[3]), using spatial data on irrigated areas and timing of growing seasons and using crop factors to represent crop growth. Effective precipitation is modelled according to USDA-SCS (1970 cited by USDA-SCS [78]), based on rainfall, crop evapotranspiration and soil water deficit. The remaining human water consumption is modelled considering annual-average population density, urbanization, livestock density (cattle, sheep, goats) and industry. The total human water consumption is then used to model total water abstraction by considering groundwater to surface water use ratios, return flows, water supply pipe leakage and irrigation efficiency (losses of abstracted irrigation water due to evaporation and percolation). It has recently been shown that GWAVA can well simulate runoff in natural catchments across Europe, compared to both measurements and other global hydrological models (Gudmundsson et al.[39]).

Graph: Fig. 1 Structure of GWAVA illustrated for a single grid cell (based on Meigh et al.[57]). Transfer of data to and from different cells is indicated with dashed arrows.

During this study, the GWAVA model has been improved as follows:

• a. the calculation of irrigation water use has been refined by distinguishing between more crops and crop categories than in previous GWAVA versions;
• b. the simulation of crop growth has been extended with the possibility to vary growth season length with climate; and
• c. the spatial resolution of GWAVA has been improved from 30′ to 5′, in order to make use of the higher-resolution data sets that are currently available for Europe.

Table 1 details the characteristics of the main inputs for GWAVA in this study.

Table 1 Input data for modelling water flows, their resolution, and source

Four parameters of the GWAVA rainfall–runoff module were calibrated using river discharge measured at 110 gauges whose catchments do not overlap and cover in total about half of the modelled area (Fig 2(a)). This resulted in different values for the four calibration parameters in each of these 110 catchments. The four calibrated parameters affect the lateral velocity of surface and sub-surface runoff, the sub-grid distribution of soil depth, and the calculation of soil field capacity and soil saturation capacity. The measured discharge data used for calibration were for the period 1980–2000 and were provided by the Global Runoff Data Centre (GRDC), Koblenz, Germany (http://www.bafg.de/GRDC) and the National Centre for Atmospheric Research (NCAR, Bodo [12]). The calibration method used was Downhill Simplex in Multidimensions (Nelder and Mead [64]), which automatically searches for the parameter values leading to the best fit, according to a user-specified measure of fit. In this study, we chose to let the calibration minimize the mean absolute error (equation (13)).

Graph: Fig. 2 Locations of gauging stations used for validation of modelled (a) water quantity and (b) water quality. Gauge locations are shown with black dots and their upstream areas are shown in grey.

Loads of TP, TN and BOD5 from land to rivers, lakes, reservoirs and wetlands are modelled as the sum of loadings from point sources and diffuse sources. Hereafter, we use "pollutant" to refer to any of TP, TN and BOD5.

Point-source loading of any pollutant is modelled on a 5′ resolution with a mass balance approach described in Williams et al. ([86], [87]). This approach distinguishes point sources arising from households, industrial discharges and runoff from paved urban areas. Pollutants from households are modelled using per-capita pollutant emissions, sewage treatment efficiencies (from Eurostat, EC [22], and Perry and Venderklein [65]), rural and urban population density (Table 1), and the fractions of the rural and urban population connected to sewage treatment works (from Eurostat, EC [22], and WHO/JMP [85]). Per-capita TN and TP emissions are taken from Grizetti and Bouraoui ([37]) after adjusting the per-capita TN emissions for international differences in protein levels of diets (from the FAOSTAT database), and after adjusting the per-capita TP emissions for international differences in consumption of detergents containing sodium tripolyphosphate (from Glennie et al.[36]). Per-capita BOD5 emission is from IPCC ([46]). Pollutants from industrial discharges are modelled using the spatial distribution of industry (from Flörke and Alcamo [32]), its return flow and the typical pollutant concentration in this return flow (from ICPDR [45], and a review by Williams et al.[87] of 45 literature references). In addition, removal of pollutants in treatment of industrial sewage is considered. Pollutants in runoff from paved urban areas are modelled using rainfall, urban area, population density and reported pollutant concentrations in urban runoff (from Mitchell [60]). We also modelled loading of pollutants from scattered settlements using the approach described in Williams et al. ([86], [87]). This is modelled similarly to loading from households, except that sewage treatment levels are different.

Agricultural diffuse loading of pollutants to surface water is derived from a calibrated export coefficient method in which measured annual average pollutant load at catchment outlets was regressed against catchment characteristics. The catchment outlets used were located at water discharge monitoring stations with nearby stations monitoring TP, TN, or BOD5 (EEA [26]). The explanatory variables considered were catchment area, cropland area, built-up area, livestock units, Köppen-Geiger climate, lake area, river channel length, runoff, temperature, slope, point source loading and fertilizer use of mineral and manure P and N and the atmospheric deposition of TN. Linear regression was performed using data from three time periods (1988–1992, 1993–1997, 1998–2002), during which only the significant (p < 0.05) explanatory variables were retained. Because there were no significant differences between the regression equations generated for the different time periods, the data for all time periods were combined to fit a single regression equation for each pollutant (number of catchments used is 79, 106 and 104 for TN, TP and BOD5, respectively). The result was an equation for area-specific loading at the catchment outlet for each of TP, TN and BOD5 having R2 = 0.79, 0.89 and 0.92, respectively:

(1)

Graph

Here, GaugeLoadj is the area-specific loading of pollutant j (TP, TN, or BOD5) at a regression gauge (kg km-2 year-1), R is runoff (mm year-1) in the gauged catchment, cR, j is the calibrated export coefficient belonging to R and ∑ci, jvari, j (i = 1, ... , m) are the remaining m terms in the linear regression equation. These m terms include the impacts of point source emissions and agricultural activities on pollutant loading at the regression gauge (Table 2). Further explanation of equation (1) is given in Malve et al. ([55]).

Table 2 Input data for modelling agricultural pollutant loading to surface water bodies

In this study we aim to model pollutant concentrations in a spatially-distributed manner. To do this, we first need an estimate of the pollutant loading to surface waters in each individual cell (as opposed to GaugeLoadj which represents a catchment average), and subsequently we route these pollutants through all downstream cells (as described in the next section) to estimate the pollutant concentrations for those cells. We estimated pollutant loading to surface waters by first applying equation (1) to each individual grid cell and then correcting the resulting value for the fact that pollutant loading to surface waters happens upstream of regression gauges. Therefore, the following correction was made to equation (1) when applying it to individual grid cells:

(2)

Graph

Here, CellLoadj is the loading of pollutant j to surface waters from diffuse and point sources in the same cell (kg km-2 year-1). To account for in-stream losses, it is assumed that pollutant loading from land to surface waters is, on average, a factor α higher than pollutant loading at the regression gauges which are further downstream. The method used to estimate αj is explained in Appendix A. The estimated values of αj for TN, TP and BOD5 are very close to 1 which indicates that the effect of αj is very small for these pollutants. However, when GWAVA is applied to much less stable pollutants with shorter residence times, then αj will be much larger than 1. Therefore, we argue that αj is an important parameter that should be retained in the model because it allows GWAVA to be used for a wide range of pollutants.

Equation (1) was further modified as follows:

(3)

Graph

Here, equals the long-term average river discharge leaving the cell (m3 s-1) if that river discharge is less than 1 m3 s-1. Otherwise is 1. Thus variable causes the pollutant loading from land to small (<1 m3 s-1) streams to become more correlated to their discharge. Variables S and account for the sewage generated in upstream neighbouring cells with low dilution capacities (<1 m3 s-1) that is discharged in the current cell because of its higher dilution capacity (kg km-2 year-1). This improves the model because the Water Framework Directive prescribes that EU countries place sewage discharges at locations where the receiving waters can dilute the effluent so that it does not harm the environment (e.g. EA [19]). We estimate that sufficient dilution capacity can be expected in rivers with a discharge above approximately 1 m3 s-1. This estimate is based on the water quality standards in Table 3, per-capita pollutant loadings to rivers estimated by Williams et al. ([87]) and the fact that most people in the UK and the USA are connected to sewage treatment plants serving, on average, around 13 800 people (derived from EA [20], and Michel et al.[58], respectively). Variable S in equation (3) is defined as follows:

(4)

Graph

Table 3 Critical levels used to indicate risk to biodiversity

Here, n is the number of upstream neighbouring cells with a river discharge below 1 m3 s-1 (0 ≤ n ≤ 8).

If the current cell and its upstream neighbours both have a discharge below 1 m3 s-1 then the current cell will discharge some of the sewage produced by its upstream neighbouring cells and at the same time it will redirect a fraction () of the sewage produced on its own area to downstream cells. Although this representation may differ from the reality in individual grid cells, in general the use of and S in equation (3) resulted in more realistic modelled pollutant loadings to water. Month-to-month variability in CellLoadj was obtained by using monthly values for R.

Pollutant concentration could be estimated by simply dividing modelled upstream pollutant loading through modelled water discharge. However, such an estimate would ignore the potentially large impact of spatial and temporal variability of drivers such as water residence time, surface water abstraction and accumulative pollutant loss. To account for this variability, we explicitly model transport of pollutants in surface waters for each individual grid cell across Europe.

The transport of pollutants after emission from point and diffuse sources is modelled by assuming that the pollutant is transported downstream with discharge through river reaches, lakes, wetlands, reservoirs and artificial water transfers. Whist being transported downstream, pollutant may leave the river network with water abstraction. In addition, it may be lost by sedimentation and transformation. The model assumes that the latter two removal processes have rates that are proportional to pollutant concentration. Modelled values of water discharge and water abstraction are calculated by the pre-existing GWAVA model. Volumes of water in lakes, wetlands, reservoirs and river reaches are used to convert pollutant loads to concentrations and to estimate the loss and travel time of pollutants in the river network. Surface water volumes are modelled using data on land surface morphology.

Here, we summarize the method used to calculate pollutant concentration. Assuming conservation of mass and complete mixing within each grid cell, the following differential equation was derived (see Appendix B for details):

(5)

Graph

Here, C is the pollutant concentration in a cell (kg m-3), t is residence time (s), Xin is the pollutant loading entering the cell (kg s-1), V is the surface water volume of a cell (m3), Qr is the river discharge leaving the cell (m3 s-1), Qa is the abstraction of water (m3 s-1), Qtr is the water outflow through artificial transfers (m3 s-1), p1 is a loss rate constant of the pollutant due to aquatic processes such as sedimentation and transformation (s-1) and p2 is a constant production rate of the pollutant (kg m-3 s-1) to ensure the net pollutant loss becomes zero when C reaches its pristine value. Equation (5) was solved analytically resulting in an estimate of C for every grid cell in each month of the modelled time period.

The pollutant loading into the grid cell (Xin) is calculated as follows:

(6)

Graph

Here, n is the number of upstream neighbouring cells (0 ≤ n ≤ 8); and are outgoing discharge (m3 s-1) and pollutant concentration (kg m-3), respectively, in neighbouring upstream cell i; is the water flux from incoming artificial transfer k (m3 s-1); and is the concentration in the grid cell where transfer k is coming from (kg m-3). Variable Xt is the pollutant loading from diffuse and point sources in the cell (kg s-1), which is modelled using equation (3).

The surface water volume of a cell (V) can comprise both river reaches (Vr) and an impoundment such as a lake, reservoir or wetland (Vl). Therefore, V is calculated as:

(7)

Graph

Volume Vl is calculated as follows:

(8)

Graph

Here, VL is the total volume of the lake, reservoir or wetland calculated for each time step by the pre-existing GWAVA; L is the total number of cells covered by the same lake, reservoir or wetland as the current cell; fland,i is the fraction of land in cell i that is not covered by the lake, reservoir, or wetland; and fland is the fraction of land in the current cell that is not covered by the lake, reservoir or wetland.

Variable Vr is calculated as:

(9)

Graph

Here d is river depth (m), w is river width (m), l is the river length without meandering (m) and fm is a meandering factor defined as actual river length divided by l. River width, w, is estimated using grid-cell discharge according to Allen et al. ([2]). Meandering factor, fm, is calculated as a function of grid-cell size according to Fekete et al. ([31]). River depth, d, is calculated according to Pistocchi and Pennington ([66]):

(10)

Graph

Here, s is river bed slope estimated from drainage network topography and the sub-grid cell (30′′ resolution) elevation distribution, 0.045 is a river bed roughness value representative for Europe (unitless) and is the average of river discharge entering and leaving the cell (m3 s-1).

The values of p1 (equation (5)) were calibrated with measured pollutant concentrations from EEA ([27]) for the period 1990–2000. This calibration was done manually by adjusting p1 to minimize the difference between the median of all modelled pollutant concentrations across the EU and the median of all measured concentrations across the EU. We calibrated p1 because literature values of pollutant loss rates are usually only validated for specific surface water types or expressed in incomparable units (Birgand et al.[11]).

Production rate p2 was estimated as:

(11)

Graph

Here, Cnat is an estimated natural background concentration for the modelled pollutant (kg m-3). Equation (11) is derived from equation (5) and the following three assumptions about cells in natural river systems: (a) inflowing water has the same natural pollutant concentration as outflowing water, (b) Xt = 0, and (c) concentration change over time is negligible (i.e. δC/δt ≈ 0). The first and second assumptions cause pollutant loading into the cell to equal the sum of the pollutant storage rate and the pollutant flux out of the cell (i.e. Xin = C(Qr + Qa + Qtr + δV/δt)). The consequence of this and the third assumption is that aquatic production of the pollutant (p2) equals aquatic loss of the pollutant (p1·Cnat). We acknowledge that the third assumption may ignore some seasonal variation. However, this is acceptable in this study, as the focus is on human-influenced systems with possible risk of eutrophication where most of the variation in C is driven by anthropogenic pollution and its dilution and transport. The Cnat that was used for BOD5 was the first percentile of the measured concentrations across Europe from the EEA Waterbase (EEA [26]) and the Cnat values for TN and TP were based on Smith et al. ([73]).

One of the aims of the present study is to indicate the human water security risk. For this we use an index that considers both the human demand for water and the availability of water. This index is GWAVA's Water Availability Index 4 (WAI4):

(12)

Graph

Here, avail10i is the multi-annual 10th percentile water availability (m3 s-1) for month i and demi is the multi-annual average human water demand (m3 s-1) for month i. Here, the term "multi-annual" indicates that the percentile or average is based on as many values as there are years in the modelled time period. WAI4 is between 0 and 1 in cells with sufficient water resources to satisfy the local demand for water. If WAI4 is between 0 and −1 then water scarcity is expected to occur regularly (more than once in 10 years) if there is no water supply infrastructure that can transfer water over large distances (further than the grid-cell radius). If such infrastructure is present, then a WAI4 between 0 and −1 indicates that future climate change or future water demand growth may result in the need to extend this infrastructure further into neighbouring regions, or to build new water supply reservoirs that can store excess water during wet months for release during dry months.

We calculated risk for aquatic biodiversity using European standards aimed at conservation of biodiversity. We assumed that biodiversity is at risk if BOD5, TN, or TP levels are above the critical levels prescribed by these standards. The critical levels that we used (Table 3) are developed to test compliance with the European Water Framework Directive (WFD) and Nitrates Directive. The UK Technical Advisory Group on the WFD (UK TAG) provides critical 90th percentile BOD5 levels for rivers and critical mean soluble reactive phosphorus (SRP) levels for lakes and rivers (WFD UK TAG [83], [84]). The UK TAG gives different critical levels for different types of lakes and rivers because the native biota living in these types have different sensitivities to BOD5 and SRP. We took the median over these different types and converted SRP levels to TP levels according to Bradford and Peters ([13]). The critical levels from the WFD UK TAG ([83], [84]), which we used for BOD5 and TP, are those that are exceeded in only 10% of "healthy" aquatic ecosystems, where even the most sensitive biological element (macro-invertebrates) is still undisturbed. Critical levels for TN were from Król and Sokół ([50]), who developed these levels for implementation of the EU Nitrates Directive. Their critical TN levels are those above which eutrophication and, thus, deterioration of biodiversity, starts to occur. The standards for stagnant waters (Table 3) were applied if more than 1% of the cell area is covered by lakes and wetlands according to the CCM2 database (Vogt et al.[81]). Otherwise, we applied the critical levels for running waters (Table 3). We only indicate risk for biodiversity where the modelled water volume per cell (equation (7)) is more than zero in all modelled months. Cells where this is not the case are unlikely to have point-source loading that is continuous, which would be a violation of one of our model assumptions.

All input data used in the modelling, except climatic input, represent the year 2000. Climatic input was monthly data from 1960 to 2000, of which the first 30 years were used for model warm-up. (Note: GWAVA downscales climate input to daily resolution for rainfall–runoff modelling.) The index of risk for biodiversity is simulated using month-to-month variability in hydrology and pollutant concentrations resulting from variability in climate input from 1990 to 2000. The human water security index (WAI4) is based on month-to-month variability in hydrology and human water demands in the period from 1970 to 2000.

Fit of modelled river discharge is indicated for 110 calibration gauges, both after and before calibration. The fit before calibration results from using default values for the four calibration parameters of the rainfall–runoff module in the 110 calibration catchments. This fit is indicative of the quality of modelled discharge in the catchments that were not calibrated. Default values for the four calibration parameters of the rainfall–runoff module are in the middle of the range of plausible values.

Fit of modelled pollutant concentrations is demonstrated in detail with plots of modelled vs measured concentrations for many locations, whereas fit of water discharge will be summarized using three indices of fit. The reasons for this more elaborate assessment of modelled pollutant concentrations is that pollutant modelling is a new feature of GWAVA that has not been published before, whereas validations of GWAVA-modelled water discharge have already been described in numerous other publications.

Model fit is indicated with three different measures, each indicating a different characteristic of the model fit. The first measure is the mean absolute error:

(13)

Graph

Here, is the mean absolute error (%), n is the number of used observations, obsi is an observation and modi is a model prediction. We use to indicate model fit because it is relatively sensitive to errors in predicted discharge during dry periods (Krauze et al.[49]). This is important, because our WAI4 index is sensitive to months with low discharge and because lower discharge can cause higher pollutant concentrations to which the biodiversity risk index is sensitive.

Secondly, we used the model bias:

(14)

Graph

Here μ denotes the arithmetic mean and the indices obs and mod indicate observed and modelled values, respectively.

Finally we quantified the model's ability to capture the temporal patterns of variability using Spearman's correlation coefficient r.

Figure 3 illustrates the capability of GWAVA to simulate river discharge for a large calibrated river basin, which is roughly in the centre of Europe: the Meuse River (near Lith, The Netherlands). We chose the Meuse to allow comparison with Fig. 4, which shows water quality in the Meuse. For other calibrated rivers in Europe, the fit of modelled discharge is given in Table C1 (Appendix C). The mean absolute error in these rivers is about 21% lower than in uncalibrated rivers (Table C1). Modelled discharge has a low fit for many rivers in Iceland, Norway and Finland. This is likely to be due to some simplifications in GWAVA's snowmelt module, especially the assumption that snowmelt is driven by temperature and but not by radiation. However, the main reason for this low fit is probably that the snowmelt module was not calibrated for this study.

Graph: Fig. 3 Comparison of modelled and observed monthly discharge in the Meuse River at GRDC station Lith.

Graph: Fig. 4 Modelled and observed annual average pollutant levels along the Meuse River around the year 1999: (a) BOD5; (b) TN; and (c) TP.

The spatial variability of modelled BOD5, TN and TP concentrations is compared to the spatial variability of measured concentrations in the Meuse River (Fig. 4). The Meuse was chosen for this comparison, as it has many suitable monitoring stations (EEA [26]), and because the values of many of its watershed properties (e.g. latitudes, slopes, point source emissions, agricultural intensities) are in the middle of the range of values that can be found in the rest of Europe. However, as in most river basins, there are not many stations with overlapping periods of measurement within the 1990–2000 period. For most stations in the Meuse, measured annual data were available for the years around 1999. Therefore, we only used measured concentrations from years close to 1999 for this validation of modelled spatial variability. The fit of modelled spatial variability in the Meuse was best for BOD5 () and worst for TP (). The general trend in the measured TN and TP concentrations along the Meuse (slight decrease from 300 to 3000 km2 upstream area followed by steep increase from 3000 to 30 000 km2 upstream area) is reproduced by the model (r = 0.91 and 0.82, respectively). The average concentration in the Meuse is modelled well for TP and BOD5 (Δμ = −6 and −7%, respectively). For all three pollutants, the spatial autocorrelation of the measured values is lower than that of the modelled values. This may be caused by the fact that measurement times do not exactly match between different stations along the Meuse, whereas the modelled concentrations for those stations represent exactly the same time period. In addition, the exaggerated spatial autocorrelation may stem from real-world local processes that are not included in the model.

The validation of modelled temporal variability of pollutant concentrations is shown for 24 stations (Fig 2(b)) selected from EEA Waterbase – Rivers version 11 (EEA [27]). The values reported in this database are measured concentrations which are aggregated seasonally or annually, although for most countries the database has only annual values. One station was selected for each country represented in the EEA database, except for countries lacking suitable stations. Within each country, the station with the largest catchment area and the highest number of measured seasons was chosen if this station had at least three years of data in the 1990–2000 period. Annual data values were only used if they were based on at least ten samples and seasonal values were only used if they were based on at least four samples. Some countries (e.g. Bosnia Herzegovina) had no station with at least three years of data values based on sufficient samples. Therefore, these countries are not used in our assessment of fit of modelled pollutant concentrations. The correlation between measured and modelled pollutant concentrations depends very much on the number of samples on which the measured concentrations were based. For example, relatively high correlations were found in station Kleve-Bimmen (Rhine River), where the number of samples per data value was very high (5–7). Concentrations near the mouths of large river basins are generally underestimated (5–7). This is likely to be due to the fact that parameter p1 was calibrated such that modelled concentrations in most grid cells have a low bias, whereas river reaches as large as those of the selected 24 water quality measurement stations (Fig. 2(b)) cover only a very specific category of the grid cells.

Graph: Fig. 5 Time series of modelled and observed BOD5 levels at the selected river water quality measurement stations. Annual average BOD5 levels are shown, unless the time-series header indicates that seasonal averages are shown. Triangles indicate modelled values and squares indicate measured values. Both modelled and measured values are aggregated over the same season or year.

Graph: Fig. 6 Time series of modelled and observed TN concentrations at the selected river water quality measurement stations. Annual average TN concentrations are shown, unless the time-series header indicates that seasonal averages are shown. Triangles indicate modelled values and squares indicate measured values. Both modelled and measured values are aggregated over the same season or year.

Graph: Fig. 7 Time series of modelled and observed TP concentrations at the selected river water quality measurement stations. Annual average TP concentrations are shown, unless the time-series header indicates that seasonal averages are shown. Triangles indicate modelled values and squares indicate measured values. Both modelled and measured values are aggregated over the same season or year.

The fit of modelled concentrations on the pan-European scale is not as good as can be expected on smaller scales (such as the scale of an individual catchment), where it is possible to parameterize more complex models using more detailed local information. However, such detailed data are not available on the scale of this study.

The fit of calibrated and validated river discharges is generally better than the fit of modelled pollutant concentrations. One reason is that modelling of pollutant concentrations relies on modelled hydrological variables, such as river discharge. Thus, any uncertainty in this modelled hydrology is added to the uncertainty in modelled pollutant concentrations. In particular, the relative uncertainty in dry-month discharge is important because pollutant concentration has a reciprocal relationship with dilution capacity. Another reason for the relatively low fit of pollutant concentrations compared to discharge is that measured concentrations in many river basins show decreasing trends during the modelled period. These decreasing trends, which are probably the result of quick-changing mitigation measures, are not reproduced by the model. For example, pollutant management has caused a marked fall in nitrate concentrations in Denmark, Germany and Latvia during the 1990–2000 period (EEA [25]). Modelled concentrations are also affected by model simplifications, such as those used to model pollutant loading from point sources (Williams et al.[87]). In this method, the most important simplifications are: (a) the use of a typical pollutant concentration in industrial return flow, (b) the assumption that this return flow has been treated in all countries, (c) the use of the same removal efficiencies for sewage treatment works of the same type in all countries (leading to errors especially because we do not know which sewage treatment works have nutrient removal), and (d) the assumption that households that are not connected to sewerage systems have good local treatment (secondary level).

Modelled concentrations are affected by uncertainties in the export coefficient equation used for diffuse pollutant loading (Malve et al.[55]). In this model, the most important uncertainties are: (a) that not all areas in Europe were equally represented when fitting the export coefficient equation, (b) that the number of water quality observations was low for some catchments used in this fitting (for some catchments as little as six measurements per year), and (c) that the export coefficient equation uses cropland area and livestock numbers as surrogates for agricultural input of mineral fertilizers and manure from livestock, respectively. The export coefficient method (equation (1)), on which pollutant loading to surface waters was based, has been validated by Malve et al. ([55]). They did this by developing an export coefficient equation for three different time periods (1988–1992, 1993–1997, 1998–2002) and three different sets of observed catchments, and showed that the resulting equations are not significantly different from each other, or from equation (1). This indicates that equation (1) is valid in time periods and locations that were not used for its calibration.

Finally, the another important source of uncertainty is the decay coefficient p1, which is spatially and temporally constant, resulting in the modelled patterns exhibiting less variability than the measured patterns (evident in 4–7).

The GWAVA model predicts that current risk for human water security is concentrated in Mediterranean, arid and densely populated parts of Europe (WAI4 < 0 in Fig. 8). The spatial pattern of water security risk appears to be mainly determined by population density and climate. In regions with warm climates having dry and warm summers (south of the Iberian peninsula, Greece, most of Turkey, Cyprus, Sardinia, Sicily, Israel and the Mediterranean coasts of Spain, France and Italy), water security risk is already indicated at moderate population densities (>25 inhabitants km-2) and in a few areas with very high irrigation water abstraction (>100 mm year, especially in the south of Spain). However, in regions further north, which have a more humid climate, water security risk is only indicated in areas that are predominantly urban (>425 inhabitants km-2). Examples of such areas can especially be found near the large cities in the west of Germany, The Netherlands and Belgium and in the southern half of the UK. In the extreme southeast of the modelled area, where arid climates dominate, water security risk is indicated even at extremely low population densities. Further, GWAVA indicates some risk for human water security in wetland-dominated parts of Scandinavia due to GWAVA's assumption that swamps cannot be used for water abstraction.

Graph: Fig. 8 Modelled human water security from 1990 to 2000 across Europe indicated by WAI4. WAI4 has a higher value if the available water quantity better fulfils the human demand for water.

The modelled human water security across Europe generally looks plausible, as it corresponds well with human water stress modelled on a similar scale by the WaterGAP model (Alcamo [1]). However, the water security risk may have been exaggerated in the wealthier regions across Europe, because our risk indicator does not include the financial and technical means of local water resources managers to extend their water supply infrastructure across cell boundaries, nor does it include their means to reduce water demand by measures such as hosepipe bans or water price increases. In the dryer parts of Europe, the water security risk may have been locally exaggerated because we did not use GWAVA's module for simulating large water transfers (e.g. transfers from large water-supply reservoirs to users in different cells) and artificial water courses across basin divides. The reason for this is the absence of a suitable data set of such structures on the pan-European scale. Moreover, this study may have exaggerated water security risk, because long-term storage of water in deep aquifers was not modelled (the model only has a shallow groundwater store) and thus some regions for which water security risk was indicated may in fact draw heavily on deep groundwater to compensate for lack of surface water. This is, for example, the case for Germany and Denmark (EEA 2005 cited by Furberg et al.[34]). However, regions that extensively exploit their groundwater run the risk of lowering their groundwater table, or of causing saltwater intrusion (as is often the case in the Mediterranean region). In such regions, the indicated water security risk is more a signal that there is an unsustainable rate of water abstraction that may cause actual water stress in the future. The last reason why water security risk may have been exaggerated is that our input data on reservoirs (Table 1) only include relatively large reservoirs. Thus, the effect of smaller reservoirs on reduction of water security risk may not have been accounted for.

Across Europe, high modelled pollutant concentrations can be found in areas with intensive agricultural activities, such as the Po Valley in Italy, and the lowlands of The Netherlands and Belgium. In contrast, low modelled concentrations occur in regions with little human influence, such as the Alps and Scotland. In addition, low pollutant concentrations are generally modelled in parts of Europe with large lakes and wetlands, such as Scandinavia, due to their long water residence times that increase decay of pollutants. Finally, the modelled pollutant concentrations may either increase or decrease in the downstream direction, depending on the degradability of the pollutant and the spatial patterns of hydrology and pollutant sources within the river basin.

Figure 9 shows that high TP concentrations are modelled especially in Israel, Belgium and Portugal. In Israel and Portugal, this is mainly due to high evapotranspiration rates leading to low river discharge, limiting the dilution of point-source pollutants, whereas the high TP concentrations in Belgium are mainly due to scattered settlements. High TP concentrations in many Eastern European rivers are modelled because of relatively low precipitation causing lower dilution capacity of these rivers. High TN concentrations are modelled along the Atlantic Ocean in Western and Southern Europe and in the Po basin. The reasons for high TN concentrations along the Atlantic Ocean in Western and Southern Europe are similar to those previously mentioned for TP, but, in addition, livestock is an important reason for high TN concentrations in The Netherlands, western France, western England and eastern Ireland. In the Po basin, high TN concentrations are mainly caused by emissions from cropland and scattered settlements. Figure 9 indicates that high BOD5 concentrations can be found in Belgium. This is due to high livestock densities and high emissions from scattered settlements. High modelled BOD5 concentrations are also common in Serbia Montenegro and Israel, due to high emissions from manufacturing.

Graph: Fig. 9 Modelled BOD5, TN and TP levels in surface waters aggregated over the period 1990–2000: (a) 90th percentile BOD5 level; (b) mean TN concentration; (c) mean TP concentration. Concentrations are only shown where surface water was present in all modelled months. Mapped pollution levels correspond to the critical levels used to indicate risk to biodiversity (Table 3).

The modelled concentrations across Europe look plausible given the spatial distribution of pollutant sources (mainly human population and livestock), dilution capacity (discharge) and sinks (mainly lakes, wetlands and big river reaches). Furberg et al. ([34]) showed the location of regions with relatively high and relatively low measured N and P concentrations throughout Europe. This matches well with some regions in Europe where GWAVA predicts that both TN and TP are relatively high (southwest of the Iberian peninsula, the Meuse and Scheldt basin, southern England). However, comparison with Furberg et al. ([34]) indicates that there are also some regions where GWAVA underestimates TN or TP concentrations (Bulgaria, Poland, Latvia). The modelled lake TP and river BOD5 levels (Fig. 9) are generally consistent with those reported in the Water Information System for Europe (WISE) (EEA [28]). The differences (higher modelled BOD in The Netherlands and lower modelled BOD in the south of Italy) can be explained by the fact that the WISE database is more representative of the situation after the year 2000, whereas our model results cover the period 1990–2000.

Modelled concentrations indicate aquatic biodiversity risk in about 35% of the European area, especially where lakes and wetlands are abundant (Fig. 10). Modelled biodiversity risk in about 67% of the affected areas is solely due to high TP concentrations. The reason why much of Scandinavia is modelled to have aquatic biodiversity risk (Fig. 10) is that lake cover exceeds 1% in about half of the Scandinavian cells. Therefore, the stricter TP standard for stagnant waters (Table 3) applies there. The number of cells where TN standards are not met is about 58% smaller than where the TP standards are not met and they largely coincide with cells where TP standards are not met. The number of cells where BOD5 standards are not met is about 87% smaller than where the TP standards are not met, although their spatial distribution is different. In contrast to TP, BOD5 standards are not met in large proportions of the Caucasus region and in Serbia-Montenegro, making BOD5 the main cause of aquatic biodiversity risk in those regions.

Graph: Fig. 10 Locations where the model indicates risk for human water security and aquatic biodiversity. Risk for aquatic biodiversity is indicated if the TP, TN, or BOD5 level is above the standard given in Table 3. Risk for human water security is indicated if WAI4 (equation (12)) is negative. Risk for aquatic biodiversity is only shown where surface water was present in all modelled months.

Our result that aquatic biodiversity risk is mostly caused by too high TP concentrations is corroborated by numerous other published studies indicating that TP is the most common driver behind freshwater eutrophication (e.g. Guildford and Hecky [40]). The combined spatial pattern of indicated risk for biodiversity and human water security agrees with local studies identifying similar risks (SCENES [70]) and with spatial modelling of similar risks done by Vörösmarty et al. ([82]).

Generally, the model results look good and plausible. Validation showed that the developed modelling method can predict the spatial variation of modelled concentrations. However, the temporal variability of concentrations is only modelled well in a proportion of the river basins. Therefore, future model improvements should focus on this modelled temporal variability. GWAVA modelled concentrations and flows have not been corrected with measurements, thus allowing model errors to propagate downstream. Such corrections were not applied in order to allow the model to be applied for prediction of the future for which measured concentrations and flows are obviously not available.

When interpreting the indicated biodiversity risk and water security risk, it should be kept in mind that no indicator can capture all factors that affect these risks. Instead the indicators are intended to provide a broad picture of risk, which is based on a few very important drivers that are likely to be dominant in many locations.

We suggest that the enhanced GWAVA version presented in this paper could benefit river basin management plans (RBMPs), which are implemented throughout the European Union as prescribed by the Water Framework Directive. The new GWAVA version can be used to assess whether the water quality objectives of RBMPs are sustainable in the long term given long-term changes in the climate, economy and population distribution. Also, it allows assessment of possible conflicts that may arise between water quality objectives in RBMPs and the water needs of the human population.

The fact that the net pollutant loss in surface waters only depends on parameters p1 and Cnat makes it feasible to parameterize the model for pollutants whose fate and behaviour in the environment is not yet very well known, as is often the case for newly-developed chemicals. This makes our model potentially useful for environmental risk assessments of new chemicals (ECA [23]).

We developed a spatially and temporally explicit model of human water security and water pollution. This model represents a substantial proportion of measured spatial and temporal variability in measured river discharge and levels of TN, TP and BOD5 in rivers on the pan European scale. However, there is scope for improvement, especially in the representation of the temporal variability of concentrations. Deviations of model results from published statistics of TN, TP and BOD5 levels can be explained by uncertainty in these statistics due to low sampling frequencies. Another reason is the relatively low complexity of the water quality module to allow for parameterization of many different pollutants on a scale as large as the pan-European scale. In general, however, the model results look good and plausible.

The developed model can account for changes in a comprehensive set of driving forces, including drivers such as climate, population, damming, cropping patterns and commitments to wastewater treatment. It thus enables integrated water resources managers to better assess the effects of anticipated changes.

The results of the developed model show that water security risk and biodiversity risk are likely where the population density is high, the agriculture is intensive, the climate is dry and lakes or wetlands are abundant.

Future use of the presented method will include modelling of future scenarios and additional pollutants.

We thank the SCENES project from the European Commission (FP6 contract 036822) for funding the research that resulted in this paper. Also, we thank GRDC for providing most of the data needed for calibration of modelled river discharge.

The values of αj (equation (2)) were estimated by first assuming that the loss of pollutants in the river network of most regression catchments is dominated by aquatic loss processes such as sedimentation and transformation. Secondly, we assumed that this river network loss can be approximated using an analogy with a long unbranched channel consisting of many well-mixed subsections, such that the loss in each of these subsections can be described as:

(A1)

Graph

Here, Cs is the pollutant concentration in a subsection (kg m-3), t is residence time (s), and p1 is a loss rate constant of the pollutant due to aquatic processes such as sedimentation and transformation (s-1). A solution of equation A1 can be expressed as follows:

(A2)

Graph

Here, Δt is the time needed for the pollutant to travel through a subsection (s), x is the number of the current subsection (starting at 1 and increasing in downstream direction until xmax), is the pollutant concentration in the river water flowing out of subsection x (kg m-3), and is the pollutant concentration in the river water flowing into subsection x (kg m-3). Variable is calculated as:

(A3)

Graph

Here, PollLoad is the pollutant loading from diffuse and point sources into each subsection (kg s-1), and Qs(x) is the river discharge in subsection x (m3 s-1). We assume that PollLand is the same for all subsections and that Qs(x) is proportional to x. Both Qs(x) and PollLand are assumed to be constant in time. Using these assumptions, we can iteratively apply equations (A2) and (A3), starting at x = 1 and ending with xmax. This results in the value of which represents the concentration at a regression gauge. The values of xmax and Δt are chosen such that:

(A4)

Graph

The value chosen for tmax is 67 000 s, which is the expected travel time of the pollutant to the corresponding regression gauge. This value is based on cell travel times estimated as the ratio of modelled V to modelled Qr (defined below equation (5)) for individual cells in the regression catchments.

Equations (A2), (A3) and (A4) allow us to estimate αj as:

(A5)

Graph

Here, the numerator represents CellLoadj and the denominator represents GaugeLoadj. The estimated value of αj does not depend on the value of PollLand and the ratio of Qs(x) to x. Thus it only depends on the value chosen for p1.

Derivation of equation (5) is based on a similar approach used in the derivation of the QUESTOR model (Eatherall et al.[21]). The following equation describes the mass balance of a pollutant in a grid cell while assuming complete mixing:

(B1)

Graph

Here, C is the pollutant concentration in a cell (kg m-3), V is the surface water volume of a cell (m3), t is residence time (s), Xin is the pollutant loading entering the cell (kg s-1), Qr is the river discharge leaving the cell (m3 s-1), Qa is the abstraction of water (m3 s-1), Qtr is the water outflow through artificial transfers (m3 s-1), p1 is a loss rate constant of the pollutant due to aquatic processes such as sedimentation and transformation (s-1) and p2 is a constant production rate of the pollutant (kg m-3 s-1) to ensure the net pollutant loss becomes zero when C reaches its pristine value.

From the product rule follows that:

(B2)

Graph

Combining equations (B1) and (B2) gives:

(B3)

Graph

Rearrangement of equation (B3) gives equation (5).

Table C1 Model fit in 110 gauges having non-overlapping catchment areas. Fit is expressed using (%), Δμ (%) and Spearman's correlation coefficient r for model runs with and without calibrated parameters (indicated with "Calibration" and "Validation", respectively)