Spatially Related Sampling Uncertainty in the Assessment of Labile Soil Carbon and Nitrogen in an Irish Forest Plantation

: The importance of labile soil carbon (C) and nitrogen (N) in soil biogeochemical processes is now well recognized. However, the quantiﬁcation of labile soil C and N in soils and the assessment of their contribution to ecosystem C and N budgets is often constrained by limited information on spatial variability. To address this, we examined spatial variability in dissolved organic carbon (DOC) and dissolved total nitrogen (DTN) in a Sitka spruce forest in central Ireland. The results showed moderate variations in the concentrations of DOC and DTN based on the mean, minimum, and maximum, as well as the coefﬁcients of variation. Residual values of DOC and DTN were shown to have moderate spatial autocorrelations, and the nugget sill ratios were 0.09% and 0.10%, respectively. Distribution maps revealed that both DOC and DTN concentrations in the study area decreased from the southeast. The variability of both DOC and DTN increased as the sampling area expanded and could be well parameterized as a power function of the sampling area. The cokriging technique performed better than the ordinary kriging for predictions of DOC and DTN, which are highly correlated. This study provides a statistically based assessment of spatial variations in DOC and DTN and identiﬁes the sampling effort required for their accurate quantiﬁcation, leading to improved assessments of forest ecosystem C and N budgets.


Introduction
Soil organic carbon (SOC) and total nitrogen (STN) play an important role in the global biogeochemical cycling of C and N within and between terrestrial ecosystems and between atmospheric and terrestrial pools [1][2][3].Whilst measurements of SOC and STN in soils are relatively straightforward, their quantification, and therefore their significance, is often hampered by marked spatial variability [4].This includes variability in both the vertical (soil depth) and horizontal directions [5].An understanding of this variability is clearly important for an improved understanding of the cycling of C and N and for increasing the precision of estimates of C and N stocks.
In soils, labile carbon and nitrogen originate from dissolved organic matter (DOM), which is defined as the entire pool of water-soluble organic matter that is either absorbed on soil or sediment particles or dissolved in interstitial pore water [6,7].Practically, it is often considered to be any organic matter fraction in solution that passes through a 0.45-µm filter [8,9].Water-extractable organic matter (WEOM), extracted from soil under various laboratory conditions, is widely considered to be a substitute for in situ DOM [10,11].Dissolved organic matter comprises only a small part of the total soil organic matter; nevertheless, the dissolved fraction is generally a reflection of the total soil organic matter, because the soluble phase tends to be in equilibrium with the solid phase [11].As soil DOM is significantly correlated with CO 2 and N 2 O emissions it can be a sensitive indicator of carbon and nitrogen mineralisation, and reflects environmental perturbations [12,13], particularly through its impact on the eutrophication of water bodies [14].
Various factors can control the dynamics of DOM, among which soil temperature and soil moisture are considered to be the two most important natural environmental drivers, often explaining trends in soil surface DOM concentrations [15].Increasing soil temperatures can also stimulate the decomposition of soil organic matter (SOM) through enhanced microbial activity [16].In general, the relationships between DOM concentrations and soil moisture vary considerably among different studies [15].Soil pH has often been shown to have a positive impact on DOM release from surface soils in laboratory studies, while these effects seem to be small in the field [15].All of the factors controlling DOM can vary with topographic position (e.g., a higher content of organic matter and water content in the lower slope positions compared to the upper slopes due to water-related movement) [17].
Due to the difficulty of soil sampling and associated costs, most studies evaluating dissolved organic matter rely on information obtained from relatively few representative sampling points or profiles [18][19][20].Given the known spatial variability in soil properties [21,22], there is often considerable uncertainty about the representativeness of this sampling approach and how the values relate to actual stocks or concentrations.Clearly, more intensive assessments are required in order to obtain a better understanding of the variability in dissolved organic matter, together with an examination of the reason(s) for any spatial variation.Quantification of ecosystem C and N stocks and their changes over time is a major challenge and requires the identification of an appropriate sampling strategy that accurately accounts for the often-large variability in soil properties.Assessment of the impact of sampling scale on the spatial variation of soil carbon and nitrogen and their relationship with environmental variables can also contribute to an improved understanding of causally related variations.The development of more robust sampling strategies will also improve the inputs to biogeochemical process-based models that are used for simulating soil carbon and nitrogen dynamics.
One of the measures that could be used to offset greenhouse gas (GHG) emissions by creating or enhancing C sinks is through afforestation and reforestation activities.Up to about 0.74 million hectares in Ireland have been forested on previously unmanaged grasslands, among which Sitka spruce is the primary species, representing 53% of total coniferous plantations [23].The majority of Stika spruce forest plantations in Ireland are associated with wet mineral gley soils.However, significant losses of carbon and nitrogen as dissolved forms from this and other ecosystems [20,[24][25][26]] could reduce their mitigation potential.Assessment of this will, however, depend on our ability to accurately capture the extent of spatial variations in soil C and N pools.In order to address this, the objectives of this study were: (1) to analyse spatial patterns in labile soil carbon and nitrogen and (2) to analyse the reason(s) for any variability.

