The changing trend in nitrate concentrations in major aquifers due to historical nitrate loading from agricultural land across England and Wales from 1925 to 2150 the

on torical nitrate loading from agricultural land on the evolution of groundwater nitrate concentrations. An addi-tional process-based component was constructed for the saturated zone of signi ﬁ cant aquifers in England and Wales. This uses a simple ﬂ ow model which requires modelled recharge values, together with published aquifer properties and thickness data. A spatially distributed and temporally variable nitrate input function was also in- troduced. The sensitivity of parameters was analysed using Monte Carlo simulations. The model was calibrated using national nitrate monitoring data. Time series of annual average nitrate concentrations along with annual spatially distributed nitrate concentration maps from 1925 to 2150 were generated for 28 selected aquifer zones. The results show that 16 aquifer zones have an increasing trend in nitrate concentration, while average nitrate concentrations in the remaining 12 are declining. The results are also indicative of the trend in the ﬂ ux of groundwater nitrate entering rivers through base ﬂ ow. The model thus enables the magnitude and timescale of groundwater nitrate response to be factored into source apportionment tools and to be taken into account alongside current planning of land-management options for reducing nitrate losses.


H I G H L I G H T S
• An approach to modelling groundwater nitrate at the national scale is presented. • The long time-lag for nitrate in the groundwater system is considered. • The impact of historical nitrate loading on groundwater quality is better understood. • Areas of high groundwater nitrate input to surface water are delineated. • The method is transferable and requires a modest parameterisation.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Nitrate (NO 3 ) is essential for living matter by acting as a source of nitrogen (N) that forms the building blocks of molecules. Man has benefited from the application of chemical nitrogen fertilisers to gain increased agricultural productivity. However, any excess nitrate applied can be leached from the soil into freshwater, depending on the timing, rate and method of fertiliser application. Nitrate concentrations in groundwater beneath agricultural land can be several-to a hundred-times higher than that under semi-natural vegetation (Nolan and Stoner, 2000). Nitrate-contaminated water can cause long-term environmental damage and threaten both economic and ecosystem health (Bryan, 2006;Defra, 2002;Mayer et al., 2002;Pretty et al., 2000;Sebilo et al., 2006;Thorburn et al., 2003;Ward, 2009).
Agricultural land is the major source of nitrate water pollution (Ferrier et al., 2004;Thorburn et al., 2003;Torrecilla et al., 2005). In England, over 70% of nitrate in groundwater and surface water has been shown to be derived from agricultural land (Foster, 2000;Hunt et al., 2004;Defra, 2006). Although legislation has been introduced to reduce nitrate water pollution, this remains an international problem (Campbell et al., 2004;European Environment Agency, 2000;Ferrier et al., 2004;Mayer et al., 2002;Rivett et al., 2007;Sebilo et al., 2006;Thorburn et al., 2003;Torrecilla et al., 2005;Wang et al., 2013;Wang and Yang, 2008;Yang and Wang, 2010;Yue et al., 2014). The EU Nitrate Directive, an integrated part of the EU Water Framework Directive (Directive 2000/60/EC), is implemented through a statutory instrument that sets rules (Action Programme rules) for best agricultural practices within Nitrate Vulnerable Zones (NVZs). These are designated areas where land drains and contributes to the excess nitrate found in contaminated groundwater and surface waters. However, the degradation of freshwater quality due to nitrate remains a problem in the UK. In many areas, nitrate concentrations are more than 50 mg NO 3 L −1 with a rising trend in both rivers (Burt et al., 2008 and aquifers (Smith, 2005;Stuart et al., 2007). However, evidence of some improvements in groundwater nitrate (e.g. Smith et al., 2010) is beginning to emerge.
Storage of nitrate in porewater and consequent slow vertical migration through the unsaturated zone of major aquifers was first recognised in the 1970s (Foster and Crease, 1974;Young et al., 1976). The importance of this for nitrate management was highlighted by Dautrebande et al. (1996). The anticipated decrease in nitrate concentrations in the aquifer following the implementation of protection measures was not observed.
In the unsaturated zone (USZ), i.e. from the base of the soil layer to the water table, pore-water pressure is sub-atmospheric, and hence the fractures and matrix of rock may be only partially saturated. Therefore, hydraulic pathways within, and between rock matrix and fractures may be restricted (Ireson et al., 2009). In addition, the water and pollutants held by capillary tension on fracture walls could be an important means of storing soluble pollutants in the USZ (Price et al., 2000;Sorensen et al., 2015). Thus, the historical nitrate loading could, in some cases, stay in the USZ for a long time before reaching the saturated zone or aquifer. Through recent research, it has become increasing clear that it could take decades for leached nitrate from the soil to discharge into freshwaters (Jackson et al., 2007;Burt et al., 2011;Howden et al., 2011;Wang et al., 2012a;Allen et al., 2014). This may cause a long time-lag between improvements in agricultural practices or targeted land-use change and the reduction of nitrate concentrations in interconnected groundwater and surface water (e.g., Burt et al., 2008;Allen et al., 2014). Current environmental management strategies rarely consider the nitrate time-lag, but rely instead on the predictions of a relatively rapid response of water quality to land-management practices (Smith and Lerner, 2008;Collins et al., 2009;Burt et al., 2011). Therefore, models or tools are needed to simulate nitrate storage and time-lag in groundwater system and thus support better-informed water resources management.
Catchment-scale models have been developed to simulate water and nitrate transport in both the USZ (e.g., Jackson et al., 2007;Ireson et al., 2009) and the saturated zone (e.g., Harbaugh et al., 2000;Jackson et al., 2006). Most of these models are catchment specific as a wide range of factors affects nitrate transport and fate in the groundwater system. At a regional or national scale, it is necessary to use a more generic methodology with an appropriate level of conceptual complexity to predict changes in average behaviour. Wang et al. (2012a) developed the national-scale parsimonious process-based "Nitrate Time Bomb" (NTB) model to produce maps showing whether the peak nitrate loading has arrived at the water table for different aquifers in Great Britain. The results helped to understand the trend in the amount of nitrate arriving at the water table and entering groundwater due to the historical nitrate loading from arable land. It was also useful for groundwater resource management to designate areas where the historical nitrate burden in the USZ is high. However, the NTB model developed for that original work needed to be extended to include a saturated zone component to make it useful for wider applications, such as estimating trends in groundwater nitrate concentrations. Another limitation of the original NTB model is that it used a single historical nitrate input function.
The purpose of this study was therefore to simulate the long-term trend in nitrate concentrations in the major aquifers of England and Wales, based on an extended NTB model which incorporates: • a national-scale hydrological conceptual model of nitrate transport and dilution in groundwater • estimates of nitrate velocity in the USZ based on readily available recharge values and aquifer properties • a spatially distributed nitrate input function which reflects the historical agricultural loadings at different locations on the ground surface.

Methodologies
The NTB model simulates the nitrate transport in the USZ and estimates the time and the amount of historical nitrate arriving at the water table (Wang et al., 2012a). This used nitrate concentrations in porewater from 300 cored boreholes, and nitrate velocity estimates in USZs from both field data for the major aquifers of the Chalk, Permo- Triassic sandstone and Jurassic limestones and literature values for other lithologies. Surface lithologies were taken from 1:625,000 scale hydrogeological mapping. This model has been extended in the following ways to investigate the impacts of historical nitrate loading from arable land on the changing trend in groundwater nitrate concentrations.

Nitrate transport velocity in the USZ
Amongst the key factors affecting pollutant velocity of transport in USZs, the groundwater recharge rate, aquifer porosity and the storage coefficient are important (Leonard and Knisel, 1988). The nitrate velocity in the USZ and hence the residence time can be expressed as (Rao and Davidson, 1985;Rao and Jessup, 1983): where Thickness USZ , i is the thickness of USZ at cell i ( Fig. 1); V USZ , i (m year −1 ) is the nitrate-transport velocity in the unsaturated zone; q i (m year −1 ) is groundwater recharge at cell i; Rf aquifer (−) is the retardation factor determined in the calibration procedure; and Sr aquifer (−) is the specific retention for the rock. The specific retention represents how much water remains in the rock after it is drained by gravity, and is the difference between porosity and specific yield.

Introducing distributed historical and projected nitrate loadings from agricultural land
The single nitrate-input-function derived by Wang et al. (2012a) reflects historical and future agricultural activity from 1925 to 2050. It has been validated using mean pore-water nitrate concentrations from 300 cored boreholes across the UK in the British Geological Survey (BGS) database (Stuart, 2005). A low nitrogen loading rate (25 kg N ha −1 year −1 ) between 1925 and 1940 reflects the pre-war low level of intensification with very limited use of non-manure-based fertilisers. The gradual intensification of agriculture during, and just after, World War II resulted in a 1 kg N ha −1 year −1 rise in nitrogen input to 40 kg N ha −1 by 1955. A more rapid rise of 1.5 kg N ha −1 year −1 between 1955 and 1975 was due to increases in the use of chemical-based fertilisers to meet the food needs of an expanding population. The nitrogen input declines from 70 kg N ha −1 in 1991 to 40 kg N ha −1 in 2020 at a rate of 1 kg N ha −1 year −1 as a result of restrictions on fertiliser application. Finally, for the future (2020-2050), it was assumed that there would be a return to nitrogen input levels similar to those associated with early intensive farming in the mid-1950s, i.e. a constant 40 kg N ha −1 nitrogen loading rate (Wang et al., 2012a). However, this single-input-function only generated a national average, rather than a spatially distributed input based on the evolution of agricultural activity across England and Wales.
NEAP-N (Anthony et al., 1996;Environment Agency, 2007, Lord andAnthony, 2000;Silgram et al., 2001) is a meta-model of the NITCAT (Lord, 1992) and NCYCLE (Scholefield et al., 1991) models, with adjustments for climate and soil type (Anthony et al., 1996). NEAP-N has been used for policy and management in the UK by regulatory agencies, such as the Environment Agency and the Department for Environment, Food and Rural Affairs (Defra). It includes a water-balance model and a leaching algorithm. Nitrate loss potential coefficients are assigned to each crop type, grassland type and livestock categories within the June Agricultural Census data to represent the short-and long-term increase in nitrate leaching risk associated with cropping, the keeping of stock and the spreading of manures. The model predicts the total annual nitrate loss from agricultural land across England and Wales and the associated water flux (hydrologically effective rainfall). For this study, the NEAP-N loss potential coefficients used were revised for each of the prediction years (1980, 1995, 2000, 2004 and 2010, corresponding to years with full agricultural census data for farms across England and Wales). This is to account for changes in nutrient applications (fertiliser and manure), crop yields and livestock yields (meat or milk) over time. The predicted NEAP-N nitrogen loadings (1 km by 1 km) for these years were used in this study to derive the temporally and spatially distributed nitrate-input-functions. NEAP-N predictions have formed the basis of much strategic policy-based work on linking agricultural emissions to water quality (e.g. Zhang et al., 2014). The trend of historical nitrogen loading in the original single nitrate-input-function of NTB was used to interpolate and extrapolate the data for the years other than the NEAP-N years. This enabled a nitrogen loading map to be calculated for each year from 1920 to 2050. Fig. 1 shows two examples of the new derived nitrate-input functions for locations in 'Chalk, Southern England' and 'Carboniferous Limestone, South Wales and South West England'.

A national-scale conceptual model for nitrate transport and dilution in aquifers
In order to simulate nitrate transport and dilution processes in aquifers at national scale for England and Wales, a simplified hydrogeological conceptual model was developed. The groundwater system in England and Wales can be conceptualised as an island system on the basis of the following assumptions: • groundwater recharge supplies water to aquifers as an input, and no lateral water flow in USZs were considered in this national-scale study • groundwater flows out of the system through rivers in the form of baseflow as an output on an annual basis. The seasonal interaction between aquifer and rivers has been ignored in this national-scale modelling study with an annual time-step, and it was assumed that nitrate in aquifers is only from the USZs • groundwater is disconnected from rivers where low permeability superficial deposits are present • the total volume of groundwater (Vol total ) for an aquifer varies from year to year due to the change of groundwater recharge. Vol total in a simulation year is the sum of the groundwater background volume (Vol background ) and the annual groundwater recharge reaching the water table (Vol rech arg e ). Groundwater recharge and baseflow reach dynamic equilibrium whereby the amount of recharge equals that of baseflow in each of aquifer zone in a simulation year. It is assumed that Vol background remains same in each simulation year.
• nitrate entering an aquifer is diluted throughout the total volume of groundwater in a simulation year • the velocity of nitrate transport in aquifers is a function of aquifer permeability, hydraulic gradient and porosity • the transport distance for groundwater and nitrate can be simplified as the total distance between the location of recharge and nitrate entering the aquifer at the water table and their discharge point on the river network (Fig. 2).
2.3.1. Groundwater available for nitrate dilution Aquifers are discretised into equal-sized cells for numerical modelling purposes. The total volume of groundwater Vol total (t) (m 3 ) in a simulation year t is calculated using the equations: where A i (m 2 ) is the area of cell i; D aquifer (m) is the depth of active groundwater (for nitrate dilution) in an aquifer; Sy aquifer (−) is the specific yield representing the aquifer drainable porosity; and n is the total number of cells in the aquifer.
Observations show that the water table responds to recharge events at the surface on a time-scale of days or months, while the residence time for pollutant fluxes in the USZ is in the order of years (Headworth, 1972;Lee et al., 2006;Wang et al., 2012a). This can be explained as a 'piston-displacement' mechanism (Headworth, 1972;Price et al., 1993). Water and nitrate are displaced downwards from the top of the USZ. Therefore, instead of travelling through the USZ, water and nitrate reaching the water table are displaced from the bottom of the USZ. On this basis, the volume of recharge entering an aquifer Vol recharge (t) (m 3 ) in the simulation year t can be expressed using the equation: where t is time (year); Rp q (year) is the water-table response time to rainfall events; and q i (t − Rp q ) (m year −1 ) is the annual recharge at cell i at time t − Rp q .

The velocity of nitrate transport in aquifers
The average velocity of nitrate transport in an aquifer VS mean (m year −1 ) can be calculated using the following equations when an aquifer cell does not overlap with a river cell: where T aquifer (m 2 day −1 ) is the transmissivity of the aquifer; G i (−) and Dist i (m) are, respectively, the hydraulic gradient and horizontal distance between cell i and the nearest point where groundwater is discharged into the river; D aquifer (m) is the depth of active groundwater in an aquifer; Φ aquifer is the porosity of an aquifer zone; GWL i (m) is the groundwater level for cell i; RL i (m) is the river level at the nearest river point to cell i; and VS i (m year −1 ) is velocity for cell i. For those aquifer  Table 1. cells that spatially overlap with river cells, the nitrate travel time in them is assigned to be zero without using the Eqs. (6)-(8).

Annual nitrate concentration in groundwater
Annual nitrate concentration Con aquifer (t) (NO 3 mg L −1 ) in year t can be calculated aquifer by aquifer, assuming there is no groundwater flow between aquifers: where RTime total, i (year) is the total residence time for nitrate to travel through both the USZ and an aquifer at cell i (Fig. 1); RTime USZ,i (year) is the nitrate residence time at cell i in USZ; RTime SZ,aquifer (year) is the average residence time for nitrate dilution and transport in the aquifer; M i (t − RTime total ,i ) (mg NO 3 ) is the amount of nitrate loading from the base of soil into USZ at cell i in the year of t − RTime total , i ; and ATT is the attenuation factor representing the percentage of nitrate mass that is attenuated in the USZs.

Model construction and data
Twenty-eight zones were used to define the important unconfined aquifers in England and Wales for this study. These were selected by focusing on those with N10-year nitrate time-lag in the USZ and a baseflow index (BFI) N 0.4. BFI is the average ratio of annual baseflow The groundwater level for cell i RL i (m) The river level for cell i ATT (−) the nitrate attenuation factor in the USZ Thickness USZ,i The thickness of USZ at cell i Monte Carlo calibration Φ aquifer (−) The porosity for an aquifer zone Sy aquifer (−) The specific yield for an aquifer zone Rf aquifer (−) The retardation factor for calculating the nitrate velocity in USZs T aquifer (m 2 day −1 ) The transmissivity for an aquifer zone D aquifer (m) Depth of active groundwater for an aquifer zone to annual river flow in a catchment and represents how well aquifers and rivers are connected. The waters in these aquifer zones ( Fig. 3 and Table 1) are, therefore, those more likely to be affected by historical nitrate loading in the longer term than other areas. The digital 1: 250,000 hydrogeological mapping of Great Britain from the BGS (2015) and BFI derived from the hydrology of soil types (HOST) classification scheme (Boorman et al., 1995) were used as the basis for identifying these key aquifer zones. The subdivision of large aquifers, such as the Chalk, into a series of areas reflects the scale at which published data on aquifer properties are available in Allen et al. (1997). The different aquifers are normally separated by aquitards and generally different zones of a specific aquifer are geographically separated. Although 'Chalk, Thames' and 'Chalk, East Anglia' are connected with each other, their groundwater-flow direction (from northwest to southeast) is almost parallel to the boundary between the two zones. This means that the groundwater interaction between them is limited. Therefore, it was assumed that there is no flow interaction between aquifer zones in this national-scale study. The study areas were discretised into 1 km by 1 km cells. Cross-correlation is a time series technique to evaluate the statistical correlation between two sets of data as a function of the lag of one relative to the other. It has been used to reveal the significance of the water-table response to rainfall after a given time; and it can also allow the time taken for the first water-table response to rainfall to be calculated (Lee et al., 2006;Mackay et al., 2014). The cross correlation analysis between the rainfall and groundwater-level time series was performed in this study to estimate the water-table response time Rp q to rainfall events. The datasets used in this calculation include the time series of monthly rainfall , from the Meteorological Office Rainfall and Evaporation Calculation System (MORECS) (Hough and Jones, 1997) and groundwater level for 57 boreholes across the study area. Rp q was set to the period of time over which there is a correlation between groundwater level and rainfall at more than 95% confidence level. The average value of Rp q was calculated in the aquifer zones where there is more than one borehole (Table 1).
A national-scale groundwater recharge model was built using the soil water balance model SLiM (Wang et al., 2012b), which objectively estimates recharge and runoff using information on rainfall intensity, potential evapotranspiration, topography, soil type, crop type, and BFI. The model was calibrated using the observed river-flow data for 102 gauging stations (http://www.ceh.ac.uk/data/nrfa/). The long-termaverage and time-variant recharge estimates were then used to simulate nitrate-transport velocity in the USZ and the groundwater volume Vol total (t), respectively.
The distance to river points for each cell was calculated in ArcMap™ using the CEH digitised (1:50,000) river network (Moore et al., 1994). Aquifers are disconnected from rivers in the areas where lowpermeability superficial deposits are present (see http://www.bgs.ac. uk/products/digitalmaps/dataInfo.html). These areas were, therefore, masked out when calculating the distance to river. The hydraulic gradient was calculated using long-term average groundwater levels (Wang et al., 2012a) and river levels derived from the gridded Digital Surface Model (NextMap DSM). The aquifer properties of active groundwater depth, porosity, transmissivity and specific yield were based on the collation of Allen et al. (1997).

Calibration
Monte Carlo (MC) simulations were undertaken to calibrate the revised model. Parameters were randomly sampled within a finite parameter range to produce one million parameter sets. The upper and lower bounds of the range for each parameter were defined based on observed results or expert judgement. Performing MC simulations is a computerintensive task especially when multiple parameters are involved. Therefore, it is good practice to reduce the number of parameters for MC simulations by fixing some parameters using available information on the aquifer zones. All parameters used in this study are summarised in Table 2. The fixed parameters were identified based on existing datasets and hydrogeological knowledge from hydrogeologists.
Two sets of MC simulations were conducted to calibrate the model against: 1) the nitrate velocity values in USZ derived from measurements of porewaters from cores from bored boreholes (Wang et al., 2012a), and; 2) the observed average nitrate concentrations for each aquifer zone calculated from monitoring data provided by the Environment Agency.
In the former, the bias (absolute difference) between simulated and observed nitrate velocity in USZs was used to evaluate the model fit. In the latter, the Nash-Sutcliffe efficiency (NSE) score (Nash and Sutcliffe, 1970) was adopted to calculate the goodness-of-fit between observed and modelled nitrate concentrations, via: where Vobs i is the observed nitrate concentration at the ith time-step; Vsim i the simulated nitrate concentration at the ith time-step; Nis the total number of simulation time-steps; and Vobs is the average value of observed nitrate concentration in N simulation times. A negative NSE score indicates that the observed mean is a better predictor than the modelled results; a value of zero denotes that modelled data are considered as accurate as the mean of the observed data; and a value of one suggests a perfect match of modelled to observed data. The model with the highest score in a set of MC simulations is deemed to have the optimum parameter set.

Sensitivity analysis
A sensitivity analysis of the model parameters was undertaken to determine which parameters contribute most to the model efficiency, and which of these parameters are identifiable within a specific range linked to known physical characteristics of an aquifer zone. Scatter plots for parameter values against the biases or NSE scores from MC simulations were produced. By representing a MC run, each dot in the scatter plots shows the model performance of this MC run in the vertical axis when using a parameter value in the horizontal axis. In each scatter plot, the grey dots represent many MC runs; the upper surface of this cloud of points represents a response surface that indicates how the model performance changes as each parameter is randomly perturbed. A flat response surface, such as the scatter plot of specific yield for 'Upper Greensand' in Fig. 4, means that the model performance does not change much when the parameter value changes; this indicates that the model is not sensitive to this parameter and the optimum parameter value is not identifiable. The black dot denotes the optimum parameter value that generates the best modelled result. Fig. 4 shows some examples of the scatter plots in estimating nitrate velocity values in USZs using specific yield, porosity and the retardation factor. Although the sensitivity of the model to these parameters differs for each aquifer zone, in general, the model is most sensitive to the retardation factor and least sensitive to specific yield. The models for 'Chalk, S England', 'Upper Greensand' and 'Corallian, Yorkshire' show clear V-shaped response surfaces for the retardation factor, indicating that this parameter is identifiable although there is more than one value with a bias close to zero. The optimum parameter values result in the minimum bias in the MC simulations. In contrast, the response surfaces for specific yield are nearly flat in these three aquifer zones and do not show a unique optimum. The response surfaces for porosity show that the model is sensitive to this parameter to some extent. Fig. 5 shows some examples of the scatter plots of the NSE scores against depth of active groundwater and transmissivity in the second set of MC simulations. The model is sensitive to both the depth of active groundwater and the transmissivity to different extents; and these parameters are identifiable for the different aquifer zones.

Estimates of nitrate concentration
The annual nitrate concentrations for the 28 selected aquifer zones in England and Wales were simulated based on the calibrated model for the period 1925 to 2150. This uses spatially distributed recharge values for 1961 to 2011. Outside this period, the long-term average recharge value  was used, resulting in less fluctuation in the modelled nitrate concentrations from year to year. Table 1 lists the values of the Root Mean Square Error (RMSE) between modelled and observed nitrate concentrations for the 28 selected aquifer zones. A smaller RMSE means better fit between modelled and observed data. RMSE values range from 0.2 mg NO 3 L −1 to 6.8 mg NO 3 L −1 with an average value of 4.1 mg NO 3 L − 1 . The results are acceptable for such a national-scale study focussing on long-term trend in nitrate  concentrations in aquifers. Fig. 6 shows some examples of the modelled time series, which had well defined trends in the observed data. Table 1

Discussion
A conceptual model was developed in this study to simulate the nitrate transport and dilution processes in groundwater. Its purpose is to assess the long-term trend of average nitrate concentrations in aquifer zones at the national scale. For zones where the observed data defined a positive upwards trend, model results show good agreement with the observed data (Fig. 6). The results showed that 12 out of 28 aquifer zones have a declining trend in nitrate concentration and the remainder will have a rising average nitrate concentration until the concentration turning-point year is reached (Table 1).
However, a number of limitations should be borne in mind when interpreting the results generated by the updated model. First, it is assumed that nitrate is transported only by intergranular movement through the matrix in the USZ. Taking the Chalk as an example, Smith et al. (1970) carried out an experiment measuring the tritium content of pore water from the Chalk USZ at a Berkshire site. Their findings suggested that 85% of the total flow in the USZ was by intergranular seepage through the matrix. Prior to this, however, it was widely believed that the flow and its solutes moved predominantly through fractures in dual-porosity media such as the Chalk (Headworth, 1972). However, there is little or no data for other dual-porosity aquifers and this has not been included in this model.
There is also the question of nitrate loss in the unsaturated zone. According to Close (2010) and Environment Agency (2005), nitrate is negatively charged and thus electrostatically repelled by media in USZs that usually have a negative charge, such as clay minerals. This means that nitrate is less likely to be sorbed within USZs. Moreover, denitrification, which is generally facilitated by the absence of oxygen, is considered to be the dominant nitrate attenuation process in the subsurface system (Rivett et al., 2007). However, denitrification was found to be relatively limited in both USZs and unconfined aquifers selected in this study (e.g., Butcher et al., 2005;Kinniburgh et al., 1994;Lawrence and Foster, 1986;Rivett et al., 2007). Therefore, the nitrate attenuation in the groundwater system has been ignored in this study, but the attenuation factor ATT can be parameterised when more information on nitrate attenuation in each aquifer zone is available. Moreover, the study included some karstic formations, i.e. the 'Carboniferous Limestone' (zones 14 (part), 15-17), where there is likely to be limited, if any, intergranular flow, with virtually all nitrate being transported through the fractures. The differences in groundwater level response time and nitrate turning point within this aquifer in different zones could possibly be due to geological differences. The limestone contains basalt horizons in Derbyshire and significant amounts of mudstone and sandstone present in the sequence in northern England.
Since agreed future climate change scenarios were unavailable for this study, the annual average recharge, which has been calculated based on annual recharge estimates from 1961 to 2011, was used to simulate nitrate concentrations in the future. Consequently, nitrate concentration values at a turning-point after 2011 should be treated with caution and are simply indicative of the most likely peak average nitrate concentrations. Better simulations can be performed when future groundwater recharge values under different climate change scenarios become available.
Long-term average groundwater levels were used in this study to calculate the depth of the USZs. The year by year groundwater volume was calculated as the sum of the time-variant groundwater recharge and the groundwater background volume Vol background . It was assumed that Vol background remains the same over the years due to the lack of information on time-variant groundwater flow at a national scale. When national time-variant groundwater-level data become available in the future, this model can be improved to consider time-variant Vol background and time-variant USZs' depth to generate better estimates of nitrate concentration.
Finally, it is assumed in this study that there is no recharge and hence nitrate entering the groundwater system where low- permeability superficial deposits are present. However, experiments show that water and nitrate can travel through these deposits, and the amounts are determined by the thickness (e.g. Butcher et al., 2009). The impacts of low permeability deposits can be considered in the revised model presented here in the future. In addition, no lateral migration of nitrate in USZs was considered in this study. All these simplifications, however, facilitated the modelling procedures for this national-scale study that aimed to update the NTB model to provide a framework for better informing land-management strategy needing to take better account of nitrate legacy issues.
The trend in nitrate concentrations of major aquifers in this study is useful to support surface water-quality management. Groundwater is essential for maintaining the flow of many rivers in the form of baseflow. The results in this study can be used to study the impacts of the long-term groundwater nitrate concentration change on river quality, especially for rivers connected with high permeability aquifers. Since the aquifer zones selected in this study have baseflow index values of N 0.4, the nitrate held up in the groundwater system can greatly affect surface water quality where rivers are connected to these aquifers, and hence the ecological quality. On this basis, some of the outputs described herein have been incorporated into ongoing updates to the SEPARATE (SEctor Pollutant AppoRtionment for the AquaTic Environment) screening tool for England and Wales (see Zhang et al., 2014 for version 1.0).

Conclusions
This paper presents an approach to modelling groundwater nitrate at the regional or national scale in England and Wales. It requires relatively modest parameterisation and runs on an annual time-step but still provides useful estimates of present and future average groundwater nitrate concentrations. These results will help decision makers to understand how the historical nitrate loading from agricultural land affects the evolution of water quality due to the long time-lag for nitrate in the groundwater system. This model will be particularly valuable for evaluating the long-term impact and timescale of land-management scenarios and programmes of measures introduced to help deliver waterquality compliance (e.g. Zhang et al., 2014). It also delineates areas of high groundwater nitrate input to surface water which would obscure the impact of, for example, farm-scale interventions. However, these modelled annual average nitrate concentrations should be carefully examined before they are used to solve localised nitrate problems. More complexities could be introduced into future NTB type models to simulate the detailed regional nitrate transport processes in groundwater.