Phosphate stable oxygen isotope variability within a temperate agricultural

Inthisstudy,weconductaspatialanalysisofsoiltotalphosphorus(TP),acidextractablephosphate(PO 4 )andthe stable oxygen (O) isotope ratio within the PO 4 molecule ( δ 18 O PO 4 ) 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 (HCl PO 4 ) solution with the two variables exhibiting a strong positive correlation.The δ 18 O PO 4 of theextractedPO 4 ranged from17.0to 21.6 ‰ with a mean of18.8 ‰ (±0.8).While the spa-tialvariabilityofTotalP hasbeenresearchedatvariousscales,thisisthe ﬁ rststudy toassess thevariability ofsoil δ 18 O PO 4 at a ﬁ eld-scale resolution. We investigate whether or not δ 18 O PO 4 variability has any signi ﬁ cant relation-shipwith:(i)itselfwithrespecttospatialautocorrelationeffects;and(ii)HCl PO 4 ,elevationandslope-both globally and locally. Results indicate that δ 18 O PO 4 was not spatially autocorrelated; and that δ 18 O PO 4 was only weakly related to HCl PO 4 , elevation and slope, when considering the study ﬁ eld as a whole. Interestingly, the latter relationships appear to vary in strength locally. In particular, the δ 18 O PO 4 to HCl PO 4 relationship may depend on the underlying soil class and/or on different ﬁ eld managements that had operated across an historical north-south ﬁ eld division of the study ﬁ eld, a division that had been removed four years prior to this study.


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 μg of phosphate (PO 4 ) can cause eutrophication, but are an order of magnitude lower than soil PO 4 concentrations required for crop growth (Heathwaite and Dils, 2000).Identifying the different pollutant sources that are impacting on a water body is critical to 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 PO 4 molecule (δ 18 O PO 4 ) to isotopically characterise PO 4 sources and transformations.However, data on the δ 18 O PO 4 of different PO 4 sources remains limited (Tamburini et al., 2014;Young et al., 2009).
One potential source of PO 4 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 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 δ 18 O PO 4 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 δ 18 O PO 4 at a field-scale.We investigate whether or not δ 18 O PO 4 variability has any significant relationship with: (i) itself with respect to spatial autocorrelation effects; and (ii) soil P (extracted using 1 M HCl (HCl PO 4 )), elevation and slopeboth globally and locally within the study field.In particular, the following four null hypotheses are tested: A. δ 18 O PO 4 does not strongly co-vary with HCl PO 4 , with elevation and with slope across the field as a whole.B. δ 18 O PO 4 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 δ 18 O PO 4 and HCl PO 4 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.

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, referred 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 between 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/clover 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).