Site Description
The study site was a Sitka spruce forest planted in 1988 in the Irish midlands (52 • 57 N, 7 • 15 W; altitude of 260 m), which is situated on the Castlecomer plateau.The mean annual temperature and precipitation values averaged over 30 years (1985-2015) are 9.3 • C and 850 mm, respectively.The studied area is characterized by a temperate climate typical of the northern temperate zone.The establishment of the forest and its management is under the control of Coillte, a semi-state company involved in forest-based land management.The Castlecomer plateau forms a distinctive geographical area not only because of its elevation and its plateau-like features, but also because it consists mainly of wet, rush-infested, agriculturally marginal land.The wetness of the land is due to the impermeable nature of a dense layer of boulder clay overlying the bedrock formations, resulting in the formation of wet gley soils.Based on the Irish soil information system (SIS) great soil group classification they are considered to be surface water gleys, which are categorized as stagnasols under the World Reference Base for Soil Resources (WRB) system.The soil is very shallow and subjected to periodic waterlogging.It is underlain by largely impermeable parent material.
The forest covers a total area of ~40 ha and is made up of two management compartments which are 25.8 ha and 16.2 ha in size, respectively.Ditches have been created to improve soil drainage and the trees were planted on a 2 m by 2 m grid at a planting density of 2300 stems ha −1 .The forest stands reach the canopy closure and are characterized by an almost total absence of understorey vegetation, except for patches of bryophytes which are mainly in proximity to the thinning lanes.The forest stand has been thinned in recent years [27] but never fertilized.The labile soil carbon and nitrogen investigations were conducted in the 25.8-ha stand, where an eddy covariance tower and associated instrumentation were installed for the purpose of providing long-term measurements of C fluxes, biomass, and climatic data [28,29].The parent material should be similar across the whole forest plantation, and the sampled area has a gently sloping topography from a northwest to a southeast direction with a slope of approximately 3% and a slope length of ~360 m.Across the site the slope varies from 1%-10%.This information, together with a more detailed description of the study site, is included in the work of Black et al. [30] and Zou et al. [27].The main soil properties of the study site are presented in Table 1.

Soil Sampling and Analysis
Soil samples were collected from an area of ~54,000 m 2 during July 2013.Each sampling location covered a grid area of 20 m × 20 m, with each separated by 20 m.In total, 111 sampling locations were chosen (Figure 1).At each location, sub-samples of soil (~7 cm 3 each) were collected using a soil auger (3 cm in diameter) at a depth of 0-10 cm at three sampling points which were approximately 100 cm apart.Manual measurements of soil temperature and moisture were taken at the same time as soil sampling using a WET sensor (Delta-T Devices Ltd., Cambridge, UK).The three sub-samples obtained at each location were subsequently mixed to produce one sample representing that location.After the soil samples were taken back to the laboratory they were air-dried and ground (<0.25 mm) prior to further physical and chemical analysis.The soil pH was determined in a water extract (soil: water, 1:2.5) using a pH meter (Thermo Fisher Scientific Inc., Waltham, MI, USA).Soils were extracted with 0.01 M CaCl 2 at a ratio of 25 mL:10 g (volume to oven dry mass equivalent), with this extraction procedure lasting approximately 30 min.The extracts were shaken for 10 min using an end-to-end shaker and then filtered through Whatman#42 filter paper into plastic tubes [11].The extracts were stored at −20 • C prior to dissolved organic carbon (DOC) and dissolved total nitrogen (DTN) measurements using a TIC-TOC analyser (TOC-V Shimadzu Corporation, Tokyo, Japan).

