OR/13/031 Material and methods

Tye, A M, Hurst, M D, and Barkwith, A K A P. 2013. Nene phosphate in sediment investigation — Environment Agency project ref: 30258. (Water Framework Directive). British Geological Survey Internal Report, OR/13/031.

The river Nene is the UK’s 10th longest river and rises close to the village of Badby near Daventry, flowing in a north-easterly direction out to the Wash via the town of Northampton and City of Peterborough (Figure 1). It is 161 km long and has a catchment area of 2270 km². The catchment for the upper six water bodies of interest in this study has an area of 1590 km². The river is navigable from the sea as far as Northampton where it connects with the Grand Union Canal. The floodplain is relatively wide (from a few hundred meters to a couple of km) and the channel frequently bifurcates and rejoins (Williams & Fawthrop, 1988; Meadows, 2007). A major characteristic of the catchment is that the river generally has very slow flow. The river course falls from ~160 m AOD at source to ~6 m AOD at Peterborough. The majority of this fall occurs in the first 9.5 km, with the channel lying at 80 m AOD at Weedon (Meadows, 2007). In the first section the river drops 1 m in 270 m and by the time it reaches Northampton it drops 1 m every 1500–2000 m and then by Thrapston it drops 1 m in every 3000 m (Meadows, 2007).

Figure 1    Map of the sampling locations in each of the six water bodies of the River Nene referred to in this study. The symbols for each water body are colour coded to represent the ‘Ecological status’ for that water body (Green = good; Orange = moderate; Red = poor).

Bedrock geology

The bedrock geology of the River Nene catchment comprises rocks of the Lias group, oolites and the Kellaways and Oxford clay formations (Figure 2). The Lias group rocks that the Nene flows over includes in order, the Charnmouth Mudstone Formation, the Dyrham Formation (interbedded siltstones and mudstones) and predominantly the Whitby Mudstone formation. Descriptions of the lithologies are from Cox et al. (1999).

Charnmouth Mudstone formation: Dark grey laminated shale, and dark to pale grey and bluish grey mudstones; occasional concretionary and tabular beds of argillaceous limestone; abundant argillaceous limestone, phosphatic or sideritic nodules in some areas; organic-rich ‘paper shales’ are found at some levels. Silty and fine sandy beds are particularly found in the upper part.

Dyrham Formation: Pale to dark grey and greenish grey, silty and sandy mudstone, with interbeds of silt or very fine sand (in some cases muddy or silty), weathering to a brown or yellow colour. There are impersistent beds or ferruginous limestone (some ooidal) and sandstone which may be very fine or laminated. There are occasional large argillaceous, silty or sandy limestone nodules.

Whitby Mudstone formation: Medium and dark grey fossiliferous mudstone and siltstone, laminated and bituminous in part, with thin siltstone or silty mudstone beds and rare fine-grained calcareous sandstone beds. Dense, smooth argillaceous limestone nodules are very common at some horizons and phosphatic nodules are found at some levels. Nodular and fossiliferous limestones occur at the base in some areas.

Figure 2    Bedrock geology map of the main six water bodies of the River Nene based on the BGS 625k scale geological map.

Superficial geology

The Nene Catchment was beyond the ice limits of the Devensian stadial but some glacial outwash is found in the catchment. These glaciogenic sediments are underlain by sands and gravels (Milton Sand) that represent earlier trunk rivers (Brown et al., 1994). Significant extraction of this sand and gravel has occurred along some reaches of the river, leaving a range of adjacent wetland environments. Major sand and gravel deposits are marked on Figure 1. On the modern floodplains alluvial soils and deposits are found.

Sampling strategy and site selection

The body of the river Nene sampled was from the source waters at Badby down to Stibbington, west of Peterborough. This section of the River Nene is non-tidal and is divided into six water bodies by the Environment Agency. The sampling strategy was designed to obtain robust estimates of sediment TP and OEP concentrations in each of the six water bodies with a reasonable estimate of variability, whilst keeping the overall analytical cost within the budget provided by the Environment Agency. Therefore, five cores were collected from each water body which were to be homogenised prior to analysis of TP and OEP, whilst a sixth core (Core_D) was extracted from each water body to allow analysis of TP and OEP variation with depth. This core was divided into 5 segments, with the top section always being 0–5 cm depth from the sediment surface. This is a standard depth on which phosphate dynamics analysis is undertaken (Jarvie et al. 2005). This 0–5 cm section was used to calculate ‘Effective phosphorous concentrations (EPC₀)’ (see Determining the effective phosphorus concentrations (EPC₀) in sediments) and the rate of adsorption and desorption (see Determining the rate of P uptake or desorption in sediments (Kinetics)) of phosphate. This was undertaken in conjunction with data collected from chemical analysis of a water sample for SRP taken in the water column above the sediment core. This water sample was filtered immediately after collection through a 0.45 µm syringe filter and kept in the dark until analysis, within 5 days of collection (Jarvie et al. 2002).

