OR/16/036 Development of the BGS unsaturated zone model
|Stuart, M E, Wang, L, Ascott, M, Ward, R S, Lewis, M A, and Hart, A J. 2016. Modelling the groundwater nitrate legacy. British Geological Survey Internal Report, OR/16/036.|
|Nitrate units of measurement|
Nitrate concentrations can either be quoted as nitrate or NO3 (mg NO3 l-1) or as the equivalent amount of N (mg NO3-N l-1). These units are converted using the factor (molecular weight NO3/molecular weight N = 62/14 =4.43. For example the drinking water limit for nitrate is 50 mg NO3 l-1 or 11.3 mg NO3-N l-1. In this report mg NO3-N l-1 is used except where existing graphical material uses other units. Nitrogen in fertilisers is quoted as kg N ha-1.
This chapter describes a series of new developments to the BGS NTB model (Wang et al, 2012a, 2012b). In the original model the distribution of nitrate arriving at the water table depended on only three functions: the nitrate input at the land surface, the rate of travel of nitrate through the unsaturated zone and the thickness of the unsaturated zone (Figure 2.1).
The unsaturated zone thickness and nitrate velocity are combined to estimate the spatial distribution of nitrate travel time in the unsaturated zone and from this the input year for nitrate reaching the water table at any defined time. A nitrate input function over time can then be used to estimate the concentration reaching the water table at any point and defined time, assuming that nitrate is conservative.
These changes described in the following section address a number of simplifications used in this first version of the model including:
- A distributed nitrate input function to replace the original single function.
- A revised unsaturated zone thickness using Ordnance Survey (OS) river data instead of assumed river locations from the Nextmap DTM.
- Travel time attribution using the 1:250 000 geological mapping rather than the 1:625 000.
- Estimating nitrate velocity in the unsaturated zone for key aquifers using a process-based model replacing single values for lithological types.
- Introducing nitrate transport processes in low permeability superficial deposits. Nitrate transport to the underlying aquifer was originally assumed to be zero.
The improvements introduced as part of this project should allow the national model to be applied more effectively at the sub-national scale. This section also includes the first estimate of the mass of nitrate stored within the unsaturated zone across the UK, made possible as a result of the improved model.
Nitrate input function
The BGS NTB model previously used a single nitrate input function (NIF) (Wang et al, 2012a, 2012b) providing a national average rather than a spatially distributed input based on agricultural activity. This reflected historical and future agricultural activity from 1925 to 2050 and was validated using mean pore-water nitrate concentrations from 300 cored boreholes across the UK in the BGS database (Figure 2.2).
A low nitrogen loading rate between 1925 and 1940 reflects a pre-World War II low level with very limited use of non-manure-based fertilizers. The gradual intensification of agriculture during and just after the war 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 fertilizers 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 with a rate of 1 kg N ha-1 year-1 as a result of restrictions on fertilizer application under the Nitrate Directive. Finally, the model assumes 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 from 2020 to 2050 (Wang et al., 2012, 2012b). The model could readily be adapted to incorporate any agreed forward look (scenarios).
The historical (and future) nitrate loading information was developed into a temporally- and spatially-distributed NIF using modelled nitrate leaching data from NEAP-N (Anthony et al., 1996; Environment Agency, 2007; Lord and Anthony, 2000; Silgram et al., 2001). NEAP-N 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). 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 agricultural census data to represent the short and long-term increase in nitrate leaching risk associated with the cropping, keeping of stock and the spreading of manures. The model predicts the total annual nitrate loss from agricultural land for 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 available — 1980, 1995, 2000, 2004 and 2010 (corresponding to years with full agricultural census data for farms across England and Wales) — 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 distributed nitrate-input-functions (NIFs). This enabled a nitrogen loading map to be calculated for each year from 1920 to 2040. Figure 2.3 shows examples of the derived nitrate-input functions for locations in the principal aquifers of the Chalk of South England and the Permo-Triassic sandstone of Yorkshire and the Coal Measures of Wales and the Forest of Dean.
Results and implications
Figure 2.4 shows the impact of using a temporally and spatially-distributed NIF. The most significant spatial differences are for non-agricultural and low productivity areas, such as moorland, in Northern England and Wales. Areas, such as the Coal Measures, which have a considerable unsaturated zone thickness due to their deeply incised nature, may be of low agricultural productivity and therefore not significant with respect to unsaturated zone nitrate.
The temporal distribution of nitrate applications is also clearly shown in Figure 2.4. The impact of high nitrate applications peaking in 1980–1990 is most obvious in East Anglia, North Yorkshire, Wessex and the Welsh borders.
The new composite NIF function provides a credible estimate of nitrogen applications that could be made back to the early 20th century to provide a source term for very long groundwater flowpaths and can also be projected forwards using agreed scenarios to assess future impacts of nitrogen management measures.
Water level mapping
Updated water level mapping
The first (original) version of the NTB model included the use of water levels derived from assumed river locations from the Nextmap DTM as well as autumn minimum groundwater level values available from hydrogeological maps. A new version of groundwater level mapping which uses OS river data has been applied in this study.
Results and implications
Figure 2.5 shows the impact of these changes on estimated unsaturated zone travel time. The unsaturated zone is considerably thinner in some areas e.g. in parts of the Chalk and this results in reduced travel times.
The use of this new version of groundwater levels is considered to be much more realistic as these levels control the extent of flushing of nitrate from the unsaturated zone in each year.
Reclassification of geology
The original model used geology attributed at the 625k scale corresponding to the most recently issued hydrogeological mapping. A project objective was to improve hydrogeological representation by using more detailed geological map data and ideally the 50k scale maps. An evaluation of the very large number of formations, which would need to be considered if the 50k maps were to be used, suggested that this would be extremely onerous and disproportionate to the anticipated benefits. Instead it was decided to use 250k scale maps. These are of a large-enough scale to allow the layering of rock strata to be addressed for example in the Coal Measures and the Jurassic limestones. There are a few instances where the 250k maps do not appear to correspond well lithologically with the 625k scale maps and in these cases a comparison with the 50k scale mapping was made to resolve the differences and incorporate the necessary refinements.
Attribution of aquifer properties
In addition to the estimated travel times used in the original model, a range of additional aquifer properties were attributed from Allen et al. (1997) and Jones et al. (2000) for selected aquifer zones to provide parameters which would be required by the processed-based model (Table 2.1).
|Aquifer zone name||Principal/
|Zone area (km2)||Water table response time to recharge (months)|
|Chalk, S England||P||4,965||8|
|Chalk, East Anglia||P||5,603||20|
|Chalk, NE England||P||3,304||13|
|Lower Greensand, Bedford-Cambridge||P||440||11|
|Lower Greensand, Weald||P||1,183||13|
|Corallian, S England||S||773||11|
|Millstone Grit, Cumbria, Durham and Northumberland||S||2,730||11|
|Permo-Triassic Sandstone, Lancashire-West Midlands||P||4,701||26|
|Permo-Triassic Sandstone, SW England||P||957||8|
|Hastings Beds, Weald||S||2,235||11|
|Carboniferous Limestone, N England||S||5112||11|
|Carboniferous Limestone, N Wales||P||442||11|
|Carboniferous Limestone, S Wales and SW England||P||893||1|
|Carboniferous Limestone, Derbyshire||P||448||4|
|Permo-Triassic Sandstone, Nottingham-N Yorkshire||P||2,953||11|
|Middle Jurassic limestone, Cotswolds-Dorset||P||2,196||7|
|Middle Jurassic limestone (excluding Lincolnshire Limestone), Lincolnshire-Oxfordshire||P||2,451||9|
|Middle Jurassic, Yorkshire||S||961||11|
|Coal Measures, Pennines||S||3,235||11|
|Coal Measures, S Wales and Forest of Dean||S||2,258||11|
|Coal Measures, Durham-Northumberland||S||1,584||11|
|Millstone Grit, S Pennines||S||4,318||11|
|Magnesian Limestone, Durham||P||635||11|
|Magnesian Limestone, Nottingham–N Yorkshire||P||888||11|
The water table response shown in Table 2.1 was calculated by cross-correlating groundwater levels and rainfall for each aquifer zone. The value was set to the period of time over which there was a correlation with a greater than 95% confidence level. Further details are given in Wang et al. (2016).
The original NTB model principally focussed on the major and minor aquifers as these were considered most likely to have a significant unsaturated travel time. These zones were addressed using\process-based model described in the next section.
Table 2.2 summarises the statistic information on the old and new datasets of the USZ thickness and nitrate velocity in the USZs. Although the maximum USZ thickness is half of the old one, the mean value is slightly smaller than the old average USZ thickness. The calculated nitrate travel time in the USZs is 1.6 times higher than the old one on average. This implies that the newly derived nitrate velocity is lower than the previous one.
The new nitrate travel time in the USZ across the England and Wales was calculated using the more detailed geological information and the new unsaturated zone thickness map (Figure 2.5) and nitrate velocity derived from the 250k hydrogeological map. Figure 2.6 shows its comparison to the old map derived in the previous study.
|Old USZ thickness (m)||0||200||13.7||21.4|
|New USZ thickness (m)||2||100||13.2||16.4|
|Old nitrate travel time in the USZs (m/year)||0||434.8||15.7||23.9|
|New nitrate travel time in the USZs (m/year)||0.2||973||25.8||40.3|
Building a national-scale conceptual model
The first stage of building a national-scale process model is to develop a simple conceptual model:
- Nitrogen from arable land and from livestock is applied to the surface of the soil; this is oxidised to nitrate and a proportion is leached from the base of the soil.
- Nitrate is transported in infiltration through the unsaturated zone at a rate controlled by recharge and individual aquifer properties.
- Nitrate reaches the saturated zone and is transported by water movement towards surface water at a rate controlled by the groundwater gradient and by aquifer properties.
Nitrate reaches surface water as baseflow or other surface receptors 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 (Figure 2.7). 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. Groundwater flows out of the system through rivers in the form of baseflow as an output. Groundwater recharge and baseflow reach dynamic equilibrium whereby the amount of recharge equals that of baseflow in a simulation year. This implies that the volume of groundwater in an aquifer remains the same at the beginning and end of each simulation year, and is called the background volume (Volbackground). The total volume of groundwater (Voltotal) for an aquifer in a simulation year is the sum of the background volume and the annual groundwater recharge reaching the water table (Volrech arg e).
- 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, and.
- The transport length 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.
Developing a national recharge model for the UK
Hydrological processes in the recharge model
A national-scale groundwater recharge model was built using the soil water balance model SLiM (Wang et al., 2012a), which objectively estimates recharge and runoff using information on rainfall intensity, potential evapotranspiration, topography, soil type, crop type, and Base Flow Index (BFI). Figure 2.8 illustrates the principal soil zone hydrological processes associated with the SLiM method.
There are many potential pathways that water can take through the system. Rainfall could, in part, be intercepted by plants, while the remaining rainfall reaches the ground surface and infiltrates the soil to reduce soil moisture deficit (SMD). Soil water is extracted by plant roots for transpiration or drawn to the bare soil surface for evaporation. When soil moisture reaches field capacity (SMD becomes zero), water drains freely from the saturated soil, and the additional water added to the soil system, called the excess water of SMD, can flow laterally overland as runoff (including surface runoff and interflow) if a slope gradient exists towards adjacent locations, or percolate downwards through saturated soil as recharge. Runoff flows to nearby areas, called run-on, can join in the soil hydrological processes at the new location. SLiM explicitly derives both runoff and recharge based on the calculated SMD and other datasets, such as BFI and slope, instead of expert judgment. Figures 2.8 and 2.9 show how these are represented in the SLiM model.
BFI, representing an average ratio of annual baseflow to annual river flow in a watershed or catchment, is strongly related to topographical, soil, hydrogeological and precipitation characteristics and less influenced by the land-cover properties of a catchment. BFI is one of the key characteristics in explaining the hydrological processes and allows separation of the total river flow into baseflow (slow flow through the groundwater system) and surface flow (fast flow through the overland process) portions (Haberlandt et al., 2001). The equations for calculating runoff and recharge in a cell are as follows:
RECH = ESMD • BFI 
Ro = ESMD • (1–BFI) 
where ESMD (mm) is the depth of SMD excess water when soil becomes free draining; RECH (mm) is the depth of water that move downwards to recharge groundwater system; and Ro (mm) is the runoff calculated from ESMD.
Many studies (see for example Haan et al., 2006; Lange et al., 2003; Tani, 1997) show that the topographic gradient is an important factor that controls runoff generation. In general, greater runoff is observed in areas with steeper slopes. As mentioned above, BFI is a long-term average ratio reflecting different catchment characteristics including average catchment slope. Therefore, the runoff and recharge calculated using equations 1 and 2 can be understood as the averages generated at the locations with an average slope in a study area. If a location has a higher than average slope, greater runoff and reduced recharge will be generated. In a flat area, where a zero slope gradient exists towards neighbouring locations, all SMD excess water becomes recharge and no runoff will be generated. The relationship between runoff and slope is shown in Figure 2.8. Equation 2 can be further formulated as:
where Slpmean (degree) is the average slope value in a catchment; and Slp (degree) is the slope value at a cell with the area.
Compared to longer less intense rainfall, high intensity short duration rainfall is more likely to exceed the capacity of the soil to infiltrate water and generate more overland flow. Tani and Abe (1987) show that rainfall intensity and antecedent soil water storage in a forested catchment affect the amount of runoff and, if the rainfall intensity is larger than a threshold (e.g. 100 mm day-1), the increase in storm runoff is almost the same as the increase in rainfall, even with dry antecedent soil water conditions. Therefore, a rainfall intensity threshold is introduced in SLiM, to represent bypass runoff, where the amount of rainfall above this threshold becomes runoff. If the rainfall intensity is less than this threshold, the SMD excess water method is used for calculating runoff. The rainfall intensity threshold needs to be calibrated.
Recharge model construction and calibration
The recharge model has the same spatial resolution as the other components of the BGS NTB model -1 km x 1 km — giving a total area of 229 619 km2. In contrast to the NTB model, it uses a daily time step over the period 1962–2011. It included the catchments for all 102 gauging stations providing observed river flows (www.ceh.ac.uk/data/nrfa/)
Monte Carlo (MC) simulations were carried out in calibrating the model against the river flow data. The Nash-Sutcliffe efficiency NSE (Nash and Sutcliffe, 1970) was adopted to calculate the goodness of fit between observed and modelled surface flow time series:
where Vobsi is the observed surface flow at the ith time step; Vsimi the simulated flow at the ith time step; N is the total number of simulation time steps; and Vobs is the average value of observed flow in N simulation times.
In general, a negative NSE indicates that the observed mean is a better predictor than the modelled results. Where NSE is zero modelled data are considered as accurate as the mean of the observed data, and NSE between zero and one can be treated as acceptable levels of model performance. The closer to 1, the more accurate the model is and an NSE of one corresponds to a perfect match of modelled to observed data (Nash and Sutcliffe, 1970).
Field experiments (e.g. Butcher et al., 2009) showed that the thickness of low permeability superficial deposit affects the amount of water and soluble pollutants (such as nitrate) entering the groundwater system. Traditionally the thickness of superficial deposits has seldom been considered in simulating groundwater recharge. Other factors such as changes in composition within the deposits and fracturing are also important but were not addressed here as it is difficult to model these at a national scale. Enhanced recharge at the margins of low permeability deposits was also not included.
A method was developed in this study to objectively estimate recharge considering the spatial distribution and thickness of low permeability superficial deposits. These were divided into five thickness classes, i.e., 0–2 m, 2 m–5 m, 5 m–10 m, 10 m–30 m, and >30 m. The reduction of recharge for each class was identified based on MC simulations (Figure 2.10).
- 0–2 m: 0.095%
- 2–5 m: 4.23%
- 5–10 m: 9.10%
- 10–30 m: 15.3%
- 30–50 m: 99.99%
Recharge model calibration results
The model showed an acceptable performance for 73% of the surface water gauging stations used for model testing/calibration with a NSE >0.5 (Figure 2.11). Twelve percent of results had NSE <0.5 which is considered reasonable performance for a simple model running at the national scale. Since 15% of the gauging stations have the problem of missing values for observed flow data, they were excluded when calibrating the national recharge model.
Two examples, i.e. Craigiehall and Llanfair, were randomly chosen from many results to demonstrate how the modelled hydrographs match with observed ones (Figure 2.12). It shows that modelled results are in line with observed hydrographs with NSE values larger than 0.5 and the model defines the episodes of higher flows very well but does not always estimate the transient high peaks.
Figure 2.13 sets out a series of recharge parameters over the annual cycle for a point randomly selected within the model. This shows rainfall to be distributed over the year but there is little or no recharge over the summer and autumn due to evapotranspiration and development of an increased soil moisture deficit (SMD). There are a few runoff events throughout the year. This is the typical pattern as recharge occurs generally during the winter when plant growth is minimal in a temperate climate (Rushton et al., 2006).
The outputs from this model, the long-term-average (Figure 2.14), daily (e.g., Figure 2.14a) and yearly (e.g., Figure 2.14b) recharge estimates, were then used to simulate nitrate-transport velocity in the USZ and the groundwater volumeVoltotal (t), respectively.
Estimating nitrate transport velocity in the unsaturated zone
In the original version of the model nitrate velocities were estimated from measured values for three of the most important aquifers in the UK; the Chalk, Permo-Triassic sandstone and the Lincolnshire Limestone together with literature values and professional judgement for the others.
In order to estimate velocities a number of aquifer properties are required. Recharge, aquifer porosity and storage coefficient/specific yield are major factors affecting pollutant transport in the USZ (Leonard and Knisel, 1988). The nitrate velocity in the USZ and hence the residence time can be expressed as (Rao et al., 1985; Rao and Jessup, 1983):
where ThicknessUSZ,i is the thickness of USZ at cell i (Figure 1); VUSZ,i (m year-1) is the nitrate-transport velocity in the unsaturated zone; qi (m year-1) is groundwater recharge at cell i; Rfaquifer (-) is the retardation factor determined in the calibration procedure, and; Sraquifer (-) 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. Model calibration is described in Estimating annual nitrate concentrations in groundwater.
Recharge can be modelled and aquifer property data, such as porosity and specific yield are available for major and minor aquifers from Allen et al. (1997) and Jones et al. (2000). However these specific yields have been estimated for the saturated zone and assume that the aquifer remains fully saturated. These may not adequately represent conditions in the unsaturated zone.
Estimating nitrate transport velocity and dilution in the saturated zone
Model calibration and estimation of errors are described in Estimating annual nitrate concentrations in groundwater.
Nitrate transport velocity can be estimated using the following equations:
where VSi (m year-1) is the nitrate velocity for cell i, Taquifer (m2 day-1) is the transmissivity of the aquifer, Gi (-) is the hydraulic gradient between a cell and its nearest river discharge point, Dist (m) is the distance between a cell and its nearest river discharge point.
Groundwater available for dilution of nitrate in the saturated zone can be estimated using the following equations:
where Ai (m2) is the area of cell i; Daquifer (m) is the depth of active groundwater (for nitrate dilution) in an aquifer; Syaquifer (-) is the specific yield representing the aquifer drainable porosity, and; n is the total number of cells in the aquifer model; Rpq (year) is the water table response time to recharge events;
It can be argued that recharge events at the bottom of the soil effectively occur at the same time as rainfall events within a monthly time step (Lee et al., 2006; Mackay et al., 2014). Therefore, Rpq was determined using the cross-correlation between the time series of monthly rainfall (1961–2011), from the Meteorological Office Rainfall and Evaporation Calculation System (MORECS) (Hough and Jones, 1997) and groundwater level for 57 boreholes across the study area. Rpq, the time taken for an instantaneous flux of recharge to reach the water table, was set to the period over which there is a significant correlation at the 95% confidence level. Figure 2.16 gives an example of cross-correlation results for the borehole at Tank Hall in the Chalk, East Anglia, which is randomly chosen from many results, and shows that groundwater level and rainfall event are correlated (larger than 95% confidence level) for 24 months. This indicates that 24 months are needed for the groundwater level to fully respond to the monthly rainfall event. The average values of Rpq were calculated for aquifer zones where there was more than one borehole and these are given in Table 2.1.
Annual nitrate concentration in groundwater
Annual nitrate concentration Conaquifer(t) (NO3 mg L-1) in the SZ in year t can be calculated aquifer by aquifer, assuming there is no groundwater flow between aquifers:
where RTimetotal, i (year) is the total residence time for nitrate to travel through both the USZ and an aquifer at cell i (Figure 2.7); RTimeUSZ, i (year) is the nitrate residence time at cell i in USZ; RTimeSZ,aquifer (year) is the average residence time for nitrate dilution and transport in the aquifer; Mi (t– RTimetotal,i) (mg NO3) is the amount of nitrate loading from the base of soil into USZ at cell i in the year of t–RTimetotal,i, and; ATT is the attenuation factor representing the percentage of nitrate mass that is attenuated in the USZs.
Nitrate transport in low permeability superficial deposits
As mentioned in Developing a national recharge model for the UK, some water and nitrate can pass through low permeability superficial deposits (drift). This process was ignored in the original NTB model to reduce the number of parameters in the national-scale model as it was previously assumed that no nitrate can be transported downwards through the low permeability superficial deposits. However, nitrate transport through (and below) these deposits has now been introduced to address the previous over simplification. This has been possible as a result of newly available datasets, such as the split between runoff and recharge derived from a national-scale groundwater recharge model.
This model was built using the soil water balance model SLiM (Wang et al., 2012a) which objectively estimates recharge and runoff using information on rainfall intensity, potential evapotranspiration, topography, soil type, crop type, and BFI (Wang, 2015; Wang et al., 2014). The model was calibrated using the observed river-flow data for 102 gauging stations (www.ceh.ac.uk/data/nrfa/). The long-term-average and time-variant recharge estimates were then used to simulate nitrate-transport velocity in the USZ and the groundwater volume, respectively.
The revised low permeability superficial deposits model addresses three aspects:
- The nitrate loading from superficial deposits now depends on their thickness and is consistent with groundwater recharge estimates.
- The nitrate travel time in superficial deposits is calculated using velocity from literature and the deposit thickness.
- The nitrate is transferred to the underlying bedrock aquifer and is then subject to unsaturated and saturated travel time and dilution process in the saturated zone.
The thicker the low permeability superficial deposits are, the less water and nitrate can pass through them and enter the groundwater system. This is consistent with the concept of the splitting of runoff and recharge from excess water of SMD in Developing a national recharge model for the UK. The reduced recharge is then diverted into surface water to increase runoff, and the same portion of nitrate in it is also flushed into surface water (Figure 2.17). Therefore, the amount of nitrate entering the groundwater system is a function of BFI, topography and the thickness of superficial deposits. All modelling cells in the selected aquifer zones were included to simulate nitrate dilution in groundwater regardless of the presence of overlying low permeability deposits. The model does not address enhanced recharge at the margins of these deposits.
The low permeability superficial deposit zones were identified using the classifications of Griffiths et al. (2011) and comprised the Low-Medium and the Low-Low zones. The average nitrate velocities derived from the work of Ross et al. (1989), Klinck et al. (1996) and Marks et al. (2004) were used for this modelling.
Based on the previously calibrated NTB model, the initial results of introducing nitrate processing in low permeability superficial deposits (before the new calibration) suggest that areas where impact is significant are the:
- Chalk of East Anglia
- Chalk, Northern
- Lower Greensand of Bedfordshire & Cambridge
- Permo-Triassic Sandstone of Lancashire–West Midlands
- Carboniferous Limestone of Northern England
- Permo-Triassic Sandstone of Nottingham–N Yorkshire
- Middle Jurassic limestone of Lincolnshire–Oxfordshire (excluding Lincolnshire Limestone)
- Coal Measures of Durham–Northumberland
- Magnesian Limestone of Nottingham–N Yorkshire.
Figure 2.18 shows examples of the impact of taking account of nitrate transport through low permeability superficial deposits above two aquifer zones. It is shown that modelled nitrate concentration peak values decreased and the time to the peak has been delayed after introducing the nitrate transport in low permeability superficial deposits before undertaking calibration (Figure 2.18b). This might be explained by longer nitrate travel time through low permeability deposits and more groundwater for dilution after considering modelling cells in aquifers overlain by low permeability deposits. After model recalibration (Figure 2.18c) a good fit with the observed data is obtained, but for the Chalk of East Anglia the modelled peak is enhanced and broadened.
Estimating annual nitrate concentrations in groundwater
Model calibration (MC) simulations were also undertaken to calibrate the improved NTB 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 judgment. Performing MC simulations is a computer-intensive task especially when many 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.3. The fixed parameters were identified based on observations and knowledge from hydrogeologists.
qi (m year-1)
NIF (kg ha-1)
|The area for cell i |
The recharge value for cell i
The water table response time to recharge events
The groundwater level for cell i
The river level for cell i
The nitrate attenuation factor in the USZ
The thickness of USZ at cell i
|Monte Carlo Calibration||Φaquifer (-)
Taquifer (m2 day-1)
|The porosity for an aquifer zone|
The specific yield for an aquifer zone
The retardation factor for calculating the nitrate
The transmissivity for an aquifer zone
Depth of active groundwater for an aquifer zone
Two sets of MC simulations were conducted to calibrate the model against: 1) the nitrate velocity values in USZs derived from measurements of porewaters from drill cores (Wang et al, 2012a, 2012b), 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 NSE score was adopted to calculate the goodness-of-fit between observed and modelled nitrate concentrations.
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, to show how the model efficiency changes as each parameter is randomly perturbed. Figure 2.19 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, Southern England, Upper Greensand and Corallian Limestone 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.
Figure 2.20 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 and these parameters are, to different extents, identifiable for the different 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. Examples are shown in Figure 2.21.
The revised model still depends on a number of assumptions and these could limit its use for some applications. These are:
- Nitrate transport is only by intergranular movement through the matrix in the USZ,
- - It does not fully take account of bypass flow in transporting nitrate rapidly to the water table. Bypass flow is likely to be important in all fractured and karst aquifers, e.g. the Corallian and Carboniferous limestones as well as the Chalk, particularly where there is very low matrix permeability.
- - The impact of diffusion will be accounted for I in measured values but is not otherwise included.
- Nitrate attenuation in the USZ was ignored. This was because:
- - Nitrate is negatively charged and will not be affected by cation exchange (Close, 2010; Environment Agency, 2005).
- - Rivett et al. (2007) suggested that denitrification only accounts for a loss of 1–2% in the unsaturated zone. In the saturated zone reaction rates are generally limited by a lack of electron donors (most often dissolved organic carbon and/or sulphide) so most denitrification occurs only in the confined aquifer, where depleted in dissolved oxygen; exceptions are where there is elevated DOC due to the presence of pollutants or surface water infiltration. However, although nitrate may diffuse into the small pore throats of the Chalk and Jurassic limestones, even if organic carbon and sulphides are present and oxygen is absent, denitrification may not occur, as nitrate reducing bacteria are excluded by the small pore throats.
- - Denitrification was found to be relatively limited in the unconfined aquifers selected in this study (e.g. Butcher et al., 2005; Kinniburgh et al., 1999).
- Climate and long-term average values for recharge in the future (2011–2150) will be the same as the recent past (1961 to 2011).
- The future NIF is extrapolated from the composite NIF based on the BGS NIF and ADAS NEAP-N for 1980, 1995, 2000, 2004 and 2010. For realistic projections a future NIF, or a series of scenarios, which reflect anticipated nitrate management needs to be agreed.
Total unsaturated zone nitrate storage
Management of legacy nitrate in groundwater in England and Wales both at national and local levels requires an understanding of the storage of nitrate in the unsaturated zone. Work by Wang et al. (2012b) has identified the peak year for nitrate to arrive at the water table. Parris (1998) and Worrall et al. (2009) show that Great Britain is a net sink of reactive nitrogen, with potentially 300 kt N stored in groundwater. However, estimations of the total mass of nitrate in the unsaturated zone have not been undertaken to date. This section details an approach that has estimated the total mass of nitrate in the unsaturated zone of England and Wales.
In this high level national scale quantification of nitrate storage in the unsaturated zone it was deemed suitable to select high and moderate productivity aquifers based on BGS 1:625 000 scale hydrogeological mapping. Areas overlain by low permeability superficial deposits (Griffiths et al., 2011) were excluded from the analysis. The time and spatially variable NIF derived from NEAP-N and the BGS NIF as discussed in section 0 was used as a nitrate input. Point source discharges of nitrate, such as slurry heaps or septic tanks, have previously been estimated as contributing <1% of the total nitrate flux to groundwater (Sutton et al., 2011) and have not been considered in this study. It is possible that these could be important for some areas or aquifers. Transport of nitrate through the unsaturated zone on a 1 km scale was derived using the approach of Wang et al. (2012) (Wang et al, 2012a, 2012b). The total mass of nitrate in the unsaturated zone was derived for each year (1925 to 2050) for each 1 km grid cell and summed across the study area by aquifer. For any year, t, the total nitrate in unsaturated zone, NUSZ for a given grid cell with an unsaturated travel time, TTUSZ and a time-variant nitrate input function, NIF, can be calculated as:
Results and discussion
Figure 2.22 shows the temporal change in total nitrate stored in the unsaturated zone divided by aquifer. The total mass of nitrate in the unsaturated zone has increased substantially through time, peaking at approximately 1400 kiloton (kt) N in 2008. From 2008 onwards, the unsaturated zone becomes a net source of nitrate to groundwater and subsequently to surface water. The temporal change in nitrate storage in the unsaturated zone in 2015 is estimated to be approximately -5 kt N year-1. In 2015, the flux from the unsaturated zone to groundwater was approximately 72 kt 5 kt N year-1.
The increase in unsaturated zone nitrate storage is dominated by the Chalk, with 70% of the total mass in 2015. Increases are also observed in other aquifers such as the Permo-Triassic Sandstones (4% total mass in 2015), Oolitic Limestones (3% total mass) and numerous other locally important formations (23%). The Chalk, Permo-Trias and Oolites have peak mass years of 2015, 1991 and 1992 respectively. The year in which the total peak mass of unsaturated zone nitrate for England and Wales occurs is significantly affected by the majority of mass being in the Chalk.
The Chalk dominates the increase in unsaturated zone storage because of its large outcrop area, extensive agricultural land use (87%) and thick unsaturated zones (Wang et al., 2012b). Thick unsaturated zones results in long travel times and consequently a large increase in nitrate storage. Figure 2.22 shows the spatial distribution of nitrate stored in the unsaturated zone in 1960 and 2015. Increases in nitrate storage in the chalk of southern and north east England can be observed. Increases are particularly large in interfluve areas where travel times are long due to thick unsaturated zones.
The estimated peak nitrate mass of 1400 kt N is substantially greater than previous first approximations of 300 kt N (Worrall et al., 2009). However, in general this study corroborates with previous work suggesting that the subsurface is a significant store of reactive nitrogen. Whilst the total nitrate storage in the unsaturated zone is now decreasing, travel times in the saturated zone can be considerable (Wang et al., 2016) and consequently the peak saturated zone mass may not have occurred yet. Further research is required to assess how this storage compares with other postulated terrestrial stores such as in-stream N retention and terrestrial N uptake in land not in production.
The approach adopted in this analysis and that of Wang et al. (2012b) is likely to be beneficial for the targeting of catchment management activities at national and regional scales. For example, Figure 2.22 illustrates that legacy nitrate in the unsaturated zone at a national scale is dominated by the Chalk. Figure 2.23 shows that within the Chalk, there a substantial historical mass of nitrate in the unsaturated zone of southern England, particularly in interfluve areas where travel times are long. Consequently, environmental managers should take into account this mass when considering the implementation of catchment mitigation measures in attempts to improve groundwater and surface water quality. This could also be important when setting environmental objectives (such as for the WFD status assessment) which involve a simple assessment of water quality metrics, e.g. measured concentrations and associated statistics, and to demonstrate their achievement.
Assessment of NTB model improvements
Table 2.4 sets out the improvements to the model and their influence on modelled results. Of these, the use of a spatially variable NIF function is probably the most important. Figure 2.24 shows the impact on nitrate concentrations in groundwater from the use of the spatially variable (composite) NIF in the 28 zones of the process-based model. This highlights the overestimate for non- agricultural areas, particularly uplands in the original NIF.
The NIF function and the revised water levels were used in the process-based modelling but the nitrate velocities from 250k scale mapping were not incorporated at this stage.
|NIF function||Spatially and temporarily variable||1. Better represent spatially distributed land-uses
2. Realistic concentrations at the water table
|National scale||Previously, when using the single NIF, it was assumed that a single arable land-use covers aquifers across the England and Wales, thus leading to over-estimated nitrate concentrations in aquifers overlain by non-agricultural land|
|Water levels||Now uses OS river data||More realistic unsaturated zone thickness can be used to derive nitrate travel time in the USZs||National scale||Allowed estimation of more reasonable parameters such as transmissivity and aquifer thickness|
|Geological mapping scale||250k||Represents layered aquifers better||Regional and catchment scale||Important in Jurassic limestones and Carboniferous strata, such as the Coal Measures. Provides potential for process-based treatment depending on data availability|
|Processed based model||Quantifiable unsaturated zone travel time for 28 aquifer areas||1. The simplified nitrate dilution conceptual model provides a way to evaluate/calibrate the NTB model
2. Make the model applicable at regions where have limited information on the nitrate velocity in the USZs
3. Allows impact of future scenarios — e.g. climate change
|Regional and catchment scale||It provides ways to evaluate the modelled results and make it possible to use alternative available datasets to simulate the nitrate transport in the groundwater system. For example, it can use nitrate velocities in the USZs or derive these values using recharge and aquifer properties.|
|Represents low permeability superficial deposits||1. Nitrate routed between surface runoff and recharge to ground
2. This makes it possible to consider the nitrate transport time in low permeability drift
3. The nitrate dilution simulation in aquifers became more realistic when involving the modelling cells of aquifers overlain by low permeability drift
|Regional and catchment scale||The final improved model has been successfully calibrated. However, Figure 2.18 shows how the introduction of low permeability drifts impact the modelled results. Table 2.5 shows how parameters changed after considering the nitrate transport in low permeability drifts taking the Chalk of East Anglia as an example.|
|Unsaturated zone storage||National scale summary||Delineates bulk unsaturated storage of nitrate and temporal changes||All scales||Targeting of catchment management issues|
|Parameter (units)||Before considering the nitrate transport in low permeability drift||After considering the low permeability drift|
|Taquifer (m2 day-1)||1280||1470|
- WANG, L, BARKWITH, A, JACKSON, C, and ELLIS, M. 2012a. SLiM: an improved soil moisture balance method to simulate runoff and potential groundwater recharge processes using spatio-temporal weather and catchment characteristics. The 12th UK CARE Annual General Meeting. Bristol, UK.
- WANG, L, STUART, M E, BLOOMFIELD, J P, BUTCHER, A S, GOODDY, D C, MCKENZIE, A A, LEWIS, M A, and WILLIAMS, A T. 2012b. Prediction of the arrival of peak nitrate concentrations at the water table at the regional scale in Great Britain. Hydrological Processes, Vol. 26, 226–239.
- ANTHONY, S, QUINN, P, and LORD, E. 1996. Catchment scale modelling of nitrate leaching. Aspects of Applied Biology, Vol. 46, 23–32.
- ENVIRONMENT AGENCY. 2007. Making Information Available for Integrated Catchment Management. Science Report, SC060035/SR2007 (Bristol, UK).
- LORD, E I, and ANTHONY, S. 2000. MAGPIE: A modelling framework for evaluating nitrate losses at national and catchment scales. Soil Use and Management, Vol. 16, 167–174.
- SILGRAM, M, WARING, R, ANTHONY, S, and WEBB, J. 2001. Intercomparison of national & IPCC methods for estimating nitrogen loss from agricultural land. Nutrient Cycling in Agroecosystems, Vol. 60, 189–195.
- LORD, E I. 1992. Modelling of nitrate leaching: Nitrate Sensitive Areas. Aspects of Applied Biology, Vol. 30, 19–28.
- SCHOLEFIELD, D, LOCKYER, D R, WHITEHEAD, D C, and TYSON, K C. 1991. A model to predict transformations and losses of nitrogen in UK pastures grazed by beef cattle. Plant Soil, Vol. 132, 165–177.
- ALLEN, D J, BREWERTON, L J, COLEBY, L M, GIBBS, B R, LEWIS, M A, MACDONALD, A M, WAGSTAFF, S J, and WILLIAMS, A T. 1997. The physical properties of major aquifers in England and Wales. British Geological Survey Technical Report WD/97/34, Environment Agency R&D Publication 8. British Geological Survey Technical Report WD/97/34, Environment Agency R&D Publication 8.
- JONES, H K, MORRIS, B L, CHENEY, C S, BREWERTON, L J, MERRIN, P D, LEWIS, M A, MACDONALD, A M, COLEBY, L M, TALBOT, J C, MCKENZIE, A A, BIRD, M J, CUNNINGHAM, J, and ROBINSON, V K. 2000. The physical properties of minor aquifers in England and Wales. British Geological Survey Technical Report WD/00/4, Environment Agency R&D Publication 68.
- WANG, L, STUART, M E, LEWIS, M A, WARD, R S, SKIRVIN, D, NADEN, P S, and COLLINS, A L. 2016. The changing trend in nitrate concentrations in the major aquifers due to historical nitrate loading from agricultural land in England and Wales. Science of the Total Environment, Vol. 542, 694–705.
- HABERLANDT, U, KLÖCKING, B, KRYSANOVA, V, and BECKER, A. 2001. Regionalisation of the base flow index from dynamically simulated flow components—a case study in the Elbe River Basin. Journal of Hydrology, Vol. 248, 35–53.
- HAAN, M M, RUSSELL, J R, POWERS, W J, KOVAR, J L, and BENNING, J L. 2006. Grazing management effects on sediment and phosphorus in surface runoff. Rangeland ecology & management, Vol. 59, 607–615.
- LANGE, J, GREENBAUM, N, HUSARY, S, GHANEM, M, LEIBUNDGUT, C, and SCHICK, A P. 2003. Runoff generation from successive simulated rainfalls on a rocky, semi-arid, Mediterranean hillslope. Hydrological Processes, Vol. 17, 279–296.
- TANI, M. 1997. Runoff generation processes estimated from hydrological observations on a steep forested hillslope with a thin soil layer. Journal of Hydrology, Vol. 200, 84–109.
- TANI, M, and ABE, T. 1987. Analysis of stormflow and its source area expansion through a simple kinematic wave equation. Forest Hydrology and Watershed Management, Proceedings of the Vancouver Symposium. SWANSON, R H, BERNIER, P Y, and WOODARD, P D (editors). 167. (Wallingford: IAHS.)
- NASH, J E, and SUTCLIFFE, J V. 1970. River flow forecasting through conceptual models part I—A discussion of principles. Journal of Hydrology, Vol. 10, 282–290.
- BUTCHER, A, GRIFFITHS, K, LAPWORTH, D, HUMPAGE, A, and BURKE, S. 2009. Investigation of rising nitrate concentrations in groundwater in the Eden Valley, Cumbria. 4, estimating recharge rates through glacial till using an applied tracer technique. British Geological Survey Open Report OR/09/059.
- LEONARD, R A, and KNISEL, W G. 1988. Evaluating groundwater contamination potential from herbicide use. Weed Technology, 207–216.
- RAO, P, HORNSBY, A, and JESSUP, R. 1985. Indices for ranking the potential for pesticide contamination of groundwater. Proceedings of the Soil and Crop Science Society of Florida, Vol. 44, 1–8.
- RAO, P, and JESSUP, R. 1983. Sorption and movement of pesticides and other toxic organic substances in soils. 183–201 in Chemical Mobility and Reactivity in Soil Systems. NELSON, D W, TANJI, K K, and ELRICK, D E (editors). (Madison, WI: Am. Soc. Agron.)
- LEE, L J E, LAWRENCE, D S L, and PRICE, M. 2006. Analysis of water-level response to rainfall and implications for recharge pathways in the Chalk aquifer, SE England. Journal of Hydrology, Vol. 330, 604–620.
- MACKAY, J D, JACKSON, C R, and WANG, L. 2014. A lumped conceptual model to simulate groundwater level time-series. Environmental Modelling & Software, Vol. 61, 229–245.
- WANG, L. 2015. Estimating groundwater recharge for Great Britain considering the heterogeneous thickness of low permeability superficial drifts. The 17th IWA International Conference on Diffuse Pollution and Eutrophication, Berlin, Germany, 13–18 Sept 2015.
- WANG, L, ZHENG, R, LIN, W, MU, X, and HU, Y. 2014. Estimating groundwater recharge for Great Britain. The 14th UK CARE Annual General Meeting, UK Chinese Association of Resources and Environment, Bristol, UK, 29 Nov 2014.
- GRIFFITHS, K J, MACDONALD, A M, ROBINS, N S, MERRITT, J, BOOTH, S J, JOHNSON, D, and MCCONVEY, P J. 2011. Improving the characterisation of Quaternary deposits for groundwater vulnerability assessments using maps of recharge and attenuation potential. Quarterly Journal of Engineering Geology and Hydrogeology, Vol. 44, 49–61.
- ROSS, C A M, BATH, A H, ENTWISLE, D C, CAVE, M R, FRY, M B, GREEN, K A, and REEDER, S. 1989. Hydrochemistry of porewaters from glacial till and chalk at the Killingholme Site, South Humberside. British Geological Survey Technical Report WE/89/25.
- KLINCK, B A, BARKER, J A, NOY, D J, and WEALTHALL, G P. 1996. Mechanisms and rates of recharge through glacial till: Experimental and modelling studies from a Norfolk site. British Geological Survey Technical Report WE/96/1.
- MARKS, R J, LAWRENCE, A R, WHITEHEAD, E J, COBBING, J E, MANSOUR, M M, DARLING, W G, and HUGHES, A G. 2004. Chalk recharge beneatth thick Till deposits. British Geological Survey Internal Report, IR/04/179.
- CLOSE, M. 2010. Critical review of contaminant transport time through the vadose zone. (Christchurch: Environment Canterbury.) ISBN 1927137543
- ENVIRONMENT AGENCY. 2005. Attenuation of nitrate in the sub-surface environment. Science Report, -030155/SR2 (Bristol, UK).
- ENVIRONMENT AGENCY. 2007. Making Information Available for Integrated Catchment Management. Science Report, SC060035/SR2007 (Bristol, UK).
- BUTCHER, A S, LAWRENCE, A R, TRIBE, E, HANNAY, S, HASAN, K, and MERRIN, P D. 2005. Investigation of rising nitrate concentrations in groundwater in the Eden Valley, Cumbria. Summary Report 1: Catchment water quality survey. British Geological Survey Commissioned Report, CR/05/258N.
- KINNIBURGH, D G, GALE, I N, GOODDY, D C, DARLING, W G, MARKS, R J, GIBBS, B R, COLEBY, L M, BIRD, M J, and WEST, J M. 1999. Denitrification in the unsaturated zones of the British Chalk and Sherwood Sandstone aquifers. British Geological Survey Technical Report WD/99/2.
- PARRIS, K. 1998. Agricultural nutrient balances as agri-environmental indicators: an OECD perspective. Environmental Pollution, Vol. 102, 219–225.
- WORRALL, F, BURT, T, HOWDEN, N, and WHELAN, M. 2009. Fluvial flux of nitrogen from Great Britain 1974–2005 in the context of the terrestrial nitrogen budget of Great Britain. Global Biogeochemical Cycles, Vol. 23.
- SUTTON, M A, HOWARD, C M, ERISMAN, J W, BILLEN, G, BLEEKER, A, GRENNFELT, P, VAN GRINSVEN, H, and GRIZZETTI, B. 2011. The European nitrogen assessment: sources, effects and policy perspectives. (Cambridge University Press.)