Data Analysis
Descriptive statistics, including the minimum, maximum, mean, skewness, kurtosis, standard deviation (SD), and coefficients of variation (CVs), and a one-sample Shapiro-Wilk test were used to explore the distribution characteristics of the different variables.We used correlation analysis to examine the relationships between DOC or DTN and the selected abiotic regulating factors.We further determined the main controlling factors using stepwise linear regressions.All the above statistical analyses were conducted using the programme SAS software v9.3 (SAS Institute Inc., Cary, NC, USA).
To examine the influence of sampling scale on soil variability, a series of sampling points using different moving windows were re-sampled using all sampling points (n = 111) within the study area [19,31] The determination of sampling number 'n' for DOC and DTN, requiring their estimation at a given accuracy level k, is given by: where µ a is a critical µ value corresponding to a certain confidence limit P 0 , of say, 99% or 95%, and k can be set to 5, 10, and 15%, depending on the accuracy required.We also used a geostatistical method to investigate the spatial variability of DOC and DTN.Before the semi-variogram calculation, a preliminary semi-variogram surface analysis was performed to detect any zonal effect or trend in direction [32].The semi-variogram of DOC and DTN were calculated using the following model: where γ (h) is the semi-variance for the lag interval h.Z(x i + h) and Z(x i ) are observations at positions x i + h and x i , respectively, and N(h) is the number of data pairs separated by a lag distance, h.The semi-variogram was fitted and generalized using least squares regression analysis, an approach also used by Zhang and Shao [19] in their study on soil organic carbon variability in the Gobi Desert of north-western China.The semi-variogram analysis was performed with GS+ 9.0 (Gamma Design Software, Plainwell, MI, USA).Ordinary kriging and cokriging techniques were used to estimate the values of DOC and DTN concentrations for unknown locations.When mapping DOC concentrations, DTN was set as an auxiliary variable, whilst DOC was set as an auxiliary variable when mapping DTN concentrations.The predictions of DOC and DTN by ordinary kriging and cokriging for the study area were evaluated based on the following parameters: mean error (ME), root mean square error (RMSE), and root mean square standardized error (RMSSE).Smaller ME and RMSE values indicate a model with fewer errors and more accurate predictions [33], whilst the RMSSE values closer to 1 indicate a model with better prediction [34].Finally, the spatial contour maps were produced with GIS software ArcGIS 10.1 and its extension module Geostatistical Analyst (ESRI Inc., Redlands, CA, USA).

Descriptive Statistics for DOC and DTN
Table 2 displays the main descriptive statistics for DOC and DTN in the study area.Mean DOC and DTN concentrations were 288.14 mg C kg −1 soil and 97.26 mg N kg −1 soil, respectively.Respective coefficients of variation (CVs) were 31.00% and 28.76% and were considered to indicate moderate variation.The variation in DOC and DTN increased with the expansion of the sampling area (Figure 2).The CVs for DOC and DTN ranged from 19.50% to 32.91% and 20.81% to 29.59%, respectively.The variations in soil temperature and pH were also shown to increase when the sampling area was expanded.The variation in soil moisture was much higher (CV ranges: 47%-52%) and was shown to have no clear relationship with the sampling area (Figure A1).

Spatial Variability of DOC and DTN
Since our data on DOC and DTN have no trend in any direction over the study area, an isotropic semi-variogram was created for DOC or DTN (Figure A2).The experimental semi-variograms showed that the best fit for the spatial variability of DOC and DTN was an exponential model (Figure A2).Table 3 lists all the essential parameters for the semi-variogram models.The nugget ratio of DOC (0.09) was similar to that of DTN (0.10).The nugget effects of DOC (740.00) and DTN (84.00) were relatively small compared to the respective sill values of 8140.00 and 812.10.The ranges of DOC and DTN were 34.2 and 35.1, respectively, and were rather similar.