Figure 1 shows the location of the cores extracted from the six water bodies. Site selection was initially based on a desk survey of the six water bodies. Consideration was given to under water infrastructure through line searches. For water bodies 1 & 2 and several samples of water body 3, the cores were collected by wading in the river. For water bodies 4, 5 and 6 a boat was used to enable sampling. Sampling positions were partially determined for these cores by the available boat launch sites and travelling distance upstream and downstream.

Collection of sediments

Overview

Sampling was delayed from the original planned date (October/November 2012) because of the very wet autumn/winter of 2012/13. Consequently the EA placed high stream alerts on the river which prevented sampling. High stream alerts were briefly lifted in the first week of January 2013, allowing water bodies 1 & 2 to be sampled. Further rain meant that high stream alerts were next cancelled on the 14-02-13 allowing boat work to begin shortly after. The subsequent sampling of the sediments therefore needs to be put into the context of unusually long periods of flooding and high stream flows experienced during the autumn and winter period of 2012–2013.

Sediment collection

Where possible, cores were taken by pushing a length of polycarbonate tube (diameter = 58 mm internal diameter) into the sediment. Core retention was through the use of a core catcher, designed by BGS workshops, at the base of the tube. However, in the upper two water bodies there was little sediment so grab samples of the pebble and sand lag were taken. The total depth of sediment at each location was examined by probing at each site. This data was used to estimate the total quantities of TP and OEP within each water body. At the site in each water body where Core_D was extracted, additional measurements were taken. These included (i) the width of river, (ii) the speed of water movement, and (iii) the depth of water. These measurements were required in the calculation of SRP desorption and adsorption fluxes to/from the sediment.

Plate 1    Sampling sediment in water body 3.

Laboratory material and methods

Sample preparation

Five core samples taken from each water body were air-dried before being sieved to <2 mm. The mass of air-dry sediment from each of these cores in each water body was measured to estimate the average bulk density of the sediment in the cores (g cm^-3). Core_D from each water body was cut into 5 sections, the top 0–5 cm being wet-sieved <2 mm for isotherm and kinetic analysis. After isotherm and kinetic analysis was completed the remainder of this sample was air-dried for other analyses. The other four sections were air dried immediately before being sieved to <2 mm.

Sediment bulk density

The bulk density of sediment is required to estimate the mass of sediment from the volume of <2 mm sediment. This will allow calculation of TP and OEP in the sediment of each water body. Bulk Density (D_b) is calculated as:

Where Md = the total mass of oven dry soil, V = total volume, Ms = mass of stones and Vs = volume of stones.

Particle size analysis

Particle size distribution (PSD) of sediment samples (0.01 to 2000 mm) was determined using a Beckman Coulter LS13 320 laser diffraction particle size analyser. Prior to analysis organic matter was removed using H₂O₂.

Loss on ignition — organic matter content

Estimates of organic matter in the samples from the sediment cores were obtained by igniting 1g samples at 475°C for 4 hours. The difference in mass before and after ignition provides an estimate of the organic matter content.

TP measurements in sediments

TP analysis was undertaken on subsamples of the dried sediment (~30g) which were ground in an agate mill to <150 µm. Milled samples were digested by accurately weighing 0.25g of soil into a Savillex™ vial and adding HF, HNO₃ and HClO₄ concentrated and analytical grade acids, with a subsequent stepped heating program up to 170°C overnight, the purpose being the digestion of silicate and oxide phases. The dry residue was re-constituted after warming with MQ water, HNO₃ and H₂O₂, to 25 mL of 5% v/v HNO₃ and stored in HDPE bottles. Reference materials (NIST SRM2710, SRM2711, GSS-6, BGS102 and BCR-2), duplicated samples and blanks were all prepared in a similar manner to check accuracy of the analytical and digestion method. In addition to the quantification of TP, a range of elements associated with P mineralogy and phosphate sorption were measured including Al, Fe, Mn and Ca by ICP-MS in the BGS UKAS and MCERTS accredited laboratory.

OEP measurements in sediments

OEP concentrations were determined using the Olsen P method (Olsen et al., 1954) on the air-dry homogenised samples from each of the cores. Olsen P is 0.5M NaHCO₃ adjusted to pH 8.5. Samples were shaken in a 1:20 sediment: solution ratio for 30 minutes before centrifuging at 2500 rpm for 15 minutes. Analysis was carried out in triplicate and the results expressed on a dry weight basis.

