Scholarly article on topic 'Phosphate stable oxygen isotope variability within a temperate agricultural soil'

Phosphate stable oxygen isotope variability within a temperate agricultural soil Academic research paper on "History and archaeology"

Share paper
Academic journal
OECD Field of science
{Phosphorus / Grassland / "Spatial analysis" / "GW model" / "North Wyke farm platform"}

Abstract of research paper on History and archaeology, author of scientific article — Steven J. Granger, Paul Harris, Sabine Peukert, Rongrong Guo, Federica Tamburini, et al.

Abstract In this study, we conduct a spatial analysis of soil total phosphorus (TP), acid extractable phosphate (PO4) and the stable oxygen (O) isotope ratio within the PO4 molecule (δ18OPO4 ) from an intensively managed agricultural grassland site. Total P in the soil was found to range from 736 to 1952mgPkg−1, of which between 12 and 48% was extractable using a 1M HCl (HClPO4 ) solution with the two variables exhibiting a strong positive correlation. The δ18OPO4 of the extracted PO4 ranged from 17.0 to 21.6‰ with a mean of 18.8‰ (±0.8). While the spatial variability of Total P has been researched at various scales, this is the first study to assess the variability of soil δ18OPO4 at a field-scale resolution. We investigate whether or not δ18OPO4 variability has any significant relationship with: (i) itself with respect to spatial autocorrelation effects; and (ii) HClPO4 , elevation and slope - both globally and locally. Results indicate that δ18OPO4 was not spatially autocorrelated; and that δ18OPO4 was only weakly related to HClPO4 , elevation and slope, when considering the study field as a whole. Interestingly, the latter relationships appear to vary in strength locally. In particular, the δ18OPO4 to HClPO4 relationship may depend on the underlying soil class and/or on different field managements that had operated across an historical north-south field division of the study field, a division that had been removed four years prior to this study.

Academic research paper on topic "Phosphate stable oxygen isotope variability within a temperate agricultural soil"


Contents lists available at ScienceDirect


journal homepage:

Phosphate stable oxygen isotope variability within a temperate agricultural soil

Steven J. Granger a* Paul Harris a, Sabine Peukert a,\ Rongrong Guo b, Federica Tamburinic, Martin S.A. Blackwella, Nicholas J.K. Howden b, Steve McGrath d

a Rothamsted Research, North Wyke, Okehampton, Devon EX20 2SB, United Kingdom

b Queen's School of Engineering, University of Bristol, Senate House, Tyndall Avenue, Bristol BS8 ITH, United Kingdom c Institute of Agricultural Sciences, ETH Zurich, Research Station Eschikon 33, 8315 Lindau, Switzerland d Rothamsted Research, West Common, Harpenden, Hertfordshire, AL5 2JQ United Kingdom


In this study, we conduct a spatial analysis of soil total phosphorus (TP), acid extractable phosphate (PO4) and the stable oxygen (O) isotope ratio within the PO4 molecule (S18OPO4) from an intensively managed agricultural grassland site. Total P in the soil was found to range from 736 to 1952 mg P kg-1, of which between 12 and 48% was extractable using a 1 M HCl (HClPO4) solution with the two variables exhibiting a strong positive correlation. The 818OPo4 of the extracted PO4 ranged from 17.0 to 21.6% with a mean of 18.8% (±0.8). While the spatial variability of Total P has been researched at various scales, this is the first study to assess the variability of soil S18OPO4 at a field-scale resolution. We investigate whether or not S18OPO4 variability has any significant relationship with: (i) itself with respect to spatial autocorrelation effects; and (ii) HClPO4, elevation and slope - both globally and locally. Results indicate that S18OPO4 was not spatially autocorrelated; and that S18OPO4 was only weakly related to HClPO4, elevation and slope, when considering the study field as a whole. Interestingly, the latter relationships appear to vary in strength locally. In particular, the S18OPO4 to HClPO4 relationship may depend on the underlying soil class and/or on different field managements that had operated across an historical north-south field division of the study field, a division that had been removed four years prior to this study.

© 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://


Article history: Received 23 July 2015

Received in revised form 13 September 2016 Accepted 21 September 2016 Available online xxxx

Keywords: Phosphorus Grassland Spatial analysis GW model

North Wyke farm platform

1. Introduction

Phosphorus (P) is an essential element for plant growth and is applied to agricultural systems, often in large quantities, to underpin intensive levels of agricultural production (Haygarth et al., 2014). It can be applied as either inorganic mineral fertilizers, or via the spreading of animal wastes or other organic materials. Such wastes can occur either directly voided by the animal in the field or applied in bulk after storage while animals are housed. However, P can have significant detrimental effects when they move from land into water bodies. These effects range from direct toxicity (Lewis and Morris, 1986) through to indirect consequences such as eutrophication (Smith et al., 1999). Surface waters are particularly sensitive to P because critical concentrations of only a few tens of |ag of phosphate (PO4) can cause eutrophication, but are an order of magnitude lower than soil PO4 concentrations required for crop growth (Heathwaite and Dils, 2000). Identifying the different pollutant sources that are impacting on a water body is critical to

* Corresponding author.

E-mail address: (S.J. Granger).

1 Present address: FWAG SouthWest, Environment Department, County Hall, Taunton, Somerset, TA1 4DY, United Kingdom.