Spatial Variability of DOC and DTN
Since our data on DOC and DTN have no trend in any direction over the study a an isotropic semi-variogram was created for DOC or DTN (Figure A2).The experime semi-variograms showed that the best fit for the spatial variability of DOC and DTN an exponential model (Figure A2).Table 3 lists all the essential parameters for the se variogram models.The nugget ratio of DOC (0.09) was similar to that of DTN (0.10).nugget effects of DOC (740.00) and DTN (84.00) were relatively small compared to respective sill values of 8140.00 and 812.10.The ranges of DOC and DTN were 34.2 35.1, respectively, and were rather similar.

Spatial Distribution Characteristics of DOC and DTN and Their Correlation with Environmental Factors
Figure 3 shows a comparison between the distribution maps of DOC and DTN generalized by ordinary kriging and those generalized by cokriging.Both methods generated similar spatial patterns of DOC and DTN, although most of the parameters (ME, RMSE, RMSSE) had smaller values using cokriging than those for ordinary kriging (Figure A4).Both DOC and DTN concentrations decreased from southeast to northwest across the study area.DOC and DTN were strongly correlated (Figure A3), with a correlation coefficient of 0.69 (p < 0.0001).However, it should be noted that there were still over 30% of DOC and DTN values that were not matched exactly, and this is why the spatial distribution maps for DOC and DTN were slightly different, but with a similar trend.The DOC and DTN densities were 23.92 and 8.07 g m −2 , respectively.
Figure 4 shows the correlation analyses of the relationships between environmental factors and DOM concentration, with DOC and DTN having similar patterns.Both DOC and DTN showed no significant correlations with environmental factors.However, elevation and soil temperature had a weak correlation (with p values falling between 0.05 and 0.1) with DOC and DTN, which indicates that the higher the elevation, the smaller the DOC and DTN concentrations, and the higher the soil temperature, the larger the DOC and DTN concentrations.
similar spatial patterns of DOC and DTN, although most of the parameters (ME, RMSE, RMSSE) had smaller values using cokriging than those for ordinary kriging (Figure A4).Both DOC and DTN concentrations decreased from southeast to northwest across the study area.DOC and DTN were strongly correlated (Figure A3), with a correlation coefficient of 0.69 (p < 0.0001).However, it should be noted that there were still over 30% of DOC and DTN values that were not matched exactly, and this is why the spatial distribution maps for DOC and DTN were slightly different, but with a similar trend.The DOC and DTN densities were 23.92 and 8.07 g m −2 , respectively.Figure 4 shows the correlation analyses of the relationships between environmental factors and DOM concentration, with DOC and DTN having similar patterns.Both DOC and DTN showed no significant correlations with environmental factors.However, elevation and soil temperature had a weak correlation (with p values falling between 0.05 and 0.1) with DOC and DTN, which indicates that the higher the elevation, the smaller the DOC and DTN concentrations, and the higher the soil temperature, the larger the DOC and DTN concentrations.

Discussion
In our study, the DOC and DTN concentrations fell well within the reported ranges of 2-400 mg kg −1 dry soil [17,35] and 1-292 mg kg −1 dry soil [35][36][37], respectively, of various forest soils.Moderate variability in dissolved carbon and nitrogen was observed in this study; however, the variability in DOC and DTN may be more complex and heterogeneous with an increase in the sampling area, because the factors regulating DOC and DTN