Calculation of effective phosphorus concentration

The potential relationship between sorbed and solution SRP was examined with isotherms using sediment from the top 5 cm of Core_D from each water body. The isotherm was based on the Freundlich model (equation 2) which was fitted using the least squares method. The Freundlich model takes the form:

Where ΔN_a = change in adsorbed P (µmol g^-1), Kf is the Freundlich constant, Ci is the concentration of SRP in solution and n is a constant. Using the fitted isotherms the Equilibrium Phosphorus Concentrations (EPC₀) was calculated for each sample (see House et al., 1995; Jarvie et al., 2005). When considering the dynamics of the interactions of SRP with river sediment, it is important to have knowledge of whether the sediment has the potential to act as a source or a sink. The EPC₀ represents the equilibrium concentration of SRP in solution (i.e. when there is no net sorption or desorption over a 24 hour period). Thus, when SRP concentrations in the overlying water are greater than EPC₀, the sediment has the potential to sorb SRP from the water column and act as a sink. By contrast when SRP concentrations are less than the EPC₀, the sediment has the potential to release SRP to the water column and act as a source. EPC₀ analysis was undertaken on wet sediment (0–5 cm) within seven days of sampling to minimize sample deterioration as suggested by Jarvie et al. (2005). The methodology for the isotherms was similar to that applied by Palmer-Felgate et al. (2011). A measured mass of sediment was placed in 6 bottles with 200 ml of a synthetic water composition roughly matching the major element chemistry of the River Nene (2 mmol CaCl₂). The bottles were spiked with different concentrations of KH₂PO₄, and placed in an orbital shaker in the dark at 10°C for 24 hours. Samples were then centrifuged and their SRP concentration determined (Murphy and Riley, 1962). In addition, the linear portion of the isotherms was used to derive K_d values (L kg^-1) (Jarvie et al., 2005).

Calculation of rate constants for SRPsorption/desorption from sediments

To calculate the rate of release or sorption of SRP from/to sediments, the kinetic methodology described by Jarvie et al. (2005) was used. For SRP release experiments (where EPC₀ > dissolved SRP), a measured mass of wet sediment (equivalent to 0.5g dry sediment) was placed in polypropylene bottles with 200 ml of CaCl₂ solution, pre-chilled to 10°C, with no additions of KH₂PO₄. For SRP uptake experiments (where SRP > EPC₀), the synthetic river solution in each bottle was spiked with KH₂PO₄ to an appropriate SRP concentration (greater than the EPC₀ and close to the measured SRP concentration in the water column at the time of sediment sampling). Bottles were then placed in an orbital incubator in the dark and shaken at 150 rpm at 10°C and aliquots removed after specific time intervals (5 mins, 15 mins, 30 mins, 1h, 3h, 6h, 15h, 24h), then centrifuged and analysed for their SRP concentrations. The sorption/desorption rates of SRP (House and Warwick, 1999; Jarvie et al. 2005) to or from the sediments were calculated based on the equation:

R is the change in amount of orthophosphate sorbed (µmol g^-1 h^-1), C_t is the orthophosphate concentration (µmol l^-1) in the overlying water, C₀ is the orthophosphate concentration in the overlying water after 24 h (µmol l^-1), K_ris a rate constant (µmol^1-n lⁿ g^-1 h^-1) and n is a power term. The Nelder-Mead algorithm as implemented in Matlab (Lagarias et al., 1998) was used to determine the parameter values by minimizing the squared difference between the observed concentrations and those predicted by the rate curve. Using the calculated values the quantity of orthophosphate released from the sediment or removed from solution at any time, d_t, is obtained from equation 4 below.

Where M is the amount of orthophosphate sorbed (mmol), C_t is the concentration in solution at any time (t), EPC₀ is the equilibrium P concentration and S is the estimated mass of fine (<2 mm) sediment (g) in the reach of the river.

The fine bed-sediment mass S(g) within each waterbody was taken from calculation of sediment volumes obtained by probing and the mean bulk density of the sediment determined from the analyses of the cores. The flux of orthophosphate to/from the bed sediment is then calculated by integrating Eq. (2) over the residence time (T_res) of river water within the length of the river reach:

where D_cs is the mean cross-sectional depth of the reach, the length of the river reach (L_riv), and Q is the mean annual discharge of the reach and W_Cs is the average width of the cross section. Bed sediment orthophosphate fluxes to the boundary layer of the river channel are expressed as g P T_res^-1. The annual flux of orthophosphate to/from the bed sediment into the boundary layer can then be estimated by multiplying the flux for this residence time by the number of residence time periods in a calendar year. However, this ignores seasonal effects.