Sample design, collection and analysis
To quantify spatial variability in δ 18 O PO 4 , a sampling pattern was chosen with the view of a geostatistical analysis not only to δ 18 O PO 4 , 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 × 75 m grid); (ii) an intermediate scale (25 × 25 m grid); and (iii) a small scale (10 × 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 δ 18 O PO 4 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 PO 4 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 PO 4 in solution is then converted to silver phosphate (Ag 3 PO 4 ) though the addition of an Ag-ammine solution and subsequent adjustment of the pH to between 7 and 8 with 0.5 M HNO 3 before incubation for two days at 50 °C in an oven.
Soil PO 4 extracted using 1 M HCl (HCl PO 4 ) represents an integration of several potential PO 4 pools (Zohar et al., 2010): (i) the most labile, water-extractable soil PO 4 , including intracellular microbial PO 4 typically released through processes such as soil drying and re-wetting and cell lysis (ii) weakly adsorbed bicarbonate-extractable soil PO 4 , and (iii) strongly fixed, calcium and iron bound P. The weak acid extraction does not include the strongly adsorbed aluminium oxide bound PO 4 , 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 18 O-labeled and unlabelled 1 M HCl.If the extracted HCl PO 4 was to contain large amounts of hydrolised organic P or condensed PO 4 species (e.g.polyphosphates, pyrophosphates, etc.) then the δ 18 O PO 4 of the sample duplicates would be markedly different.For all tested samples, the 18 O of the sample extracted with unlabelled and 18 O-labeled acid was no N1%.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 b6.Therefore HCl PO 4 is assumed to extract PO 4 from the same soil pool as is released by water, anion resins, and NaHCO 3 extraction, namely extracellular labile PO 4 and metabolic intracellular microbial PO 4 and also some iron bound PO 4 .Using the 1 M HCl as an extractant enables far greater quantities PO 4 to become available, quantities that are required to allow the successful precipitation of Ag 3 PO 4 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 Ag 3 PO 4 were weighted into silver capsules together with a small amount of fine carbon black powder to promote the reaction between the Ag 3 PO 4 and carbon to produce CO.The produced reaction gases are carried through a temperature-controlled chromatography column and ultimately to the IRMS.Calibration and corrections for instrumental drifts were done by repeated measurements of an internal standard (ACROS Ag 3 PO 4 N 97.5%, δ 18 O = + 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 δ 18 O 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‰.

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 δ 18 O PO 4 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 δ 18 O PO 4 is taken as some function of HCl PO 4 , elevation, slope, soil class and the historical north-south division.In addition to these analyses, that are centred on δ 18 O PO 4 , complementary spatial analyses were conducted on TP and HCl PO 4 to provide useful context.All study hypothesis test results are reported at the 95% significance level; and all statistical analyses were implemented in R (http://www.r-project.org).

Standard exploratory analysis
Conditional boxplots and conditional scatterplots were used to assess δ 18 O PO 4 distributions and relationships.This included evidence for multiple populations for δ 18 O PO 4 , according to the historical northsouth field division or to soil class; and whether or not the relationship for δ 18 O PO 4 with HCl PO 4 , 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 Information 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.

Variographic analysis
To assess spatial dependence in δ 18 O PO 4 , 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 HCl PO 4 .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 twopoint 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).

Geographically weighted correlations
Spatial heterogeneity in δ 18 O PO 4 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 of bandwidth; (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.In 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).

Regression analysis
The following three regressions were considered, where δ 18 O PO 4 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 AIC procedure (Fotheringham et al., 2002).To compare the fit of chosen regressions, AIC and R 2 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 tvalues (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.

Distributions of TP, HCl PO 4 and δ 18 O PO 4 with the raw clustered data
The TP of the soil across Great Field ranged from 736 to 1952 mg P kg −1 , while the HCl PO 4 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 HCl PO 4 (r = 0.91), where the proportion of TP as HCl PO 4 was found to increase with increasing soil TP content.The δ 18 O PO 4 of the HCl PO 4 ranged from 17.0 to 21.6‰ with a mean of 18.8‰ (±0.8).
The spatial distributions of TP, HCl PO 4 and δ 18 O PO 4 are presented in Fig. 1a-c.The historic field divide appears an important discriminating variable in terms of both TP and HCl PO 4 with higher values to the north of the divide.Soil class also appears a useful discriminator of TP and HCl PO 4 .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 HCl PO 4 .Conversely, there does not appear to be any spatial trend in δ 18 O PO 4 or environmental factors that influence its variability.Note that all data descriptions in this section are naïve given that any bias due to the sample design, are not (as yet) considered.

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, HCl PO 4 and δ 18 O PO 4 , 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 HCl PO 4 ; and (iii) very slightly under-estimate it for δ 18 O PO 4 .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 δ 18 O PO 4 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.

Complementary analyses for TP and HCl PO 4
Although our focus is the spatial analysis of δ 18 O PO 4 , it is useful to provide an insight into how TP and HCl PO 4 vary spatially.This provides context to the analysis of δ 18 O PO 4 , especially as we choose to investigate its relationship to HCl PO 4 in subsequent sections.As opposed to TP, HCl PO 4 more directly relates to our understanding of δ 18 O PO 4 .Thus in Fig. 3, prediction surfaces are given for TP and HCl PO 4 ; both of which were constructed using a kriging with external drift (KED) model, 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 HCl PO 4 its trend was informed by elevation, slope, soil class and the north/south division.The R 2 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 HCl PO 4 were not used to help predict each other, as this would not inform predictions on a grid (δ 18 O PO 4 was similarly not used in this respect).Variography for both KED models was reasonably wellbehaved, 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 HCl PO 4 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 (Chilès and Delfiner, 1999).

Exploratory analysis for δ 18 O PO 4
Using both the clustered and strongly declustered data, the relationship matrices for δ 18 O PO 4 , HCl PO 4 , 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 δ 18 O PO 4 , the strongest relationship is with HCl PO 4 , but this is weak with a correlation of only 0.30 (for the weighted correlation).Clearly, δ 18 O PO 4 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 δ 18 O PO 4 with HCl PO 4 .
Conditional boxplots are used in Fig. 4c & d to relate the distribution of δ 18 O PO 4 to both the historic field division and to soil class; again using both clustered and strongly declustered data.Marginally higher δ 18 O PO 4 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 δ 18 O PO 4 .In general, there is little to choose between the clustered and declustered data analyses.Thus the sample design does not appear to strongly bias δ 18 O PO 4 in this respect.Regardless of any biasing effects, all relationships for δ 18 O PO 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 δ 18 O PO 4 as the response.Table 3 reports the results using all predictors (HCl PO 4 , 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 R 2 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 naïve 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 R 2 values result when directly using the clustered data).None of the predictor variables appear particularly worthy predictors of δ 18 O PO 4 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.

Variographic analysis for δ 18 O PO 4
We next investigate for spatial dependence in δ 18 O PO 4 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 δ 18 O PO 4 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 δ 18 O PO 4 , cross-dependence with HCl PO 4 was also absent (even though spatial dependence is present in  HCl PO 4 , from Section 3.3).Given these results, it is unnecessary to formally test for spatial dependence in δ 18 O PO 4 and study hypothesis (B) is accepted.
From Fig. 5, the effects of data clustering on the variograms are quite evident, where semi-variances at the lower lags of the clustered data variograms are commonly heightened, resulting in very poor small-  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 δ 18 O PO 4 has a higher variance in the north than that found in the south.The behaviour of the withinclass 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.

Local relationships for δ 18 O PO 4
Given that spatial autocorrelation effects are absent for the δ 18 O PO 4 process, out next objective is to assess for any spatial heterogenic effects with respect to the relationships between δ 18 O PO 4 with HCl PO 4 , 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 δ 18 O PO 4 and HCl PO 4 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 of these correlations (using standard t-tests with correlations) with relative assurance.Here the global correlation of 0.29 between HCl PO 4 and δ 18 O PO 4 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 spatiallyvarying 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 δ 18 O PO 4 with HCl PO 4 , with elevation, and with slope; each specified with a bandwidth of 40%.Co-variability for δ 18 O PO 4 with HCl PO 4 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 HCl PO 4 and δ 18 O PO 4 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 δ 18 O PO 4 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.

Regression analysis for δ 18 O PO 4
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 R 2 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 N0.05.Given these test results, study hypothesis (C) is accepted.
The localised regression coefficients from GWR associated with HCl PO 4 are mapped in Fig. 7a, where it appears that the relationship for δ 18 O PO 4 with HCl PO 4 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 HCl PO 4 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.

Variability of HCl PO 4 and TP
Analyses for HCl PO 4 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 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 HCl PO 4 indicate that N50% 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 PO 4 it is assumed that the bulk of the HCl PO 4 is comprised of PO 4 adsorbed to soil particles and microbial PO 4 .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 HCl PO 4 measured are very similar to the levels of microbial P reported by Blackwell et al. (2013) and therefore we assume that HCl PO 4 in these soils represents, in large part, adsorbed soil PO 4 and intracellular microbial PO 4 , the latter normally exceeding the former (Blackwell et al., 2010).

Variability of δ 18 O PO 4
In order to discuss δ 18 O PO 4 values observed, we need to estimate the theoretical equilibrium δ 18 O PO 4 value that might be expected for PO 4 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 PO 4 oxygen with the oxygen in H 2 O 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 PO 4 and water, catalyzed by pyrophosphatase, over a typical environmental temperatures (3-37 °C): Eδ 18 O PO 4 = − 0.18 T + 26.3 + δ 18 O H 2 O where Eδ 18 O PO 4 is the stable oxygen isotope ratio of PO 4 at equilibrium in ‰, T is the temperature in degrees Celsius and δ 18 O H 2 O is the stable oxygen isotope ratio of H 2 O in ‰.The intracellular phosphate, already at equilibrium, is released to the soil after cell lysis.
Although measurements of soil water δ 18 O H 2 O were not directly made it has been suggested that soil water at the location is very similar δ 18 O H 2 O to ground water, in that it is an integrated value of many rainfall events (Granger et al., 2010).As such the δ 18 O H 2 O should be similar to that of the global meteoric water line for this area with δ 18 O H 2 O 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 δ 18 O PO 4 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 δ 18 O PO 4 and temperature values to understand better the measured δ 18 O PO 4 data compared to the theoretical equilibrium δ 18 O PO 4 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 δ 18 O H 2 O is more pronounced and subject to strong influence from individual rainfall events and also from enrichment in 18 O 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 δ 18 O H 2 O is enriched by up to 4‰ compared to groundwater, then the theoretical δ 18 O PO 4 equilibrium values range between 18.0 and 23.0‰.This range of values better matches the δ 18 O PO 4 of the soil and suggests that the HCl PO 4 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 PO 4 and HCl PO 4 were broadly in equilibrium with soil water and proposed, amongst other reasons, that a flux of intracellular equilibrated PO 4 from the soil microbial biomass may be causing this.Tamburini et al. (2012) also conclude that in young and developing alpine soils, regardless of the contribution of PO 4 from minerals or vegetation, or from the activity of  extracellular enzymes, microbial biomass was the main controller of the δ 18 O PO 4 , 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 δ 18 O PO 4 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 δ 18 O PO 4 signatures similar to that of the equilibrium δ 18 O PO 4 value.The δ 18 O PO 4 of inorganic fertilizers have been reported to show a very wide range of values from 14.8 to 27.0‰ while virtually no δ 18 O PO 4 values have been published for animal excreta (Young et al., 2009) or for stored and managed animal wastes, however, water-extractable PO 4 from dairy farm slurries in this area have been found to range from 12.0 and 15.0‰ (Granger, unpublished data).Analysis results for δ 18 O PO 4 relationships within the data are enigmatic.Considering the field as a whole, all relationships for δ 18 O PO 4 are weak or indistinct.HCl PO 4 , elevation, slope, soil class and the historic field divide are not particularly good drivers of δ 18 O PO 4 variability; and δ 18 O PO 4 is not spatially correlated with itself.For within-field  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, δ 18 O PO 4 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 δ 18 O PO 4 with HCl PO 4 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 HCl PO 4 and δ 18 O PO 4 that is present in the north.The scenario described by the data seems to suggest that, in the north, elevated HCl PO 4 occurs as 'hotspots' with associated higher δ 18 O PO 4 , and that as the HCl PO 4 of these hotspots declines so the δ 18 O PO 4 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 et al., 2005;Penn et al., 2007).However, the available data on δ 18 O PO 4 in excreta, although limited, indicates that it is not elevated.If however metabolic PO 4 released from leaf litter has an elevated δ 18 O PO 4 (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 δ 18 O PO 4 and HCl PO 4 it is not reflected in the soil total carbon content which shows no relationship with δ 18 O PO 4 .

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 HCl PO 4 , TP and δ 18 O PO 4 .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 HCl PO 4 , TP or δ 18 O PO 4 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.

Analysis limitations
Although the statistical analyses for δ 18 O PO 4 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 δ 18 O PO 4 .In this respect, future δ 18 O PO 4 (and HCl PO 4 ) 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 collinearity (Páez 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 δ 18 O PO 4 may be driven by environmental covariates not considered in this study, as such, future studies should consider identifying these missing covariates.Measurement error in δ 18 O PO 4 could also be mitigating factor in this study's null results and could be incorporated in future work, noting that the δ 18 O PO 4 values of this study represent an average of a three replicates.

Conclusions
This study has shown that, despite differences in elevation, slope, the under lying soil class and management, the δ 18 O PO 4 has no relationship with these field variables even though the HCl PO 4 can have.The soil δ 18 O PO 4 values within the field were found to be within the range predicted for that has been microbially cycled and is at equilibrium with δ 18 O H 2 O .However given the lack of information of PO 4 source δ 18 O PO 4 signatures, it is not clear whether this is due to the rapid microbial cycling of PO 4 or because PO 4 sources have a similar δ 18 O PO 4 to that of the equilibrium value.Although no relationships were found between the δ 18 O PO 4 , and itself (with respect to spatial autocorrelation) and the field variables investigated, there was a suggestion of a positive correlation between δ 18 O PO 4 and HCl PO 4 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 PO 4 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.

Fig. 1 .
Fig. 1.The spatial distribution of (a) TP, (b) HCl PO 4 , and (c) δ 18 O PO 4 of the HCl PO 4 within 'Great Field' on the North Wyke Farm Platform.

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

Fig. 4 .
Fig. 4. Scatterplots and correlation coefficients for the (a) clustered and (b) strongly declustered data, respectively.These display the pairwise relationships between δ 18 O PO 4 , HCl PO 4 , 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 δ 18 O PO 4 , with respect to the historic field division or the three soil classes for (c) clustered and (d) strongly declustered data, respectively.

Fig. 6 .
Fig. 6.(a) Strongly declustered data conditional scatterplot for δ 18 O PO 4 with HCl PO 4 .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 δ 18 O PO 4 with (b) HCl PO 4 , (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.

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

Table 1
Global means of clustered and declustered data sets.

Table 2
Weighted correlations using the cell-declustering weights.

Table 3
MLR fits using the full predictor variable set.

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

Table 5
Bootstrap test q-statistics, associated p-values and randomisation test p-values. 18O PO 4 and HCl PO 4 , elevation and slope; only the local relationship for δ 18 O PO 4 with HCl PO 4 appears to have any value,