Discussion
In our study, the DOC and DTN concentrations fell well within the reported ranges of 2-400 mg kg −1 dry soil [17,35] and 1-292 mg kg −1 dry soil [35][36][37], respectively, of various forest soils.Moderate variability in dissolved carbon and nitrogen was observed in this study; however, the variability in DOC and DTN may be more complex and heterogeneous with an increase in the sampling area, because the factors regulating DOC and DTN variations change at different sampling scales [38].The value of the fractal power parameter for DOC (0.12) was slightly larger than that for DTN (0.087) indicating that the sampling scale may have a higher impact on the variability in DOC concentration.
A high spatial resolution over a large area is the ideal case for soil sampling [39] that, in practice, is hardly achievable.As the DOC and DTN data followed a normal distribution (p > 0.05), we could obtain a number for the samples required to estimate the mean soil DOC and DTN concentrations under different confidence levels and varying precision using Equation (1).The required number of samples was generally larger for higher confidence levels and greater precision.For instance, the sample number was generally larger for the 99% confidence level than for the 95% level, and further increased when the precision increased from 15% to 5% (Table A1).For the estimation of DOC and DTN, considerably more samples would have to be taken in a larger study area because of increased heterogeneity.Consequently, the sample number for labile C and N estimation with a precision of 5% for the current study was around 40, which is much smaller than the sample number that was actually used (111 locations).Details on the sample numbers required for estimating DOC and DTN at a given accuracy level with different sampling areas in this study can be found in Table A1.Therefore, careful consideration of the appropriate sample number to be used in the field needs to be undertaken both for practical reasons as well as for accuracy.
A nugget sill ratio smaller than 25% indicates strong spatial dependence for a variable of interest [40].The nugget ratios for both DOC and DTN were smaller 0.25, indicating a strong spatial dependence.The variability range determines the spatial autocorrelation.When a variable is within the range values, it is spatially autocorrelated, and when it is outside of the range values, it is not.This information provides guidelines for designing appropriate sampling schemes [41].For DOC samples with a distance smaller than 34.2 they were spatially correlated, whilst those with a distance larger than 34.2 were not spatially correlated.For DTN samples with a distance smaller than 35.1 these were spatially correlated, and for those with a distance larger than 35.1 they were not spatially correlated.As the sampling grid interval for this study (20 m) was smaller than the minimum range of 34.2, this validates our spatial analysis by meeting the requirements of sampling intervals.
The mean error (ME) and root mean square error (RMSE) were relatively small, and the root mean square standardized errors (RMSSEs) were close to 1 (Figure 4), all of which indicated that the predictions of DOC and DTN by ordinary kriging and cokriging for the study area were reliable [33,34].Most of the parameters had smaller values using cokriging than those using ordinary kriging, which suggests that the accuracy of the cokriging technique is higher than that of the ordinary kriging technique.This may be due to the fact that DOC and DTN have similar variability structures, and this may help to enhance the prediction accuracy when one variable is set as an auxiliary variable.An investigation on the spatial variability of soil carbon stocks in the Urucu river basin, Central Amazon (Brazil) also reported a better performance using cokriging interpolation for mapping SOC stocks when the compound topographic index (CTI, which is inversely correlated with SOC stock) was incorporated as an auxiliary variable [42].
In this plantation, disturbances due to human or animal activities are rare and are largely influenced by thinning-related forest management practices.Thus, other factors such as parent material, climate, topography, etc., might play a more important role in the spatial characteristics of DOM distribution at this site [19,43,44].Microbial activity is considered to be an important factor that directly impacts the production and consumption of soil DOM, though we did not measure any associated variations in soil microbial activity [15].Temperature, soil moisture, and pH can all affect microbial activity and thus the dynamics of soil DOM production [15].In this study, pH was shown to have a negative relationship with elevation; however, no significant relationships were found between elevation and soil moisture and soil temperature.Surprisingly, these factors were not significantly correlated with DOC and DTN.However, the lower values found at the higher elevations might indicate that the small-scale variations in topography that influence the movement of DOC and DON in solution could contribute to the observed spatial variability.Nevertheless, the study area used in this investigation was relatively small and the climate, soil type, parent materials, and terrain are likely to be broadly similar among the locations examined.This would suggest that the variability in DOC and DTN shown in this study is likely to depend more on soil biological characteristics, including microbial activity, which require further investigation.Given that a substantial proportion of the Irish forest estate is associated with wet mineral gley soils, the results may have some general applicability to other similar sites in Ireland and elsewhere.However, it may be difficult to extrapolate these results to other forestry plantations with different tree species, soils, or climatic characteristics and this also requires further study.