Molybdenum blue method for analysis of orthophosphate (PO₄-P)

Phosphate concentrations in (a) river water samples (b) Olsen P extracts, (c) isotherm analysis and (d) kinetic analysis were undertaken using the molybdenum blue method (Murphy & Riley, 1962). All samples used a calibration curve with concentrations of 0, 0.125. 0.25, 0.5, 0.75 and 1 mM P in 100 ml volumetric flasks; these being equivalent to ~3.7, 7.5, 15, 22.5 and 30 µg P. Calibration curves and analysis were undertaken using a Perkin Elmer Lambada 35 uv/vis spectrometer.

Figure 3    Photo showing (a) calibration and (b) Olsen P extract solutions.

A combined calibration curve from all the phosphate analysis runs is shown in Figure 4. Concentrations were converted from µM L^-1 to µg L^-1 to generate the relationship in equation 6 from values shown in figure 4 which relates PO₄-P concentrations to absorbance readings:

This combined calibration equation was used to determine all P-PO₄ concentrations in solutions from absorbance readings. In addition, ‘Limits of Detection (LOD)’ were calculated for each of the different matrices used. These included water, 2 mMol CaCl₂, and the 0.5 M NaHCO₃ Olsen P extractions (Table 1). Limits of detection were calculated by making up 10 samples of each of the blank matrices, measuring the absorbance and calculating the standard deviation. The LOD was then calculated as the standard deviation multiplied by 3.

Figure 4    Calibration curve for molybdenum blue method for determining PO₄-P (µm L^-1). Standard deviations for each point are shown but are generally contained within the symbol.

Table 1    Calculated Limits of detection for each of the matrices used in PO₄-P analysis

Sample Matrix	Mean of blank (n=10) (µg P)	Standard deviation	LOD µg P (3 x SD)
Water	0.08	0.07	0.22
2 mM 20 ml	0.11	0.11	0.34
2 mM 40 ml	0.22	0.07	0.22
2 mM 60 ml	0.72	0.06	0.18
2 mM 80 ml	0.74	0.12	0.37
Olsen P	1.61	0.16	0.59

Calculating sediment fluxes using a landscape evolution model

Background

The modelling phase of this project utilises a complex, distributed landscape evolution platform to derive sediment transport, erosion and deposition over the last century. During this period, the widespread heavy usage of phosphate fertiliser was introduced into the arable farming industry as a means of increasing production. Through processes of runoff and fluvial sediment deposition, phosphates sorbed to sediment during this period can accumulate in the river bed. Phosphates stored in this manner can later be remobilised and transported during a storm event. The CDP (CAESAR-DESC Platform), based on the well established CAESAR landscape evolution model, uses distributed rainfall, potential evaporation and soil properties to drive an integrated surface-subsurface hydrological model. Sediment transport in the model is dictated by the regional hydrology and terrain characteristics. Continuous monitoring of erosion and deposition allows the spatial distribution of sediments to be mapped.

Model description

Introduction

This CDP was created by combining algorithms from various backgrounds to produce a modelling system from which a variety of earth system interactions can be explored. A modified version of the well established Cellular Automaton Evolutionary Slope and River (CAESAR) model (Coulthard and Van De Wiel, 2006), CAESAR-Listflood (Bates et al., 2010) is used as the platform kernel. CAESAR is quasi three-dimensional and incorporates a modular design, with great versatility in the range of simulated spatio-temporal scales to which it can be applied. CAESAR has been used to investigate a variety of sediment transport, erosional and depositional processes under differing climatic and land use pressures in river reaches and catchments around the world. Improvement of the surface hydrological representation prior to this study allowed bespoke distributed fluvial soil pathway and confined groundwater models to be added to the CAESAR model, forming the CDP. The recently updated components have previously been validated against analytical solutions and observed data, and the coupled model has been applied at a regional scale in the Eden Valley catchment, Cumbria (UK), where the role that subsurface water fluxes play in shaping catchment scale geomorphology was examined (Barkwith et al., 2012).

Water partitioning