understand its ecosystem health. However apportioning a pollutant to any given source or sources is fraught with difficulty and in recent years techniques have been developed to try to elucidate pollutant origins (e.g. Baker et al., 2002, Collins et al., 1997; Old et al., 2012) and these include the use of natural abundance stable isotope ratios (e.g. Amberger et al., 1987; Granger et al., 2008). More recently still the stable isotope approach has been applied to P and although P only has one stable isotope, the technique uses the stable oxygen (O) isotope ratio within the PO4 molecule (S18OPOJ to isotopically characterise PO4 sources and transformations. However, data on the ô18OPO4 of different PO4 sources remains limited (Tamburini et al., 2014; Young et al., 2009).

One potential source of PO4 in water is from soil, which often receives P inputs in excess of requirements resulting in P accumulation within the soil (Haygarth et al., 1998a). Therefore, as with many other soil properties, an understanding of soil P variability is essential for designing sampling strategies or the evaluation of the effectiveness of diffuse water pollution mitigation measures (Goovaerts, 1998; Rivero et al., 2007). Despite this, few studies describe spatial variability of soil properties and their inter-relationships at a landscape scale (i.e. Marriott et al., 1997; Page et al., 2005). The spatial variability of a given soil property may be related to the combined action of several physical, chemical or biological processes that act at different spatial

0016-7061/© 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (

scales depending on the soil property and process of interest (Goovaerts, 1998). Soil spatial variability is present in both natural and agricultural systems, even if the latter have a long-term uniform management history (Goovaerts, 1998; Marriott et al., 1997). Understanding soil spatial variability is essential to land-based experiments at all scales and its omission is detrimental to the conclusions drawn from such experimental data.

In this study, we conduct a spatial analysis of soil 618OPO4 from an intensively managed agricultural grassland site. While the spatial variability of total P has been described at different scales by other researchers (e.g. Glendell et al., 2014; Page et al., 2005; Penn et al., 2007), this is the first study to assess the variability of soil 618OPO4 at a field-scale. We investigate whether or not 618OPO4 variability has any significant relationship with: (i) itself with respect to spatial autocorrelation effects; and (ii) soil P (extracted using 1 M HCl (HClPOJ), elevation and slope -both globally and locally within the study field. In particular, the following four null hypotheses are tested:

A. 618OPo4 does not strongly co-vary with HClPO4, with elevation and with slope across the field as a whole.

B. 618OPO4 is not a spatially-autocorrelated process.

C. The relationships of (A) do not significantly change in different areas of the field.

D. The relationship between 618OPO4 and HClPO4 is not conditional on the the under-lying soil class or different management histories.

Furthermore, in order to efficiently test the given set of null hypotheses, all statistical analyses are conducted in a manner that accounts for a certain sub-optimality in the study sample design. In particular, statistical methods are chosen to cater for significant areas of under- and over-sampling.

2. Methods

2.1. Study site

To characterise soil P spatial variability, a series of soil samples were

collected from one field of the Rothamsted Research 'North Wyke Farm

Platform', in south-west England (50.8°N, 3.0°W). The field sampled, re-

ferred to as 'Great Field', was located on a north-west facing hillslope

and comprises clay loam soil overlying the shales and thin subsidiary

sandstone bands of the Crackington formation (Harrod and Hogan,

2008). The soils can be further divided in three main types; Hallsworth (USDA Aaerichaplaquept, FAO Stagni-verticcambisol), Halstow (USDA

typichaplaquepts, FAO dysticgleysol), and Denbigh (USDA Dysticeutrochrept, FAO Stagni-eutriccambisol) (Harrod and Hogan,

2008). The long-term annual temperature and rainfall are 9.6 °C and 1056 mm, respectively, with a high proportion of rainfall occurring be-

tween October and March resulting in waterlogged soils.

Records of the historic farm management show that prior to 2010

the field comprised two separate areas (north 1.5 ha and south 5.6 ha)

with contrasting management histories. The northern part had been managed as permanent grassland for at least 25 years, whereas the southern part has been ploughed three times in the last 25 years, most recently in September 2007 when it was re-seeded with a ryegrass/clo-ver mixture following a previous winter barley crop. At the time of the study, both areas of the field had the same vegetation cover being classified under the National Vegetation Classification (NVC) category: "MG7 Lolium perenne Poa trivialis and related grasslands". The southern region, however, had a higher clover content (Trifolium repens), less

dense vegetation cover and approximately double the sward height

compared to the northern part. These management histories are representative of normal management cycles of intensive grassland management (Bilotta et al., 2007).

22. Sample design, collection and analysis

To quantify spatial variability in 618OPC,4, a sampling pattern was chosen with the view of a geostatistical analysis not only to 618OPO4, but to other soil variables, as presented in Peukert et al. (2012). In order to assess spatial variability across three different spatial scales, 78 soil samples were taken in total with samples at: (i) a broad scale (75 x 75 m grid); (ii) an intermediate scale (25 x 25 m grid); and (iii) a small scale (10 x 10 m grid). A hand-held GPS (Nomad Trimble, Sunnyvale, USA) was used to map and mark the sampling points.

All samples were collected in May 2011 to a soil depth of 7.5 cm and were oven dried at 105 °C for 24 h. Dried soils were then sieved through a 2 mm mesh. Total P (TP) was determined at an external laboratory (Natural Resource Management, Berkshire, U.K.) through digestion of the soils ground to 0.5 mm in aqua-regia followed by subsequent analysis on an ICP-AES. The S18OPO4 of the soil was determined on the 1 M HCl extractant described by Tamburini et al. (2010) but with a few modifications. Briefly, between 10 and 20 g dry soil was added to 100 ml 1 M HCl and shaken overnight. The supernatants were collected after separation from the residual solids by centrifugation and filtration. Phosphate concentrations were determined colourimetrically on an Aquachem 250 analyser using a molybdenum blue reaction (Murphy and Riley, 1962) after they were diluted by at least 1/10 to avoid acid interference with the molybdenum chemistry. The extracted PO4 is precipitated and dissolved as firstly ammonium phospho-molybdate and then magnesium ammonium phosphate before excess magnesium and chloride is removed through the addition of a cation resin and a small dose of silver nitrate crystals respectively. The resultant PO4 in solution is then converted to silver phosphate (Ag3PO4) though the addition of an Ag-ammine solution and subsequent adjustment of the pH to between 7 and 8 with 0.5 M HNO3 before incubation for two days at 50 °C in an oven.

Soil PO4 extracted using 1 M HCl (HClPO4) represents an integration of several potential PO4 pools (Zohar et al., 2010): (i) the most labile, water-extractable soil PO4, including intracellular microbial PO4 typically released through processes such as soil drying and re-wetting and cell lysis (ii) weakly adsorbed bicarbonate-extractable soil PO4, and (iii) strongly fixed, calcium and iron bound P. The weak acid extraction does not include the strongly adsorbed aluminium oxide bound PO4, which requires a NaOH-based extraction (Tiessen and Moir, 1993), nor does it include P in organic matter, the hydrolysis of which with 1 M HCl has been shown to be negligible during the development of the extraction protocol by Tamburini et al. (2010). To confirm that organic P forms were not being hydrolysed, duplicate soil samples were extracted using 18O-labeled and unlabelled 1 M HCl. If the extracted HClPO4 was to contain large amounts of hydrolised organic P or condensed PO4 species (e.g. polyphosphates, pyrophosphates, etc.) then the S18OPO4 of the sample duplicates would be markedly different. For all tested samples, the 18O of the sample extracted with unlabelled and 18O-labeled acid was no > 1%. The contribution of calcium-bound P is also considered to be negligible as both the soil parent material was neither igneous nor calcareous, and the soil pH is generally <6. Therefore HClPO4 is assumed to extract PO4 from the same soil pool as is released by water, anion resins, and NaHCO3 extraction, namely extracellular labile PO4 and metabolic intracellular microbial PO4 and also some iron bound PO4. Using the 1 M HCl as an extractant enables far greater quantities PO4 to become available, quantities that are required to allow the successful precipitation of Ag3PO4 from small amount of soil.

Oxygen isotope analysis was carried out ETH Zurich on a Vario Pyro Cube (Elementar GmbH, Hanau, Germany) coupled in continuous flow to an Isoprime 100 isotopic ratio mass spectrometer IRMS (Isoprime, Manchester, UK). Triplicate samples of ~0.3 mg of Ag3PO4 were weighted into silver capsules together with a small amount of fine carbon black powder to promote the reaction between the Ag3PO4 and carbon to produce CO. The produced reaction gases are carried through a

temperature-controlled chromatography column and ultimately to the [RMS. Calibration and corrections for instrumental drifts were done by repeated measurements of an internal standard (ACROS Ag3PO4 > 97.5%, 618O = +14.2%»; measured and calibrated at the University of Lausanne by T. Vennemann), and of benzoic acids IAEA 601 and IAEA 602 ( + 23.3% and +71.4%, respectively). The 618O values are given in the standard % notation with respect to VSMOW (Vienna Standard Mean Oceanic Water). Silver phosphate standards are routinely analyzed, and standard deviation on analytical replicates was better than 0.3%.

23. Statistical analyses

Unfortunately, the described sampling design of Peukert et al. (2012) was flawed due to a mis-interpretation of the requirements needed for an efficient statistical analysis. The design resulted in: (i) large areas of gross under-sampling or voids and (ii) areas of preferential or clustered sampling due to the chosen sampling scheme. Both flaws also interact and compound each other. Given this, the direct application of many statistical methods is likely to be inefficient and suboptimal, resulting in biased outputs (Diggle et al., 2010; Gelfand et al., 2012; Olea, 2007). As a consequence, all of this study's statistical analyses had to be conducted in a manner that would account for this sub-optimal sample design.

As a first step to mitigate against possible biasing effects, the following data pre-processing actions were taken: (1) the limits of the study area were set to lie significantly within the field boundary using a 10 m buffer around the all sample points that lie to the edge; (2) the data were declustered by removing observations that contribute the most to the preferential sampling; and (3) to complement Action 2, one-point declustering weights were found using the algorithm of Deutsch and Journel (1998), that enables the full data set to be used but where the clustered data are down-weighted (and so negate potential bias).

The statistical analyses for assessing the variation in 618OPO4 accorded to the following four linked stages: (i) a standard exploratory analysis; (ii) a variographic analysis to assess spatial dependence (e.g. Goovaerts, 1997); (iii) a geographically weighted (GW) correlation analysis (Fotheringham et al., 2002; Harris and Brunsdon, 2010) to investigate spatial heterogeneity in paired relationships; and (iv) a confirmatory regression analysis based on that observed in stages (i) to (iii), where S18OPO4 is taken as some function of HClPO4, elevation, slope, soil class and the historical north-south division. [n addition to these analyses, that are centred on 618OPO4, complementary spatial analyses were conducted on TP and HClPO4 to provide useful context. All study hypothesis test results are reported at the 95% significance level; and all statistical analyses were implemented in R (

23.1. Standard exploratory analysis

Conditional boxplots and conditional scatterplots were used to assess S18OPO4 distributions and relationships. This included evidence for multiple populations for 618OPO4, according to the historical north-south field division or to soil class; and whether or not the relationship for S18OPO4 with HClPO4, with elevation and with slope; changed according to the north-south division. Multiple linear regression (MLR) fits were conducted using a step-wise, ordinary least squares (OLS) estimation that finds the best predictor variable subset according to the Akaike [nformation Criterion (AIC). Considering the data is spatial, the OLS fits were viewed as exploratory; definitive spatial regression fits are described in Section 2.3.2, below. Hypothesis tests relating to this stage, follow the usual parametric approach via t-tests.

2.3.2. Variographic analysis

To assess spatial dependence in 618OPO4, we limited our investigations to: (i) raw data variograms; (ii) residual data variograms from the OLS MLR trend fits of Section 2.3.1; (iii) (outlier-resistant) robust

variograms (Hawkins and Cressie, 1984); (iv) within-class variograms (e.g. Goovaerts, 1997); (v) local variograms (where (iv) and (v) accord to the historical field division); (vi) normal-scores transformed data variograms; and (vii) cross-variograms with HClPO4. Lags for all empirical variograms were chosen to reflect the study's sample design; and only omni-directional variograms were found. To counter any biasing effects on variogram estimation caused by the sample design, a two-point declustering algorithm was also used to provide weighted variograms (Richmond, 2002). For variograms that displayed spatial structure, they were modelled using a weighted least squares (WLS) approach specified with an exponential variogram model-type.

For investigation (ii), an OLS MLR fit and its corresponding WLS residual variogram model fit are sub-optimal (but often informative in an explorative context). As such (and when required), both sets of parameters (i.e. those for the MLR and those for the variogram) were optimally re-estimated using restricted maximum likelihood (REML) (e.g. Schabenberger and Gotway, 2004). REML is also useful in that it can similarly account for variogram bias due to the sample design, as can the two-point declustering algorithm, above (Marchant et al., 2013).

2.3.3. Geographically weighted correlations

Spatial heterogeneity in S18OPO4 relationships was investigated via the mapping of GW correlations. These localised correlations are found at the sample sites of the study area, where they are calculated in the same manner as their standard (global) counterpart, but only use data that are spatially nearby to the sample locations. Nearby data are spatially weighted (via a kernel weighting function), so that the closest data points are given more influence in the local statistic. Geographically weighted correlations form one of a suite of GW methods (Fotheringham et al., 2002; Lu et al., 2014; Gollini et al., 2015), another of which, GW regression (GWR) is described below.

A GW method can be viewed as a moving window weighting technique, where the size of the window over which any localised statistic/model might apply is controlled by the kernel's bandwidth. Commonly, this exploration of spatial heterogeneity involves: (i) a carefully judged choice ofbandwidth; (ii) a test for significance of the observed heterogeneity; together with (iii) a mapping of the outputs. For this study, we specified the GW correlations using an adaptive Gaussian kernel. [n the absence of an objective procedure for bandwidth selection, we experimented with three user-specified bandwidths of 30%, 40% and 50%. Hypothesis tests for this stage, followed a randomisation approach where the true local correlation is compared to that found from 999 random permutations of the data (e.g. Fotheringham et al., 2002).

2.3.4. Regression analysis

The following three regressions were considered, where S18OPO4 is the response variable: (i) MLR; (ii) MLR with spatially-autocorrelated error; and (iii) GWR (Brunsdon et al., 1996). The first two regressions assume data relationships are globally-fixed, while the third regression allows data relationships to locally-vary. The parameters of the two MLR models can be estimated by OLS and REML, respectively. Similar to GW correlations, GWR enables the exploration of data relationships, where localised regressions are fitted using spatially weighted data and the resultant regression coefficients are mapped. GWR parameters are estimated using a WLS procedure, where we again specify an adaptive Gaussian kernel weighting scheme. For GWR, an objective function exists for bandwidth selection; and here we use an automatic A[C procedure (Fotheringham et al., 2002). To compare the fit of chosen regressions, A[C and R2 values are reported.

Hypothesis tests for significant regression coefficient heterogeneity follow both: (i) a randomisation approach (Brunsdon et al., 1998) and (ii) a parametric bootstrap approach (Harris et al., 2015). For the former, GWR is applied to 999 random permulations of the data and the variance of a given coefficient is found. The actual variance of the same coefficient is then included in the ranked distribution of variances, where

its position (i.e. its p-value) relates to whether there is significant spatial variation in this coefficient.

For the latter, bootstrap samples of the response variable are found for a null (fixed coefficient) model (e.g. MLR), where the simulations are based on the estimated parameters of the null model fit to the sample data. The predictor variables are not considered random and are the same as the sample data. A test statistic q that measures the spatial variability in the GWR coefficients is then used to test against the null hypothesis. Thus a GWR model is fitted to each bootstrap sample and a q-statistic is found for each regression coefficient. As the null model is a random process, even when coefficients do not vary spatially, one would not expect the GWR coefficients to be identical in different locations. The bootstrap analysis determines how much coefficient variability one might expect to encounter due to random variation in a model, and to compare the level of variability in the observed data set, against this.

The bootstrap tests are run with the number of simulations set at 999. For each regression coefficient, the 95% points of the bootstrap samples are computed, and significance levels are found for upper single-tailed hypothesis tests. In addition, at every sample (i.e regression point) location, the localised set of regression coefficients are tested for significant difference to their corresponding fixed coefficient estimates. Here at each sample location, a bootstrap sample of pseudo t-values (e.g. Harris et al., 2010) are found for each coefficient, enabing a bootstrap p-value to be found accordingly. Mapping these p-values identifies where local coefficients significantly differ to their fixed (or global) coefficient counterpart.

3. Results

3.1. Distributions of TP, HClpo4 and 818OPo4 with the raw clustered data

The TP of the soil across Great Field ranged from 736 to 1952 mg P kg-1, while the HClPO4 ranged from 93 to 821 mg P kg-1 extracting between 12 and 48% of the TP present in the sample. As would be expected there was a strong positive correlation between TP and HClPO4 (r = 0.91), where the proportion of TP as HClPO4 was found to increase with increasing soil TP content. The S18OPO4 of the HClPO4 ranged from 17.0 to 21.6% with a mean of 18.8% (±0.8).

The spatial distributions of TP, HClPO4 and S18OPO4 are presented in Fig. 1a-c. The historic field divide appears an important discriminating variable in terms of both TP and HClPO4 with higher values to the north of the divide. Soil class also appears a useful discriminator of TP and HClPO4. The study field slopes downwards from its south-east corner to broadly where the historical divide starts to the west; and similarly

slopes downwards from the north to the same point. Thus the field is broadly concave in shape, and this also appears to influence the distribution of TP and HClPO4. Conversely, there does not appear to be any spatial trend in S18OPO4 or environmental factors that influence its variability. Note that all data descriptions in this section are naive given that any bias due to the sample design, are not (as yet) considered.

3.2. Actions taken to address sub-optimal sample design

As a first action to address potential sub-optimality (Action 1, from Section 2.3), a 10 m buffer was used around the data to limit all analyses to only a sub-region within the study field (Fig. 2). For Action 2 the data were both moderately and strongly declustered by removing 7 and 22 observations, respectively through expert judgement. Actions 1 and 2 are depicted in Fig. 2. Both declustered data sets simply lessen the impact of the three main areas of clustering, where data have been manually removed according to their location (and not by their measurements). The main drawback to the use of declustered data is the reduced information, which is already limited to 78 locations.

Next, we assessed for bias in the (global) means of TP, HClPO4 and 618OPO4, according to the three different data sets depicted in Fig. 2. An alternative set of weighted means were also calculated using weights found from a cell-declustering to a 25 m grid cell (a natural choice given the design); which is Action 3 from before. Results are presented in Table 1, where the clustered data tends to: (i) slightly under-estimate the mean for TP; (ii) slightly over-estimate it for HClPO4; and (iii) very slightly under-estimate it for 618OPO4. Results suggest that continuing with the clustered data is reasonable, with a proviso that only models that cater for possible bias are applied. Weighted correlations and weighted regressions (WLS MLR) can also be found using the clustered data, where we assume that the cell-declustering weights found for an unbiased global mean (for 618OPO4 in Table 1), are also suitable to down-weight data relationships in areas of clustering. It is also prudent to calibrate models with the strongly declustered data; and their outputs compared with models calibrated with the clustered data.

3.3. Complementary analyses for TP and HClPO4

Although our focus is the spatial analysis of 618OPO4, it is useful to provide an insight into how TP and HClPO4 vary spatially. This provides context to the analysis of 618OPO4, especially as we choose to investigate its relationship to HClPO4 in subsequent sections. As opposed to TP, HClPO4 more directly relates to our understanding of 618OPO4. Thus in Fig. 3, prediction surfaces are given for TP and HClPO4; both of which were constructed using a kriging with external drift (KED) model,

Fig. 1. The spatial distribution of (a) TP, (b) HClPO4, and (c) 81sOPO4 of the HClPO4 within 'Great Field' on the North Wyke Farm Platform.

Fig. 2. The distribution of: (a) clustered data, (b) moderately declustered data and (c) strongly declustered data.

with all parameters estimated optimally via REML, and a global neighbourhood specified. For TP, its KED trend component was informed by elevation, slope and the north/south division; whereas for HClPO4 its trend was informed by elevation, slope, soil class and the north/south division. The R2 values for the trend fits (i.e. via OLS) were strong at 0.69, in both cases.

Observe that for the trend component of each KED model, TP and HClPO4 were not used to help predict each other, as this would not inform predictions on a grid (618OPO4 was similarly not used in this respect). Variography for both KED models was reasonably well-behaved, with a clear single structure depicting spatial dependence up to distances of around 100 m. As both variables are highly correlated, their parameterisation and surfaces are broadly similar. The highest TP and HClPO4 values clearly lie to the north of the historical division, while the lowest values lie in a broad swathe south and east of the historical division.

All analyses in this section were conducted on the clustered data, but with consideration to potential bias. Observe that the clustered data can be used when kriging (i.e. after its parameterisation), as it inherently accounts for such configurations via its information, screening and relay effects (Chiles and Delfiner, 1999).

3.4. Exploratory analysis for S18OPO4

Using both the clustered and strongly declustered data, the relationship matrices for 618OPO4, HClPO4, elevation and slope are presented in Fig. 4a & b. Correlation coefficents are given for the complete data sets, while scatterplots depict relationships that are conditional to the historic field division. The latter of which, provides a first insight into possible local relationships. Weighted correlations using the cell-declustering weights are also given in Table 2. By focusing only on those relationships with 618OPO4, the strongest relationship is with HClPO4, but this is weak with a correlation of only 0.30 (for the weighted correlation). Clearly, S18OPo4 poorly correlates with elevation and with slope. It appears that data relationships may be conditional to the historic field division; and in particular, the relationship of S18OPO4 with HClPO4.

Table 1

Global means of clustered and declustered data sets.

Clustered Moderately Strongly Clustered data

Variable data declustered data declustered data with weights

TP 1134.82 1146.49 1151.94 1145.61

HClPo4 317.77 322.37 314.73 311.79

S18Opo4 18.75 18.81 18.89 18.62

Conditional boxplots are used in Fig. 4c & d to relate the distribution of S18OPO4 to both the historic field division and to soil class; again using both clustered and strongly declustered data. Marginally higher S18OPO4 values are generally found in the northern part of the field, but in general, discrimination is poor. Similarly, the three soil classes do not appear to be a strong discriminator of 618OPO4. In general, there is little to choose between the clustered and declustered data analyses. Thus the sample design does not appear to strongly bias S18OPO4 in this respect. Regardless of any biasing effects, all relationships for S18OPC>4 are weak or indistinct. Locally however, this may not be the case, and we investigate this further in Sections 3.6 and 3.7.

Tables 3 and 4 present the results from a series of MLR fits to the clustered and strongly declustered data, with S18OPO4 as the response. Table 3 reports the results using all predictors (HClPO4, elevation, slope, soil class and north/south division), while Table 4 reports the results using predictor subsets chosen via step-wise AIC. In both cases, WLS MLR fits are also reported using the cell-decustering weights from before. Clearly, all MLR fits are poor, where the highest R2 is only 0.35. There is also no consistency in the make-up of the predictor subset or in the significance of those variables chosen (note that significance tests are naive in that any spatial dependence in the data is not as yet considered). For the former, this is not surprising given that reductions in AIC are consistently small. Unlike the previous analyses, it now appears that the sample design is detrimental to fitting regressions (as the poorest R2 values result when directly using the clustered data). None of the predictor variables appear particularly worthy predictors of 618OPO4 in this global sense, but it is unclear which, if any varables, could be safely removed, before we proceed to more detailed analyses. Given these analyses, study hypothesis (A) is accepted.

3.5. Variographic analysis for 61SOPO4

We next investigate for spatial dependence in 618OPO4 using raw data, residual data, robust, local and within-class empirical variograms. Again, we compare results using the clustered and strongly declustered data. For the residual variograms, OLS MLR trend fits with all predictor variables are used. We also provide a weighted variogram to the clustered data and REML model fits to both data sets using all available predictors to inform the trend. All variograms are given in Fig. 5, where only very weak evidence for spatial dependence is found and the S18OPO4 data is essentially random. This is endorsed by the REML results, where the AIC for the spatial model was 4 units higher than the corresponding non-spatial model for the both data sets. Unsurprisingly, as spatial dependence was absent in 618OPO4, cross-dependence with HClPO4 was also absent (even though spatial dependence is present in

Fig. 3. The spatial distribution of (a) TP and (b) HClPO4; each found using a KED model.

HCiPo4, from Section 3.3). Given these results, it is unnecessary to for- From Fig. 5, the effects of data clustering on the variograms are quite

mally test for spatial dependence in S18OPO4 and study hypothesis (B) is evident, where semi-variances at the lower lags of the clustered data accepted. variograms are commonly heightened, resulting in very poor small-

(â) Conditional scatterplots (red = south; blue = north) - clustered data

140 145 160 155

17 18 19 20 21

Slope -0.26 0.19 0.15

• •• • j • >;>./-• : Elevation -0.28 -0.0014

• • • * • . t / • s • * »V • , • •» . • • • ; - - HCIpo4 0.24

.vv-'viv S18Op04

(b) Conditional scatterplots (red = south; blue = north) - declustered data

140 145 150 155

17 10 19 20 21

100 300 500 700

Slope -0.23 0.11 0.17

• • • » • .1- Elevation -0.25 0.026

• * • •. • • •. » • ••• •.....» . HCIp04 0.29

•• ••¿i.'^v. • • •• ♦ • »7 • • • 518Opo4

(c) Conditional boxplots - clustered

(d) Conditional boxplots - declustered

Fig. 4. Scatterplots and correlation coefficients for the (a) clustered and (b) strongly declustered data, respectively. These display the pairwise relationships between S1sOPO4, HClPO4, elevation and slope. Points coloured red and blue relate to samples to the south and north of the historic field division, respectively. Conditional boxplots for S1sOPO4, with respect to the historic field division or the three soil classes for (c) clustered and (d) strongly declustered data, respectively.

Table 2

Weighted correlations using the cell-declustering weights.

Slope Elevation HC1po4 S18Opo4

Slope 1 -0.30 0.14 0.19

Elevation - 1 - 0.26 - 0.019

HClpo4 - - 1 0.30

818Opo4 - - - 1

scale structure. Conversely, most of the strongly declustered data variograms depict small-scale structure, but no longer have any semivariance values at the lowest lag distances. The differences between the raw data variograms (using complete data sets) and the south data variograms indicate that 618OPo4 has a higher variance in the north than that found in the south. The behaviour of the within-class variograms somewhat endorses this (noting that they mimic the south variograms at higher lags as no such semivariance pairs exist in the north). North data variograms were not found as there were too few data points to compute a valid variogram. All residual variograms clearly depict pure nugget effects, as do the REML models. Finally, normal-scores variograms were also found, but did not help in structure identification.

3.6. Local relationships for S18OPO4

Given that spatial autocorrelation effects are absent for the S18OPO4 process, out next objective is to assess for any spatial heterogenic effects with respect to the relationships between 618OPO< with HClPO4, elevation, and slope. Furthermore, if such heterogeneities exist, do they depend on the historical field divison or alternatively, the under-lying soil class?

Evidence for such heterogeneities have already been observed in Fig. 4, where the relationship between 618OPO< and HClPO4 appears conditional to the historic field division. In order to explore this particular relationship further, the same conditional scatterplot is given in Fig. 6a, using the strongly declustered data. Here we calculate correlations for the northern and southern parts of the field. Given that spatial autocorrelation effects are absent (both globally and within each partition), and that the declustered data is used, we can report the significance ofthese correlations (using standard t-tests with correlations) with relative assurance. Here the global correlation of 0.29 between HClPO4 and S18OPO4 is significantly different to a zero correlation, whereas the local correlation in the south is very weak at 0.02 and not significant. In the north, the local correlation was relatively strong at 0.55, but given that it is found using only eleven data points, this correlation was also insignificant. Observe that we have shown the data with their linear fits, where the intercept and slope vary locally. This concept of spatially-varying regression coefficients is re-visited in Section 3.7. Note that as we are investigating locally, it is assumed that an unbiased analysis will result by only using the strongly declustered data (see Section 4).

To further our investigations of relationship heterogeneity, Fig. 6b-d presents GW correlation maps for S18OPO4 with HClPO4, with elevation, and with slope; each specified with a bandwidth of 40%. Co-variability for S18OPO4 with HClPO4 tends to be stronger in the northern part of the field (Fig. 6b), and these findings are endorsed by the associated randomisation tests that indicate areas of unusually strong correlation, as well as unusually weak correlation in the south. The distribution of GW correlations between HClPO4 and 618OPO4 largely confirm that

found in Fig. 6a, but with more detail. As a continuous (Gaussian) weighting scheme is specified, all GW correlations are actually informed by all 58 observations of the strongly declustered data (a bandwidth of 40% entails that the nearest 22 observations exert the greatest influence on each localised correlation). This weighting specification is considered crucial to a GW analysis to such a relatively small data set. Thus the GW correlations in the north and the south of the field are better informed than the two partitioned correlations from before. From Fig. 6c-d, correlations for S18OPO4 with elevation, and with slope, only marginally vary across the field. It is possible that all such heterogeneities are dependent on the field division (and therefore different managements) or the under-lying soil class.

3.7. Regression analysis for S18OPO4

Given the analyses of the preceeding sections, only MLR (estimated using OLS) and GWR models need to be considered. Extending MLR to account for spatial autocorrelation in the residual term is unnecessary given the residual data variograms and REML results of Section 3.5. Models are applied, using the strongly declustered data only and all predictor variables are used. Thus the MLR model is the same as that given in Table 3. The bandwidth for the GWR model is optimally found at 88%, indicating a shallow weighting and suggesting weak relationship heterogeneity. The model fit results for MLR and GWR give AIC values of 125.4 and 124.6, respectively; together with R2 values of 0.35 and 0.38, respectively. Clearly, there is little to be gained from applying GWR in preference to MLR. Table 5 tests for coefficient heterogeneity using: (i) the parametric bootstrap approach with MLR as the null model and (ii) the randomisation approach. Clearly, and as expected, there is no significant evidence for coefficient (i.e. relationship) heterogeneity in this data, as all p-values are >0.05. Given these test results, study hypothesis ( C) is accepted.

The localised regression coefficients from GWR associated with HClPO4 are mapped in Fig. 7a, where it appears that the relationship for S18OPO4 with HClPO4 varies spatially. Furthermore, this relationship appears to depend on the historical field divison and/or the under-lying soil class. There is also a north-west to south-east trend in the coefficients, reflecting the shallow weighting function that was specified. In Fig. 7b, the bootstrap p-values for the local sets of pseudo t-values indicate where the coefficients associated with HClPO4 are significantly smaller (to the north) and significantly larger (to the south) than the global value. However from the map legend, this only relates to a 90% significance level and not the stated study significance level of 95%. Thus, the local coefficients are not significantly different to their global counterpart, and as such, study hypothesis (D) is accepted. Observe that given the results of Table 5, we do not further investigate the localised regression coefficients for the intercept or the other predictor variables.

4. Discussion

4.1. Variability ofHClPC,4 and TP

Analyses for HClPO4 and TP indicate that both variables are spatially autocorrelated; where elevation, slope, soil class and the historic field divide can influence this variability. However the historic field divide is considered the most important driver, with higher values to the north of the divide than to the south. This is not unexpected given the management differences between the two sectors and is because the

Table 3

MLR fits using the full predictor variable set.

Data set Estimation R2 Predictor set Significant predictors

Clustered OLS 0.23 All predictors Intercept, Hallsworth soil class

Clustered WLS 0.31 All predictors Intercept, Hallsworth soil class

Strongly declustered OLS 0.35 All predictors Intercept, HClp04, Hallsworth soil class

Table 4

MLR fits using predictor variable subsets chosen by the step-wise AIC procedure.

Data set Estimation R2 AIC reduction Predictor subset Significant predictors

Clustered OLS 0.22 4.6 Slope, soil class Intercept, slope, soil class

Clustered WLS 0.31 3.3 HClPO4, slope, soil class Intercept, slope, Hallsworth soil class

Strongly declustered OLS 0.34 1.2 HClPO4, elevation, north/south division, soil class Intercept, HClpo4, elevation, Hallsworth soil class

northern part of the field has not been ploughed for a significant time, while the southern part had been ploughed and most recently in 2007. The surface enrichment of P is normal in agricultural soils, especially grasslands which may not be ploughed frequently, and reflects the accumulation of P from inorganic fertilizers and manures over time. The decrease in P with depth is often marked over a few cm with the bulk of the soil profile considerably lower in P. Where ploughing has occurred, such profiles are destroyed within the plough zone and the high P surface soil diluted or replaced with lower P soil brought up from depth (e.g. Haygarth et al., 1998a, b; Watson and Matthews, 2008). The amounts of the HClPO4 indicate that > 50% of P in the soil is in 1 M HCl recalcitrant forms such as organic P or aluminium oxides. Given the negligible contribution of bedrock apatite PO4 it is assumed that the bulk of the HClPO4 is comprised of PO4 adsorbed to soil particles and microbial PO4. The soils in this region have been shown to have a high microbial biomass P content, which in turn has been shown to be released upon drying and rewetting through cell lysis (Blackwell et al., 2013). The amounts of HClPO4 measured are very similar to the levels of microbial P reported by Blackwell et al. (2013) and therefore we assume that HClPO4 in these soils represents, in large part, adsorbed soil PO4 and intracellular microbial PO4, the latter normally exceeding the former (Blackwell et al., 2010).

42. Variability of S18OPO4

In order to discuss S18OPO4 values observed, we need to estimate the theoretical equilibrium S18OPO4 value that might be expected for PO4 in equilibrium with soil water which is believed to be mediated by the ubiquitous intracellular enzyme pyrophosphatase (Blake et al., 2005). This causes the exchange of PO4 oxygen with the oxygen in H2O and results in a temperature dependent relationship initially described by Longinelli and Nuti (1973). However recent work by Chang and Blake (2015) has developed a refined, rigorous and controlled laboratory calibration of the temperature-dependence of equilibrium PO4 and water, catalyzed by pyrophosphatase, over a typical environmental temperatures (3-37 °C):

ES18OPO4= - 0.18 T +26.3 + S18OH2Owhere ES18OPO4 is the stable oxygen isotope ratio of PO4 at equilibrium in %», T is the temperature in degrees Celsius and S18OH2O is the stable oxygen isotope ratio of H2O in %». The intracellular phosphate, already at equilibrium, is released to the soil after cell lysis.

Although measurements of soil water S18OH2O were not directly made it has been suggested that soil water at the location is very similar S18OH2O to ground water, in that it is an integrated value of many rainfall events (Granger et al., 2010). As such the S18OH2O should be similar to that of the global meteoric water line for this area with S18OH2O ranging between — 5.5 to — 6.0% (Darling et al., 2003). The average soil temperature measured at 10 cm depth nearby was 13 °C for the month of May and this would give estimated vales of S18OPO4 at equilibrium with soil water of between 18.0 and 18.5%. However given the uncertainties in these assumptions it is wise to examine a window of potential soil S18OPO4 and temperature values to understand better the measured S18OPO4 data compared to the theoretical equilibrium 618OPO4 values. For example, it might be expected that the integrated soil temperature of the soil profile from surface to 7.5 cm depth is slightly warmer than that measured at 10 cm. Further, although ground water is an analogue for soil water, the variability of soil water S18OH2O is more pronounced and subject to strong influence from individual rainfall events and also from enrichment in 18O through evapotranspiration, especially in the surface of the soil (e.g. Hsieh et al., 1998; Treydte et al., 2014). Therefore if it is assumed that the integrated 7.5 cm profile soil temperature ranges between 13 and 15 °C, and that the S18OH2O is enriched by up to 4% compared to groundwater, then the theoretical S18OPO4 equilibrium values range between 18.0 and 23.0%. This range of values better matches the S18OPO4 of the soil and suggests that the HClPO4 within the soil is at or around equilibrium. This finding is in common with other researchers. Angert et al. (2012) found that across a Mediterranean bed rock and rainfall gradient, both soil resin extractable PO4 and HClPO4 were broadly in equilibrium with soil water and proposed, amongst other reasons, that a flux of intracellular equilibrated PO4 from the soil microbial biomass may be causing this. Tamburini et al. (2012) also conclude that in young and developing alpine soils, regardless ofthe contribution of PO4 from minerals or vegetation, or from the activity of

Fig. 5. Variograms for 81SOPO4 using (a) clustered data and (b) strongly declustered data.

Fig. 6. (a) Strongly declustered data conditional scatterplot for 81SOPO4 with HClPO4. Points are separated on their location in the field relative to the historic field division. Lines of best fit, using: all of the data, north data only, and south data only are shown. Geographically weighted correlations for 81SOPO4 with (b) HClPO4, (c) elevation, and (d) slope - at the strongly declustered data locations. Associated randomisation test results are circled, indicating locations of unusual correlation. Maps are given with the historical field division and the three soil class divisions.

extracellular enzymes, microbial biomass was the main controller ofthe S18OPO4, keeping it in equilibrium with soil water. In both these examples however the soils were generally low in P, and where P is limiting it is expected that microbial cycling would be rapid. Within this current study soil P should not be limiting and it is therefore interesting that S18OPO4 would still appear to be relatively rapidly cycled and have no obvious trace of any P source signatures. One explanation for this, would be if the predominant P sources had S18OPO4 signatures similar to that of the equilibrium 618OPO4 value. The 618OPO4 of inorganic fertilizers have been reported to show a very wide range of values from 14.8 to

27.0% while virtually no S18OPO4 values have been published for animal excreta (Young et al., 2009) or for stored and managed animal wastes, however, water-extractable PO4 from dairy farm slurries in this area have been found to range from 12.0 and 15.0% (Granger, unpublished data).

Analysis results for 618OPO4 relationships within the data are enigmatic. Considering the field as a whole, all relationships for S18OPO4 are weak or indistinct. HClPO4, elevation, slope, soil class and the historic field divide are not particularly good drivers of S18OPO4 variability; and S18OPO4 is not spatially correlated with itself. For within-field


Bootstrap test q-statistics, associated p-values and randomisation test p-values.

Intercept HClP04 Elevation Slope N/S division Hallsworth soil class Halstow soil class

Actual 0.726 0.000 0.006 0.006 0.091 0.026 0.059

MLR 95% 2.475 0.001 0.017 0.026 0.242 0.078 0.128

MLR p-value 0.266 0.140 0.220 0.400 0.208 0.193 0.158

Randomisation test p-value 0.275 0.096 0.200 0.575 0.443 0.900 0.584

Fig. 7. The (a) GWR coefficients for HClPO4 with a bandwidth of 88% and (b) the associated results from parametric bootstrap test.

relationships between S18OPO4 and HClPO4, elevation and slope; only the local relationship for 618OPO4 with HClPO4 appears to have any value, showing a positive relationship north of the historic field divide which does not exist to the south. However, such localised relationships are not significantly different to that found for the field as a whole. Thus in essence, S18OPO4 is not only randomly distributed with respect to itself, but also to the particular environmental factors that this study has considered.

If it had been found that the positive relationship between 618OPOi with HClPO4 was significant in the northern part of the field, which this study can only allude to; this may imply that the process of ploughing, which has destroyed the soil P profile in the south, has also decoupled any link between HClPO4 and 618OPO4 that is present in the north. The scenario described by the data seems to suggest that, in the north, elevated HClPO4 occurs as 'hotspots' with associated higher 618OPO4, and that as the HClPO4 of these hotspots declines so the 618OPO4 becomes lower. Such P hotspots have been described previously under cattle grazing regimes which not only leads to an increasing the TP in the surface of the soil but also to a heterogeneous distribution of TP 'hotspots' as a result of livestock excreta (Page etal., 2005; Penn etal., 2007). However, the available data on S18OPO4 in excreta, although limited, indicates that it is not elevated. If however metabolic PO4 released from leaf litter has an elevated S18OPO4 (Pfahler et al., 2013) this might explain the difference between the north and south of the field. In the south, ploughing has buried and mixed soil organic matter, which in the north has had time to accumulate. Potentially, this accumulation may not be homogenous and may still be affected by the hotspots of excreta described previously or it might just reflect a heterogeneous accumulation of non-excretal plant material within the soil profile. However if this is the cause of the relationship between 618OPO4 and HClPO4 it is not reflected in the soil total carbon content which shows no relationship with S18OPO4.

4.3. Controling variables

Throughout this study, we have not attempted to de-couple the effects of the historic field divide from the underlying soil class, in how these categorical variables influence variation and co-variation in HClPO4, TP and 618OPO4. It is clear from Fig. 1 that both categorical variables act as similar discriminators, as data only coincides with the Denbigh class in the northern part of the field, while the Halstow class primarily relates to data in the south. Such similarity in discriminating power is then clearly evident in the subsequent analyses. Given that

sample size is small and many results are null, it was considered appropriate to maintain the controlling interaction between these categorical variables in all regression fits. However, as the distinction between the three soil classes is considered minimal in practise, any HClPO4, TP or 618OPO4 difference between the soil classes is considered more likely a reflection of the historic field division rather than a characteristic of the soil classes themselves.

4.4. Analysis limitations

Although the statistical analyses for S18OPO4 point to a series of null results, it should be borne in mind that the study data is limited in size and that the sampling configuration may miss the key scales of spatial dependence and co-dependence in 618OPO4. In this respect, future 618OPO4 (and HClPO4) studies should consider sampling on a finer grid with increased sampling, with nested sampling strategies (e.g. Atteia et al., 1994; Corstanje et al., 2008). Increased sampling should also guard against results simply reflecting sampling variation rather than the true properties of the spatial process.

Issues of sample design necessitated the use of the strongly declustered data for all localised analysis; and the associated loss of information. The only way to have proceeded with the clustered data would have required the calculation of localised sets of de-clustering weights, as it was not viable to use the globally-found ones of Section 3.2. As this would have presented a challenge for any partitioned analysis, and more so for any GW analyses (where the nature of the clustering bias would locally-vary according to the kernel and its bandwidth), the use of localised de-clustering weights was not considered.

Given the poor model fits throughout, the small sample size and the known issues for GWR with respect to local predictor variable collinear-ity (Paez et al., 2011), the results of Section 3.7 are viewed with caution. Although here, the GWR fit was assessed for adverse collinearity effects following the procedures outlined in Gollini et al. (2015), and although collinearity was present, its effects were not considered serious. Further work could extend the GWR analyses in this respect, where in addition to the fit of a penalised form of GWR, a mixed GWR (Brunsdon et al., 1999) could be considered where the categorical predictors are globally fixed, while the other predictors are allowed to locally-vary.

It is also worth noting that variation in S18OPO4 may be driven by environmental covariates not considered in this study, as such, future studies should consider identifying these missing covariates. Measurement error in 618OPO4 could also be mitigating factor in this study's

null results and could be incorporated in future work, noting that the S18OpO4 values of this study represent an average of a three replicates.

5. Conclusions

This study has shown that, despite differences in elevation, slope, the under lying soil class and management, the 618OPO4 has no relationship with these field variables even though the HClPO4 can have. The soil 618OPO4 values within the field were found to be within the range predicted for that has been microbially cycled and is at equilibrium with 618OHjO. However given the lack of information of PO4 source 618OPO4 signatures, it is not clear whether this is due to the rapid microbial cycling of PO4 or because PO4 sources have a similar 618OPO4 to that of the equilibrium value. Although no relationships were found between the 618OPo4, and itself (with respect to spatial autocorrelation) and the field variables investigated, there was a suggestion of a positive correlation between 618OPO4 and HClPO4 in the northern unploughed, field sector which was not present in the southern ploughed part of the field. While no conclusive evidence has been found to explain this, it has been suggested that this might be related to the metabolic PO4 within plant material heterogeneously accumulating in the soil surface in the unploughed part of the field.

This study provides an important advance to understanding spatial dependencies and spatial relationships in phosphate stable oxygen isotopes within a temperate agricultural soil. Some variability information was as hypothesised while intriguingly others were not. Given this, further work is ear-marked to learn and benefit from this study's results, via the implementation of a coherent multivariate sampling design, to a different field, with a known long-term management history.


Research is supported by the Biotechnology and Biological Sciences Research Council (BBSRC BB/J004308/1) for the North Wyke Farm Platform. The authors also wish to thank the anonymous reviewers and editor for their time improving this manuscript.


Amberger, A., Schmidt, H.L., 1987. The natural isotope content of nitrate as an indicator of its origin. Geochim Cosmochim. Acta 51 (10), 2699-2705.

Angert, A., Weiner, T., Mazeh, S., Sternberg, M., 2012. Soil phosphate stable oxygen isotopes across rainfall and bedrock gradients. Environ. Sci. Technol. 46 (4), 2156-2162.

Atteia, O., Dubois, J.-P., Webster, R., 1994. Geostatistical analysis of soil contamination in the Swiss Jura. Environ. Pollut 86,315-327.

Baker, A., 2002. Fluorescence properties of some farm wastes: implications for water quality monitoring. Water Res. 36 (1), 189-195.

Bilotta, G.S., Brazier, R.E., Haygarth, P.M., 2007. The impacts of grazing animals on the quality of soils, vegetation, and surface waters in intensively managed grasslands. Advances in AgronomyAdvances in Agronomy Vol. 94. Elsevier Academic Press Inc., San Diego, pp. 237-280.

Blackwell, M.S.A., Brookes, R.C., de la Fuente-Martinez, N., Gordon, H., Murray, P.J., Snars, K.E., Williams, J.K., Bol, R., Haygarth, P.M., 2010. Phosphorus solubilization and potential transfer to surface waters from the soil microbial biomass following drying-rewetting and freezing-thawing. In: Sparks, D.L. (Ed.), Advances in AgronomyAdvances in Agronomy Vol 106. Elsevier Academic Press Inc., San Diego, pp. 1-35.

Blackwell, M.S.A., Carswell, A.M., Bol, R., 2013. Variations in concentrations of N and P forms in leachates from dried soils rewetted at different rates. Biol. Fertil. Soils 49 (1), 79-87.

Blake, R.E., O'Neil, J.R., Surkov, A.V., 2005. Biogeochemical cycling of phosphorus: insights from oxygen isotope effects of phosphoenzymes. Am. J. Sci. 305 (6-8), 596-620.

Brunsdon, C., Fotheringham, A.S., Charlton, M.E., 1996. Geographically weighted regression: a method for exploring spatial nonstationarity. Geogr. Anal. 28 (4), 281-298.

Brunsdon, C., Fotheringham, A.S., Charlton, M.E., 1998. Geographically weighted regression - modelling spatial non-stationarity. J. R. Stat. Soc. Ser. D Stat. 47 (3), 431-443.

Brunsdon, C., Fotheringham, A.S., Charlton, M.E., 1999. Some notes on parametric signficance tests for geographically weighted regression. J. Reg. Sci. 39 (3), 497-524.

Chang, S.J., Blake, R.E., 2015. Precise calibration of equilibrium oxygen isotope fraction-ations between dissolved phosphate and water from 3 to 37°C. Geochim. Cosmochim. Acta 150, 314-329.

Chiles, J.P., Delfiner, P., 1999. Geostatistics - Modelling Spatial Uncertainty. Wiley, New York.

Collins, A.L., Walling, D.E., Leeks, G.J.L., 1997. Sediment sources in the Upper Severn catchment: a fingerprinting approach. Hydrol. Earth Syst. Sci. 1 (3), 509-521.

Corstanje, R., Kirk, G.J.D., Pawlett, M., Lark, R.M., 2008. Spatial variation of ammonia volatization from soil and its scale-dependent correlation with soil properties. Eur. J. Soil Sci. 59,1260-1270.

Darling, W.G., Bath, A.H., Talbot, J.C., 2003. The O & H stable isotopic composition of fresh waters in the British Isles. 2. Surface waters and groundwater. Hydrol. Earth Syst. Sci. 7 ( 2), 183-195.

Deutsch, C.V., Journel, A.G., 1998. GSLIB: Geostatistical Software Library and User's Guide. Oxford University Press, New York.

Diggle, P.J., Menezes, R., Su, T., 2010. Geostatistical inference under preferential sampling. J. R. Stat. Soc. 59,191-232.

Fotheringham, A.S., Brunsdon, C., Charlton, M., 2002. Geographically Weighted Regression - The Analysis of Spatially Varying Relationships. Wiley-Blackwell, Chichester.

Gelfand, A.E., Sahu, S.K., Holland, D.M., 2012. On the effect ofpreferential sampling in spatial prediction. Environmetrics 23 (7), 565-578.

Glendell, M., Granger, S.J., Bol, R., Brazier, R.E., 2014. Quantifying the spatial variability of soil physical and chemical properties in relation to mitigation of diffuse water pollution. Geoderma 214,25-41.

Gollini, I., Lu, B., Charlton, M., Brunsdon, C., Harris, P., 2015. GWmodel: an R package for exploring spatial heterogeneity using geographically weighted models. J. Stat. Softw.63 (17), 1-50.

Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Applied Geostatistics SeriesOxford University Press, New York.

Goovaerts, P., 1998. Geostatistical tools for characterizing the spatial variability of microbiological and physico-chemical soil properties. Biol. Fertil. Soils 27 (4), 315-334.

Granger, S.J., Bol, R., Meier-Augenstein, W., Leng, M.J., Kemp, H.F., Heaton, T.H.E., White, S.M., 2010. The hydrological response of heavy clay grassland soils to rainfall in south-west England using 82H. Rapid Commun. Mass Spectrom. 24 (5), 475-482.

Granger, S.J., Heaton, T.H.E., Bol, R., Bilotta, G.S., Butler, P., Haygarth, P.M., Owens, P.N., 2008. Using 815N and 81sO to evaluate the sources and pathways of NO3- in rainfall event discharge from drained agricultural grassland lysimeters at high temporal resolutions. Rapid Commun. Mass Spectrom. 22 (11), 1681-1689.

Harris, P., Brunsdon, C., 2010. Exploring spatial variation and spatial relationships in a freshwater acidification critical load data set for Great Britain using geographically weighted summary statistics. Comput Geosci. 36 (1), 54-70.

Harris, P., Brunsdon, C., Gollini, I., Nakaya, T., Charlton, M., 2015. Using bootstrap methods to investigate coefficient non-stationarity in regression models: an empirical case study. Procedia Environ. Sci. 27,112-115.

Harris, P., Fotheringham, A.S., Juggins, S., 2010. Robust geographically weighed regression: a technique for quantifying spatial relationships between freshwater acidification critical loads and catchment attributes. Ann. Assoc. Am. Geogr. 100 (2), 286-306.

Harrod, T.R., Hogan, D.V., 2008. The soils of North Wyke and Rowden Unpublished report to North Wyke Research, revised edition of original report by In: Harrod, T.R (Ed.), Soil Survey of England and Wales (1981). North Wyke Research, Okehampton, Devon.

Hawkins, D.M., Cressie, N., 1984. Robust Kriging - a proposal. J. Int. Assoc. Math. Geol. 16 (1), 3-18.

Haygarth, P.M., Chapman, P.J., Jarvis, S.C., Smith, R.V., 1998a. Phosphorus budgets for two contrasting grassland farming systems in the UK. Soil Use Manag. 14,160-167.

Haygarth, P.M., Hepworth, L, Jarvis, S.C., 1998b. Forms of phosphorus transfer in hydro-logical pathways from soil under grazed grassland. Eur. J. Soil Sci. 49 (1), 65-72.

Haygarth, P.M., Jarvie, H.P., Powers, S.M., Sharpley, A.N., Elser, J.J., Shen, J.B., Peterson, H.M., Chan, N.I., Howden, N.J.K., Burt, T., Worrall, F., Zhang, F.S., Liu, X.J., 2014. Sustainable phosphorus management and the need for a long-term perspective: the legacy hypothesis. Environ. Sci. Technol. 48 (15), 8417-8419.

Heathwaite, A.L., Dils, R.M., 2000. Characterising phosphorus loss in surface and subsurface hydrological pathways. Sci. Total Environ. 251, 523-538.

Hsieh, J.C.C., Chadwick OA, Kelly, E.F., Savin, S.M., 1998. Oxygen isotopic composition of soil water: quantifying evaporation and transpiration. Geoderma 82 (1-3), 269-293.

Lewis, W.M., Morris, D.P., 1986. Toxicity of nitrite to fish - a review. Trans. Am. Fish. Soc. 115 (2), 183-195.

Longinelli, A., Nuti, S., 1973. Oxygen isotope measurments of phosphate from fish teeth and bones. Earth Planet. Sci. Lett. 20 (3), 337-340.

Lu, B., Harris, P., Charlton, M., Brunsdon, C., 2014. The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models. Geo-spat. Inf. Sci. 85 (17), 85-101.

Marchant, B., Viscarra Rossel, R.A., Webster, R., 2013. Fluctuations in method-of-moments variograms caused by clustered sampling and their elimination by declustering and residual maximum likelihood estimation. Eur. J. Soil Sci. 64, 401-409.

Marriott, C.A., Hudson, G., Hamilton, D., Neilson, R., Boag, B., Handley, L.L., Wishart, J., Scrimgeour, C.M., Robinson, D., 1997. Spatial variability of soil total C and N and their stable isotopes in an upland Scottish grassland. Plant Soil 196 (1), 151-162.

Murphy, J., Riley, J.P., 1962. A modified single solution method for determination of phosphate in natural waters. Anal. Chim. Acta 26 (1), 31 -36.

Old, G.H., Naden, P.S., Granger, S.J., Bilotta, G.S., Brazier, R.E., Macleod, C.J.A., Krueger, T., Bol, R., Hawkins, J.M.B., Haygarth, P., Freer, J., 2012. A novel application of natural fluorescence to understand the sources and transport pathways of pollutants from livestock farming in small headwater catchments. Sci. Total Environ. 417,169-182.

Olea, R.A., 2007. Declustering of clustered preferential sampling for histogram and semivariogram inference. Math. Geol. 39, 453-467.

Paez, A., Farber, S., Wheeler, D., 2011. A simulation-based study of geographically weighted regression as a method for investigating spatially varying relationships. Environ. Plan. A 43 (12), 2992-3010.

Page, T., Haygarth, P.M., Beven, K.J., Joynes, A., Butler, T., Keeler, C., Freer, J., Owens, P.N., Wood, G.A., 2005. Spatial variability ofsoil phosphorus in relation to the topographic

index and critical source areas: sampling for assessing risk to water quality. J. Environ. Qual. 34 (6), 2263-2277.

Penn, C.J., Bryant, R.B., Needelman, B., Kleinman, P., 2007. Spatial distribution of soil phosphorus across selected New York dairy farm pastures and hay fields. Soil Sci. 172 (10), 797-810.

Peukert, S., Bol, R., Roberts, W., Macleod, C.JA, Murray, P.J., Dixon, E.R., Brazier, R.E., 2012. Understanding spatial variability of soil properties: a key step in establishing field- to farm-scale agro-ecosystem experiments. Rapid Commun. Mass Spectrom. 26 (20), 2413-2421.

Pfahler, V., Durr-Auster, T., Tamburini, F., Bernasconi, S., Frossard, E., 2013.18O enrichment in phosphorus pools extracted from soybean leaves. New Phytol. 197 (1), 186-193.

Richmond, A., 2002. Two-point seclustering for weighted data pairs in experimental variogram calculation. Comput. Geosci. 28 (2), 231-241.

Rivero, R.G., Grunwald, S., Osborne, T.Z., Reddy, K.R., Newman, S., 2007. Characterization of the spatial distribution of soil properties in water conservation area 2A, Everglades, Florida. Soil Sci. 172 (2), 149-166.

Schabenberger, O., Gotway, CA, 2004. Statistical Methods for Spatial Data Analysis. Taylor & Francis.

Smith, V.H., Tilman, G.D., Nekola,J.C., 1999. Eutrophication: impacts of excess nutrient inputs on freshwater, marine, and terrestrial ecosystems. Environ. Pollut. 100 (1-3), 179-196.

Tamburini, F., Bernasconi, S.M., Angert, A., Weiner, T., Frossard, E., 2010. A method for the analysis of the delta O-18 of inorganic phosphate extracted from soils with HCl. Eur. J. Soil Sci. 61 (6), 1025-1032.

Tamburini, F., Pfahler, V., Bunemann, E.K., Guelland, K., Bernasconi, S.M., Frossard, E., 2012. Oxygen isotopes unravel the role of microorganisms in phosphate cycling in soils. Environ. Sci. Technol. 46 (11), 5956-5962.

Tamburini, F., Pfahler, V., von Sperber, C., Frossard, E., Bernasconi, S.M., 2014. Oxygen isotopes for unraveling phosphorus transformations in the soil-plant system: a review. Soil Sci. Soc. Am. J. 78 (1), 38-46.

Tiessen, H., Moir, J.O., 1993. Characterization of available P by sequential extraction. In: Carter, M.R. (Ed.), Soil Sampling and Methods of Analysis. Lewis, Boca Raton, Florida, pp. 75-86.

Treydte, K., Boda, S., Pannatier, E.G., Fonti, P., Frank, D., Ullrich, B., Saurer, M., Siegwolf, R., Battipaglia, G., Werner, W., Gessler, A., 2014. Seasonal transfer of oxygen isotopes from precipitation and soil to the tree ring: source water versus needle water enrichment. New Phytol. 202 (3), 772-783.

Watson, C.J., Matthews, D.I., 2008. A 10-year study of phosphorus balances and the impact of grazed grassland on total P redistribution within the soil profile. Eur. J. Soil Sci. 59 (6), 1171-1176.

Young, M.B., McLaughlin, K., Kendall, C., Stringfellow, W., Rollog, M., Elsbury, K., Donald, E., Paytan, A., 2009. Characterizing the oxygen isotopic composition of phosphate sources to aquatic ecosystems. Environ. Sci. Technol. 43 (14), 5190-5196.

Zohar, I., Shaviv, A., Klass, T., Roberts, K., Paytan, A., 2010. Method for the analysis of oxygen isotopic composition of soil phosphate fractions. Environ. Sci. Technol. 44 (19), 7583-7588.