Conclusions
The following conclusions can be drawn from our investigations: (1) There was a moderate variation in DOC and DTN concentrations in the predominantly wet gley mineral soil of the Sitka spruce plantation, based on statistical and geostatistical analyses.(2) The variations in DOC and DTN concentration are scale-dependent and rise as the sampling area increases.Based on the analysis in this study, 40 sample locations would be required for adequately assessing labile C and N in this or other similar study areas.(3) Geostatistical analysis proved to be a powerful interpolation method for analysis of soil spatial variability and was useful for estimating DOC and DTN concentrations.(4) Cokriging techniques may have an advantage over ordinary kriging for predictions of highly corelated variables such as DOC and DTN in that they have higher prediction accuracy.
Appendix A

Figure 1 .
Figure 1.The site at the Dooary forest in central Ireland (left) where the soil samples were collected and the soil sampling locations within the study area (right).
. The chosen window sizes of 20 m × 20 m, 40 m × 40 m, 60 m × 60 m, 80 m × 80 m, 100 m × 100 m, 120 m × 120 m, and 140 m × 140 m within the initial sampling area were used to re-sample all the original sampling locations in a southeast to northwest direction.The metric CV (a mean CV value was calculated for each of the resampled areas listed below) was used to quantify the variability in DOC and DTN and other soil properties (pH, soil moisture, and soil temperature).

Figure 2 .
Figure 2. Relationship between the CV of DOC and DTN with increasing sampling area.Error bars represent the standard error of the mean.

Figure 2 .
Figure 2. Relationship between the CV of DOC and DTN with increasing sampling area.Error bars represent the standard error of the mean.

Figure 3 .
Figure 3. Spatial distribution of DOC and DTN estimated by ordinary kriging (a, b) and cokriging (c, d) across the study area.

Figure 3 .
Figure 3. Spatial distribution of DOC and DTN estimated by ordinary kriging (a,b) and cokriging (c,d) across the study area.

Figure 4 .
Figure 4. Correlation matrix (with individual p values shown in the upper left of the panel) of DOC and DTN with environmental factors.VW: soil moisture; T: soil temperature.

Figure 4 .
Figure 4. Correlation matrix (with individual p values shown in the upper left of the panel) of DOC and DTN with environmental factors.VW: soil moisture; T: soil temperature.

Figure A1 .
Figure A1.Relationship between the CV of pH, soil temperature, and moisture with increasing sampling areas.Error bars represent the standard error of the mean.

Figure A1 .Figure A2 .
Figure A1.Relationship between the CV of pH, soil temperature, and moisture with increasing sampling areas.Error bars represent the standard error of the mean.

Figure A3 .Figure A3 .
Figure A3.Correlation between DOC and DTN in the study area.

14 Figure A1 .Figure A2 .
Figure A1.Relationship between the CV of pH, soil temperature, and moisture with increasing sampling areas.Error bars represent the standard error of the mean.

Figure A3 .Figure A4 .
Figure A3.Correlation between DOC and DTN in the study area.

Table 1 .
Bulk soil properties at the study site.

Table 2 .
Descriptive statistics of dissolved organic carbon (DOC) and dissolved total nitrogen (DTN).SD: standard deviation; CV: coefficient of variation.

Table 3 .
Parameters of semi-variogram models for DOC and DTN.Nugget: the undetectable measurement error, inherent variability, or variation within the minimum sampling distance.Sill: total variation.Nugget sill ratio: the degree of spatial variability.RSS: residual sum of squares.

Table 3 .
Parameters of semi-variogram models for DOC and DTN.Nugget: the undetectable measurement error, inherent variability, or variation within the minimum sampling distance.Sill: total variation.Nugget sill ratio: the degree of spatial variability.RSS: residual sum of squares.

Table A1 .
Sample numbers (n) of DOC (DOC_n) and DTN (DTN_n) required for their estimation at a given accuracy level in different sampling areas.DOC_CV: coefficient of variation of DOC.DTN_CV: coefficient of variation of DTN.