The partitioning of rainfall between, evaporation, runoff and recharge to groundwater in the CDP is achieved using a soil water balance model (Wang et al., 2012). We implement a simple technique that represents potential groundwater recharge and runoff processes based on spatio- temporally distributed soil moisture conditions. Soil moisture is a function of; rainfall, potential evapotranspiration, soil moisture condition, topography, soil types, crop type and base flow index (BFI). The method we use represents these soil water processes responding to variable soil water storage properties (see Rushton, 2003) and vegetation growth stages (see Allen et al., 1998). When soil moisture reaches field capacity and is unable to store further additions of water, water drains freely in the saturated soil. Additional water inputs to the soil result in lateral runoff (routed by Lisflood), if a gradient exists towards adjacent locations, and if required as percolation downwards through the saturated soil (groundwater recharge). If the rainfall intensity is greater than the capacity of soil to maintain infiltration, water accumulates on the soil surface and becomes surface runoff. Water not accounted for in soil storage, evapotranspiration or uptake by vegetation is termed excess water. Excess water is divided between runoff and recharge to groundwater based on a base flow index parameter. Base flow index is an average surface to subsurface water partitioning ratio reflecting the permeable nature of the catchment in addition to other catchment characteristics. The base flow index parameter is estimated, by performing a baseflow separation on a river flow time-series, and further refined through an iterative calibration process. In general, greater runoff and reduced recharge is observed in areas with steeper slopes. Consequently, average and nodal terrain gradient are factored into the calculation of recharge and runoff.

Surface water routing

Routing of surface runoff and channel flow is controlled by Lisflood (Van der Knijff et al., 2010), a two-dimensional hydrodynamic flow model. Routing is based on the one dimensional Saint-Venant equations, as modified by Bates et al (2010), where surface flow takes into account; acceleration, water surface gradient and friction properties and has enhanced stability due to increased friction forcing water flow towards zero (Liang et al., 2006). If the flow depth, as calculated using the surface partitioning component, is above a defined threshold and stability criteria are met, the flux between cells at the proceeding time-step depends on the water surface gradient between the adjacent cells, Manning’s coefficient, the time-step length and the gravitational constant. For each time-step flow between neighbouring cells is computed and flow depths updated simultaneously at all points within the model domain. A full description of Lisflood is provided by Van der Knijff et al. (2010).

Groundwater flow

Groundwater flow is simulated using a two dimensional lattice of square cells interacting according to the von Neumann type neighbourhood in a technique similar to that utilised by Rothman (1988). Cells consist of distributed hydraulic conductivity and specific yield, and a datum referenced aquifer head that is modified through time. Aquifer head is updated at each time-step using Darcy’s law to calculate water flux between cells. Transmissivity can be approximated by multiplying hydraulic conductivity by aquifer depth. The total flux for a central cell is a combination of; fluxes to neighbouring cells and additional source and sink terms. Recharge from the surface hydrological component provides the source term and baseflow to the surface acts as a sink. Total fluxes are used to update aquifer heads at each point in the domain simultaneously at each time step using a discrete mass balance equation (see Ravazanni, 2011). Two user-defined lateral boundary condition types have been implemented into the CDP code, allowing a variety of hydrological situations to be represented. Specified (Dirichlet) boundary conditions fix aquifer head at the boundary to the initial condition and a no-flow (Neumann B) condition sets flux across the boundary to zero. The base of the aquifer is defined as impermeable; however, leakage to and from the base of the modelled aquifer is included in the flux algorithm as a secondary source and sink terms. The surface boundary allows water flux to be returned to the surface component as baseflow where aquifer head is greater than terrain height.

Sediment transport

Following the derivation of flow depth at each node, fluvial erosion and deposition are calculated, where fluvial transport is determined using the Wilcock and Crowe (2003) method, to calculate mixed-size particle movement. The Wilcock and Crowe formulation (equation 7) derives the sediment flux of a grain size fraction using the fractional volume of the particular sediment size (F_i), shear velocity (u_*), ratio of sediment to water density (r_ws), gravity (g_i), and the function W_i; which allows the total transport to be calculated from a shear stress factor (Φ) derived from shear stress ( ) of the fractional grain size and t_ri, a reference shear stress:

The Wilcock and Crowe method can transport sediment within a catchment as a suspended load or a bedload (Equation 10). Both depend on the volume of sediment (q_sed) transported per time- step (d_t), where the sediment transported from a central cell to a neighbouring cell (k) is calculated by equating the coefficient X to either; slope, for bedload transport, or velocity for suspended load transport. The calculation of suspended load transport is simplified, as it assumes an equal distribution of sediment throughout the water column. Sediment deposition also differs between the two transport types, with bedload sediment deposited (and subsequently re-entrained) at every time-step and suspended sediment deposition based on fall velocities.

Input parameters

Initialisation of the CDP requires a number of essential input items in the form of either gridded or list based ASCII files. Spatially discretised initial model inputs are entered as Cartesian grids, with header information containing cell size, domain dimensions and geo-referencing details.

