## 2 Materials and Methods
### 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 chosen window sizes of 20 m \(\times\) 20 m, 40 m \(\times\) 40 m, 60 m \(\times\) 60 m, 80 m \(\times\) 80 m, 100 m \(\times\) 100 m, 120 m \(\times\) 120 m, and 140 m \(\times\) 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).  
The determination of sampling number '\(n\)' for DOC and DTN, requiring their estimation at a given accuracy level k, is given by:  
\[n={\mu_{\text{a}}}^{2}{\left(\text{CV}\%/\text{k}\right)}^{2} \tag{1}\]  
where \(\mu_{\text{a}}\) is a critical \(\mu\) 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  
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**).  
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:  
\[\gamma(\mathrm{h})=\frac{1}{2\mathrm{N}(\mathrm{h})}\sum_{i=1}^{\mathrm{N}( \mathrm{h})}\left[Z(\mathrm{x}_{i})-Z(\mathrm{x}_{i}+\mathrm{h})\right]^{2} \tag{2}\]  
where \(\gamma\) (h) is the semi-variance for the lag interval h. \(Z(\mathrm{x}_{i}+\mathrm{h})\) and \(Z(\mathrm{x}_{i})\) are observations at positions \(\mathrm{x}_{i}+\mathrm{h}\) and \(\mathrm{x}_{i}\), respectively, and \(\mathrm{N}(\mathrm{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 [PERSON] and [PERSON] [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 D 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).