Distributed daily rainfall, land use, DEM, soil type hydrology and evapotranspiration datasets are used to initialise and force the surface hydrology model. Field capacity, wilting point, maximum root depth and crop coefficient are used in conjunction with the reference evapotranspiration to calculate the PE values for different crop types. As the groundwater module is enabled for this study, the distributed hydraulic conductivity and specific yield need to be defined, initial groundwater levels need to be derived and the boundary conditions specified. A description of this process is contained in the ‘Setup’ section of this report. Platform sediment transport is based upon relative density of submerged sediment, grain size, flow depth and the water surface gradient. The latter two are derived from the platform hydrological model, while the distribution of grain size and density is determined through analysis of the cored sediment samples collected during the fieldwork phase of this study.

Output

Simulated changes to sediment volume, representing the past century, are presented for the Nene catchment and the six sub-catchment water bodies as defined by the Environment Agency in the project scope. Through analysis of the changes in volume during the simulation over the latter half of the century, and through a projection of these rates back in time, the time averaged and total twentieth century sediment flux rates and phosphate levels can be derived. Using modelled suspended sediment and bed load volumes calculated by the simulation at the exit of each water body of the Nene, the annual phosphorus and phosphate volumes released into the river can be estimated.

Model application

Conceptual model and study area

The overall Nene catchment area is approximately 2270 km². Despite the presence of some large population centres (Northampton, Wellingborough, Kettering, Corby and Peterborough) the catchment is largely rural. Areas surrounding the Fens and to the East of Peterborough, where arable crop is produced, are intensively farmed. To the west of Peterborough, the land is farmed less intensively.

The River Nene and its tributaries, the Kislingbury Branch, the Brampton Branch and Wootton Brook, meet close to Northampton and flow across gently undulating rural country to the flat plains around Peterborough, before entering the embanked tidal reach across the Fens. Much of the Fens lies below sea level relying on pumping stations for drainage.

The underlying geology of the catchment broadly comprises mudstones to the west of Northampton, limestones between Northampton and Peterborough, and clays to the east of Peterborough. Where the underlying rock is composed of non-porous mudstones, there are higher rates of rainfall runoff, which flows directly into the surrounding watercourses. In the areas where limestone or sandstone bedrock is present, runoff may infiltrate the rock, lagging the response of rivers to rainfall and reducing peak flood flows. Excess storm water contained in the groundwater store can make these areas prone to flooding with little subsequent rainfall. In the lower fenland areas downstream of Peterborough the predominance of peat soils and the low gradients means that water moves slowly to the river channels.

For this study we are not interested in the area to the east of Peterborough, and therefore, to reduce the computational expense, this region of the catchment is discarded from the numerical modelling.

Discretisation

Model boundaries and internal sub-catchment boundaries are defined in Figure 5. The groundwater and surface water catchments are bounded by the Neumann B condition, where no water can flow across the boundary. The exception to this setting is along the eastern boundary, where surface water can be discharged from the catchment in the form of a river. The Nene catchment is discretised into 200 m square grid cells, with a variable sub-daily time-step utilised for calculations of surface flow and sediment transport. The sub-daily time-step limits the amount of sediment that can be passed between cells, maintaining stability within the model during simulation. The selected grid spacing allows a viable execution time, whilst retaining enough detail to adequately represent the sediment transport system. In addition to a catchment wide assessment, the grid spacing also facilitates the sediment and water flux properties of the six sub-catchments to be assessed post-simulation.

Figure 5     Digital Elevation Model of the area used for landscape evolution modelling. The model boundaries are no flux type boundaries except for the eastern edge of the model domain. Water body outlets and catchment areas are also shown.

Parameterisation

Distributed daily rainfall, landuse (LCM2000), DEM, the hydrology of soil types (HOST), river flow and evapotranspiration (MORECS) datasets are used to initialise and force the surface hydrology model. The HOST dataset contains 29 soil classes of similar hydrological response and the BFI value for each class. Field capacity and wilting point of each soil type, and the maximum root depth and crop coefficient (which is multiplied by the reference evapotranspiration to calculate the PE values for different crop types) for each landuse type which was gathered from the experimental data of Allen et al. (1998) and the hydrological modelling work of Griffiths et al. (2008) and Young et al. (2008). The sediment transport components of the platform utilise the same DEM as the hydrology components, but also require bedrock elevation referenced to the same datum to determine the thickness of the sediment store layer. Grain size distributions as determined by field measurement of fluvial sediment are used to initialise the model.

Setup

Although it would be preferable to run the model for the entire century, access to rainfall and potential evapotranspiration data for the region is limited to the forty years 1962–2001. Representation of sediment flux rates for each of the water bodies over the periods not covered by the simulation will be projected back to 1910 and forward to 2010. The input data required by the model, and a short description of each, is contained in Table 2.

Table 2    Input datasets for the calibration, ‘spin-up’ and transient simulations. The ‘spin-up’ is the period where the model is allowed to reach a ‘steady state’ and accounts for initial conditions that are not correct (i.e. no water in the system and heterogeneous sediment distribution)

Input	Type	Description
Hydraulic Conductivity	Single	Value determined through calibration of groundwater levels and base flow return to river
Specific Yield	Single	Value determined through calibration of groundwater levels and base flow return to river
Potential Evapotranspiration	Distributed	1962–2002 at a 40 km daily resolution, as determined by the Met. Office MORECS dataset
Precipitation	Distributed	1962–2002 at a 25 km daily resolution, as determined by CEH dataset
Landuse	Distributed	LCM2000 land cover map, as determined by CEH dataset based on satellite data
HOST	Distributed	Fixed, 1 km national soil map dataset, integrating river catchment behaviour
Initial NSSS	Distributed	Determined through spin‐up to steady‐state value
Initial SMD	Distributed	Determined through spin‐up to steady‐state value
Initial Sediment Distr.	Heterogeneous	Nine grain size fractions determined by averaging field observations
Hydrological Boundary	Distributed	No‐flow at all boundaries, except at the river outlet
Debris Flow Method	Single	Sand‐pile method, with critical failure angle set to 30°

The percentage distribution by weight of nine grain size fractions was incorporated into the model domain (Table 3), which were determined by dry-sieving stream sediment from two core samples, taken in Water Bodies 1 and 2 respectively, and taking the mean of the two. This distribution is assumed to be fully mixed and applied throughout the catchment in the initial model setup. As the model evolves through time, sediment grain sizes are spatially redistributed across the active parts of the catchment. This distribution initially occurs as a strongly transient volumetric sediment flux that recedes as grain size distributions adjust to the local hydrological conditions. During the initial period, sediment transport fluxes will be much greater than the period that they represent. By negating sediment flux rates during this ‘spin-up’ period, which approaches a steady-state after approximately ten simulated years for the Nene catchment at the desired resolution, uncertainty in the modelled output is reduced.

The main simulation will comprise a forty year transient run representing 1962–2002, where sediment erosion and deposition are monitored on a 10-day basis and suspended sediment and river flows are monitored continuously.

Table 3     Catchment average initial grain size distribution as determined by field observation

Grain Size (mm)	<0.06	<0.16	<0.36	<0.55	<0.8	<1.5	<3.0	<7.0	<30.0
Percentage (%)	5.4	11.8	45.6	50.1	62.1	78.1	88.4	98.0	100.0

Modelling caveats

As with all numerical modelling techniques, there are several caveats associated with the datasets utilised by the CDP. Firstly, discretisation of the grid reduces all observed sub-grid variation to a nodal average. As we are working at the catchment scale, this should have little influence on the model output. The spatial resolution of the meteorological and soil based datasets is much greater than that of the grid, reducing the accuracy of the modelling. It is currently not possible to obtain higher resolution datasets without extensive field campaigns, which are beyond the scope of this study. Again, as the resolution of these datasets is an order of magnitude less than the catchment size, they do not impart much error in the final model output. The ability of the model to represent gauged river flow record, suggests that, although improvements are possible, the use of these datasets is justified. The homogenisation of sediment distribution across the catchment during initialisation is rectified by the model during attenuation to steady state during simulation. During this process, sediment is sorted in a distributed manner across the catchment, with upper reaches consisting of coarser materials and the lower reaches having a greater percentage of fines.

Meteorological and gauged river flow data

The distributed average 1962–2001 precipitation (mm d^-1) is presented in Figure 6. Due to the interaction of orography and the atmosphere boundary layer under the prevalent wind flow direction, precipitation is up to 33% greater in the elevated sections of the catchment in comparison to the lowest sections.

The Upton (west) and Orton (east) river gauging stations are used to calibrate the model hydrology between 1970 and 1979 (Figure 7). Although the gauging stations are at opposite ends of the catchment and have a flow values at differing orders of magnitude, they exhibit the same pattern of flow fluctuation. During the 70 years of available flow data (1939–2010), the highest flows were found on the 18th May 1947, remaining well above average for over two weeks. This peak coincided with widespread flooding associated with the rapid thaw of deep snow, which had accumulated during the previous winter. These river floods were the largest floods for over 200 years (RMS, 2007) although other major floods have occurred in 1976/77 and 1998. The peak flows during the 1947 flood are more than double any subsequent observations. The effect of this level of flooding on sediment fluxes for the Nene catchment is unknown and, therefore, there is an elevated level of uncertainty associated with projecting results to beyond this date.

Figure 6    Distributed mean daily precipitation in the Nene Catchment, averaged over the period 1962–2001.

Figure 7     Gauged river flow data for the Orton (Blue) and Upton (red) stations. Note the differing flow scales on the vertical axis.

For the calibration, a daily cycle of model output was specified for the calculation of groundwater level (mAOD), base flow (m³ day^-1), recharge (m³ day^-1) and surface water (m³ day^-1). To assess the errors associated with differing water flow routes, and to simplify the process, time series of simulated groundwater and surface flow were calibrated to observed data separately. The first phase of calibration removed the simulated base flow component from surface flows. These where compared to a base flow separated flow records from the Orton and Upton gauging stations for the surface component analysis. Once an acceptable agreement was achieved, the second phase of calibration was implemented, where the base flow component was re-instated in the simulation and a comparison to the non-separated gauged river flow dataset used for analysis. The same river gauging stations were used for both phases of the calibration.

Following the first phase of calibration, there is an acceptable match between simulated and observed river flows with the base flow component removed (see Figure 8 for the comparison at the Orton gauging station). The timing of large flow discharge events shows excellent agreement; however, the representation of peak discharge values is not as good. The disparity most likely arises from the misrepresentation of true rainfall and surface frictional characteristics under the current discretisation scheme, and would be improved with input data at a higher spatial resolution. Although the peak discharges are not fully represented in the simulation, the average peak discharge value is very similar to the observed value. Following the second phase of calibration, where base flow is introduced, the match between simulated and observed river flows is in some places poor (Figure 9). The worst representation occurs during recession of the river, where groundwater recession rates are not fully representative of the system. This suggests that the parameterisation of the groundwater model is not accurate. However, the simulated peak flow values show close agreement with the observed values. For this study the latter is the most important attribute of the calibrated output, as it is during these periods where the majority of sediment is transported.

Figure 8    Base flow separated surface runoff derived from the Orton gauged river flow dataset (blue) compared to the simulated runoff (red) from 1970–1979. These datasets were used to calibrate the surface flow characteristics of the CDP.

Figure 9    Surface flow (m³ s^-1) from the Orton gauged river flow dataset (blue) compared to the simulated runoff (red) from 1970–1979. These datasets were used to calibrate the groundwater surface flow characteristics of the CDP.

Simulation projection

The projection of results forwards and backwards in time through a post-processing analysis forms an important part of this project. It is necessary as we do not have meteorological data required to drive the hydrological model available for this catchment prior to 1962 or after 2002. Therefore, based on a statistical analysis of the simulated sediment flux rates, we projected the rates forwards and backwards in time.

Calculations of total phosphorus and phosphate in river sediments and waters

Total mass of P in sediment of the six water bodies

For each of the six water bodies we calculated the mean TP content from the five cores (VPC):

Where TP_i and BD_i are respectively, the TP concentration and bulk density of the ith core. The total mass or stock of P (Sed P_tot) in the sediment of each water body can then be computed by:

where S_vol is the estimate of the total volume of sediment in the channel (from work package 3). The total stock of P in the sediment of all six water bodies is the sum of the stock in each water body.

Extractable P stock in sediments of the six water bodies

The approach to estimate OEP stock in sediment will be the same as for the TP shown above, but in this case the data for OEP will be used instead of the TP concentrations in Equation 11. So, here mean volumetric OEP (VêPC) will be estimated using OEP P concentrations (TeP):

and the total mass or stock of OEP (Sed eP_tot) in the sediment of each water body can then be computed by:

The total stock of OEP in the sediment of all six water bodies is the sum of the mass in each water body.

Estimating the annual mass of orthophosphate adsorbed/desorbed to each water body and across all water bodies

The river system can be considered as a series of inter-connected reaches, each of which has a known mass of fine (<2 mm) sediment in the river bed which interacts with the overlying water column to sorb or release orthophosphate. The approach used is the same as that of Jarvie et al. (2005) where they assumed that the boundary layer of water that interacts with the bed sediment is 0.1 m thick, and orthophosphate is exchanged between the river bed and the volume of water within the boundary layer during the residence time of the river water within the reach. This exchange is characterised by: (i) the measured initial SRP concentration in the water column, (ii) the EPC₀ of the bed sediment and (iii) kinetic parameters for P-sorption, which were measured in the laboratory (see Calculation of rate constants for SRPsorption/desorption from sediments). Therefore orthophosphate flux estimates to/from bed sediments can be calculated for each of the 6 water bodies and for the Nene as a whole using the outputs of mean annual river flux from the CDP model which is calibrated to data from flow gauging stations. These fluxes provide a quantitative estimate of the orthophosphate uptake/release potential of the river bed sediments to the river water. This analysis will not provide predictions of river water orthophosphate concentrations, but it will provide an assessment of the magnitude of potential release or sorption of orthophosphate to water in each water body.