Spatial variability in the amount of forest litter at the local scale in northeastern China: Kriging and cokriging approaches to interpolation

Abstract Litter is essential to promote nutrient cycling and maintain the sustainability of forest resources. However, its variability has not been sufficiently studied at the local scale. The prediction of litter amount using ordinary cokriging with Pearson correlation analysis (COKP) and ordinary cokriging with principal component analysis (COKPCA) was compared with that using ordinary kriging (OK) based on cross‐validation at the local scale of a 1‐ha plot over natural spruce–fir mixed forest in Jilin Province, China. Litter samples in semidecomposed (F) and complete decomposed (H) horizons were collected using an equidistant grid point sampling of 10 m × 10 m. Pearson correlation analysis and principal component analysis (PCA) were used to confirm auxiliary variables. The results showed that the amount of litter was 19.65 t/ha in the F horizon and 10.37 t/ha in the H horizon. The spatial structure variance ratio in the H horizon was smaller than that in the F horizon, indicative of its stronger spatial autocorrelation. Spatial distributions of litter amount in both horizons exhibited a patchy and heterogeneous pattern. Of the selected stand characteristics and litter properties, litter moisture content indicated the strongest relationship with litter amount. Cross‐validation revealed that COKPCA using the comprehensive score as an auxiliary variable produced the most accurate map. The average standard error and root‐mean‐square error between the predicted and measured values were always smaller, the mean error and mean standardized error were much closer to 0, and the root‐mean‐square standardized error was closer to 1 than COKP using litter moisture and OK. Therefore, a clear advantage of cokriging based on principal component analysis was observed and COKPCA was found to be a very useful approach for further interpolation prediction.

Spatial variation refers to the difference and diversity in different spatial locations within a certain region, affecting the distribution pattern of the community or landscape, and it is common but very important feature in forest ecosystems (Burt & Butcher, 2010;Li & Reynolds, 1995;Urban, O'Neill, & Shugart, 1987). Given that it is essential to understand biogeochemical cycles and other ecosystem functions, and because many biological and nonbiological processes that occur at a large scale are rooted in smaller-scale changes, there is a need to study the spatial heterogeneity of the amount of forest litter at the local scale.
Litter amount is linked to litter input and decomposition, as well as other factors such as topography, microclimate, vegetation composition, geological substrates, and soil physicochemical properties (Aceñolaza, Zamboni, Rodriguez, & Gallardo, 2010;Hermansah, Masunaga, Wakatsuki, & Aflizar, 2002;Kitayama & Aiba, 2002;Proctor, Nderson, Ogden, & Vallack, 1983;Röderstein, Hertel, & Leuschner, 2005;Staelens et al., 2011). There are particular controlling factors that are important in different research areas or at different scales. For example, similar latitude, altitude, slope, temperature, and precipitation have relatively weak effects on the spatial variability of the amount of litter at the local scale, and thus particular attention has been paid to the relationships between litter amount and stand characteristics and the properties of litter itself.
Previous studies dealt with the effects of stand characteristics on litter amount and decomposition (Cseresnyés, Csontos, & Bózsing, 2006;Jonard, 2008;Prescott, 2002;Sariyildiz, Anderson, & Kucuk, 2005;Trogisch, He, Hector, & Scherer-Lorenzen, 2016). These studies suggest that modifying stand conditions could result in complex influences on the amount of litter. Litter amount can affect the rates of decomposition by changing physicochemical characteristics and resource variability in litter horizons (Fang, Zhao, Zhou, Huang, & Liu, 2015;Gripp et al., 2018;Parsons et al., 2014). Conversely, the discrepancies in litter properties at different decomposition stages may also result in variability of litter amount. Recently, Petraglia et al. (2018) reported that litter quality not only shows the strongest relationship with early litter decomposition dynamics, but it also regulates the impact of water availability on mass loss. Owing to all the correlative factors, considerable uncertainties remain in the prediction of spatial variability in litter amount. This is notably true for mixed stands and natural forests, which have higher biodiversity and variability in canopy composition, and a higher rate of localized disturbances (Nickmans, Jonard, Verheyen, & Ponette, 2019;Wang, Wang, & Huang, 2008).
Geostatistics is a typical method that is used for investigating spatial patterns and interpolations based on measured values and distance between detecting points (Liu et al., 2014). Ordinary kriging (OK) is most often used for unbiased optimal interpolation of regionalized variables (Ferguson, Newbury, Maxted, Ford-Lloyd, & Robertson, 1998). Its capability to account for spatial dependence directly relies on the quantity and quality of the sample data (Miller, Franklin, & Aspinall, 2007), but its predictive ability is limited in regions with complicated environmental factors. Ordinary cokriging (COK) makes use of the correlation between different regionalized variables (Odeh, McBratney, & Chittleborough, 1995) and is an extension of kriging for the optimal estimation of regionalized variables developing from a single attribute to two or more coregionalized attributes (Goovaerts, 1998). The initial comparison between kriging and cokriging was based on theoretical analysis that indicated the superiority of cokriging and that it ensured consistency between an estimation of a sum and the separate estimation of each of its terms (Wackernagel, 1994). Recently, the enhancement of cokriging for spatial interpolation has been investigated in practical application (Adhikary, Muttil, & Yilmaz, 2017;Basaran et al., 2011;Kuntz & Helbich, 2014;Yang, Luo, Jiang, Li, & Yuan, 2016). Basaran et al. (2011) concluded that cokriging was superior to kriging in predicting soil hydraulic conductivity with limited available data on a 12.8 km 2 area. In their study, the water-stable aggregate was significantly related to soil hydraulic conductivity, and thus it was used as an auxiliary variable for the estimation of soil hydraulic conductivity. Yang et al. (2016) found that cokriging improved the result by 45.5% compared with kriging, with soil bulk density as a primary variable and soil water content as an auxiliary variable. Moreover, a few studies combined geostatistics with principal component analysis (PCA) in order to extract auxiliary variables (Borůvka, Mládková, Penížek, Drábek, & Vašát, 2007;Levi & Rasmussen, 2014;Markhvida, Ceferino, & Baker, 2018). Borůvka et al. (2007) applied the combination of geostatistics and PCA to study how stand factors affect the spatial distribution of soil characteristics. All these interpolation methods have mainly been applied to the studies of soil properties, but how to assess litter amount using cross-validation and select comprehensive auxiliary variables is not fully understood.
Our study was based on a 100 m × 100 m permanent plot over natural spruce-fir mixed forest that was established for litter sample collection in the semidecomposed (F) and complete decomposed (H) horizons and the survey of stand characteristics before the peak of fall (late August) in Jilin Province, China. Spatial variability of litter amount was evaluated using geostatistical analysis, and three spatial interpolation methods, including OK, COK P (using the strongest correlation with litter amount as an auxiliary variable), and COK PCA (using the comprehensive score as an auxiliary variable), were compared. The present study attempted to (a) examine spatial variability of litter amount in the F and H horizons at the local scale, (b) test the relationships of litter amount with selected factors of stand characteristics and litter properties, (c) investigate the possibility of using different auxiliary variables to predict litter amount, and (d) compare kriging with cokriging in litter amount estimation at the local scale based on cross-validation.

| Study site and sample plot
This study was conducted at the Jingouling Forest Farm, Wangqing County, Jilin Province (43°17′-43°25′N, 130°05′-130°20′E, an altitude of 300-1,200 m, a slope of 5-25°, see Figure 1). The climate is influenced by north temperature monsoon with a mean annual temperature of 3.9°C and an average annual precipitation ranging from 600 to 700 mm. Soils are developed over basalt, gneiss, or granite and are predominantly Bor-Udic Cambisols (Gong, 1999), which is equivalent to Humic Cambisols (IUSS Working Group WRB, 2014). The forest types that are present include coniferous forest, broad-leaved forest, and mixed coniferous and broad-leaved forest.
In July 2012, we established a 100 m × 100 m plot ( Figure 1

| Sample collection and stand survey
In August 2017, within the center of each 10 m × 10 m subplot, litter was collected in a quadrat of 50 cm × 50 cm, and it was separated into F and H horizons as there was little or even no litter in the undecomposed (L) horizon at this time. The F horizon consists of medium to strongly fragmented material with many mycelia and thin roots, whereas the H horizon is made up of a humified amorphous material (Papamichos, 1990). Subsequently, digital images of the canopy were obtained using Nikon fish-eye lenses. Tree species and DBH for the individual trees and stem number for each subplot were recorded if DBH was larger than 5 cm.

| Determination of litter properties and calculation of stand characteristics
The determination of litter amount, moisture content, organic carbon (OC), total nitrogen (TN) and total phosphorous (TP) concentrations, and calculation of stem number, species number, the proportion of conifer stems and species, Gleason (D), Pielous (J) and Shannon Wiener (H′) indices, and basal area were conducted following descriptions by Qin et al. (2019). In addition, canopy density was analyzed using the digital images of the canopy, which were processed using Adobe Photoshop (version CS6, Adobe Systems Software Ireland Ltd) and calculated as follows (Qi, Luo, & Zhao, 2009): where, ε is the canopy density, e and E are the pixel value of sky portion in the selected area and that of the selected area, respectively.
The selected area is encompassed by the fish-eye lense images that we processed to remove forest land, undergrowth bush and trees outside the plot.

| Descriptive statistical analyses
The mean, maximum, minimum, and standard deviation (SD) were calculated for the items of each horizon. The Kolmogorov-Smirnov statistical test, together with skewness and kurtosis values, was performed to check whether the data were normally distributed, and if not, logarithm transformation was used (Limpert & Stahel, 2011).
Pearson correlation analysis was used to test the correlations of litter amount with affecting factors. The PCA is a statistical analysis method that converts multiple indices into a few principal components (PCs) to reduce dimensionality, and it is calculated using the eigenvalue decomposition of a data covariance (or correlation) matrix or singular value decomposition of a data matrix after normalization of the original data (Armenise, Redmile-Gordon, Stellacci, Ciccarese, & Rubino, 2013;Everitt & Dunn, 1992). The cumulative contribution rates of the first several PCs > 80% are considered, and the PCA results are usually analyzed according to component and comprehensive scores as Equations 2 and 3 (Abdi & Williams, 2010;Shaw, 2009).
where, Y i is the component score of the selected PCs, j is the normalized data of factors related to litter amount, u ij is the factor loading, and l is the number of factors.
where, W is the comprehensive score, v k is the eigenvalue of PCs, and m is the number of the selected PCs.
All the descriptive parameters, Pearson correlation analysis, and PCA were obtained using SPSS (version 21.0, International Business Machines Corporation) and R (version 3.3.1, https ://www.r-proje ct.org/) for Windows.

| Geostatistical analysis
The spatial variability of litter amount was conducted using ArcGIS (version 10.2, http://devel opers.arcgis.com/). Based on a semivariogram and a cross-variogram, appropriate models were fitted and the input parameters for the spatial interpolation of kriging and cokriging were provided (Cressie, 1993;Goovaerts, 1997). Three parameters of nugget (C 0 ), sill (C + C 0 ), and range (A) can show the spatial dependence of litter amount (Burgos, Madejón, Pérez-de-Mora, & Cabrera, 2006), and the degree of spatial autocorrelation can be evaluated using structure variance ratio, that is, the value of C 0 /(C 0 + C) (Jiang et al., 2017).
In a spatial context, the kriging method can use the data in the neighborhood to predict the value of the litter amount at a location where it has not been measured (Wackernagel, 1994). The ordinary kriging estimates were calculated as Equation 4: where, Z OK (x 0 ) is the predicted value at the location of x 0 ; Z(x i ) is the measured value at the sampling site x i , λ i is a weighting coefficient, which can be calculated according to the unbiasedness, optimality and the Lagrange minimization principle (Sakata, Ashida, & Zako, 2004); and a is the number of sites within the search neighborhood around x 0 used for the prediction.
Cokriging is an interpolation method used to combine the litter amount and auxiliary variable while accounting for the spatial pattern and sampling procedure (Goovaerts, 2019). The ordinary cokriging estimates were calculated as Equation 5: where, Z COK (x 0 ) is the predicted value at the location of x 0 ; Z 1 (x j ), target variable, and Z 2 (x k ), auxiliary variable, are the measured values at the sampling site x j and x k ; λ 1j and λ 2k are weighting coefficients, and b and c are the number of Z 1 (x j ) and Z 2 (x k ) within the search neighborhood around x 0 used for the prediction.

| Accuracy evaluation
The accuracy of the spatial distribution map was checked using the cross-validation approach. Five indicators, including mean error (ME), mean standardized error (MSE), average standard error (ASE), root-mean-square error (RMSE), and root-meansquare standardized error (RMSSE), were used to evaluate the performance of ordinary kriging and cokriging predictions and calculated as follows (Ali & Othman, 2017;Johnston, Hoef, Krivoruchko, & Lucas, 2001;Zapatarios, Rivero, Naja, & Goovaert, 2012): and n is the number of data pairs. The ME provides a bias of predicted value, the ASE and RMSE provide the accuracy between the predicted value and measured value, and the MSE and RMSSE provide the accuracy of the standard error.

| Prediction procedure
Prediction of the litter amount was carried out following five steps: (1) The original datasets were processed using a normal distribution test and logarithm transformation if necessary; (2)   OC:TP, and OC:TP concentrations exhibited a reverse trend of being lower by 9.50%, 0.02%, 6.15%, 7.22%, 18.29%, and 9.92%, respectively, in the F horizon than in the H horizon. The CVs of stand characteristics were within the range of 6.17%-56.72%, indicating a small or moderate spatial discretization of stand characteristics in natural spruce-fir mixed forest (Table 2).

| Pearson correlation analysis
The correlation analysis exhibited some differences between the two decomposed horizons (Figure 2). In the F horizon, litter amount had positive correlations with tree species number (r = .237, p < .05) and Gleason index (r = .237, p < .05), but it had a negative correlation with moisture content (r = −.252, p < .05). In

| Principal component analysis
n and the increasing positive effect of TP (0.486) on PC1 was shown.

| Spatial distribution pattern
Based on the semivariogram concept and models, OK, COK P , and COK PCA interpolation methods were evaluated for the spatial distribution of litter amount in the two decomposed horizons ( Figure 3). Litter amount in the study area exhibited a patchy and heterogeneous distribution for both horizons. In the F horizon, litter amount showed several high-value centers but was discrete ranging between 29.99 and 54.77 t/ha. The spatial distribution differed along the south-north direction, but not for the west-east direction. In the H horizon, litter amount in the northwestern part had lower values of 3.86-7.94 t/ha, and they were mainly concentrated on the middle and southeastern aspect in the range of 11.31-24.86 t/ha. In both horizons, the spatial prediction map of litter amount using COK P and COK PCA showed that the distribution trend of litter amount was similar to that estimated using OK, but it was less spatially detailed or more uniform. Furthermore, COK PCA interpolation had the most spatial heterogeneity, revealing the randomness and dependence of litter amount spatial distribution in certain local areas.

| Estimation accuracy comparison
The estimation accuracy was listed for OK, COK P , and COK PCA to emphasize differences among the three datasets (Table 6). An ME close to zero indicates a lack of bias, and an MSE close to zero and an RMSSE close to one indicate smaller estimation errors, while the ASE and RMSE should be as small as possible. In the F horizon, where minimum ASE was 8.681, higher accuracy between the predicted value and the measured value was found using COK PCA . The RMSE decreased from 9.719 using OK to 9.439 and 9.200 using COK P and COK PCA , respectively.
For COK PCA , the ME and MSE were close to zero and RMSSE was close to one compared with the OK and COK P datasets. The five indicators also showed better predictions using COK PCA for the H horizon. In both horizons, the estimation accuracy showed COK PCA > COK P >OK.

F I G U R E 3
Distribution of litter amount based on the three spatial interpolation methods (n = 100). OK, ordinary kriging; COK P , ordinary cokriging with litter moisture content as an auxiliary variable; COK PCA , ordinary cokriging with the comprehensive score as an auxiliary variable TA B L E 6 Accuracy comparison of prediction (n = 100)   (Kavvadias, Alifragis, Tsiontsis, Brofas, & Stamatelos, 2001) is lower than our data. This could be attributed to the fact that spruce and fir with higher fiber and secondary metabolite concentrations are not easily decomposed (Hendricks & Boring, 1992;Monk & Gabrieison, 1985). The litter amount in our study was higher than that in similar spruce-fir mixed forest in Liaoning Province, China (13.52 t/ha) reported by Liu et al. (2017). The study of spatial heterogeneity requires a higher density of sampling points (100 subplots per hectare).

| D ISCUSS I ON
However, some studies have been conducted with fewer points or have been carried out at larger scales to calculate average values instead. Moreover, it is of great importance to focus on uniform standards for research methods.
The factors of stand characteristics and litter properties that were related to litter amount varied between the F and H horizons. The litter amount in the F horizon was significantly related to moisture content, species number, and Gleason index. This suggests that litter quality may be different between tree species (Staelens et al., 2011). We found that litter amount in the H horizon was significantly correlated with moisture content, TN, TP, OC:TP, TN:TP, and canopy density. The litter TN, TP, and TN:TP can affect the rate of litter decay and the leaching and release of N (Berg & McClaugherty, 2003). However, the finding that litter amount in the H horizon was significantly and negatively related to TN is in contrast with Scott and Binkley (1997). A positive or negative effect of site quality on litter amount in the F and H horizons beneath different tree species was noticed by Kavvadias et al. (2001). Canopy structure was likely to influence the important factors affecting the litter decomposition processes, such as the temperature and moisture content of the forest floor (Sariyildiz, 2008). Therefore, canopy density may affect the spatial variability of litter amount.
Our results indicated that litter amount had the highest correlation with litter moisture content in both horizons. Previous studies also found effects of litter moisture content on litter decomposition and animal activities (Halupa & Howes, 1995;Kavvadias et al., 2001;Levings & Windsor, 1984). Traditionally, much slower litter decomposition rates and fewer groups of arthropods in these moisture-limited forests resulted in litter accumulation.
The PCA summarized the relationships between litter amount and multiple factors. In our study, two PCs in the F horizon and three PCs in the H horizon were selected explaining more than 80% of total variance. In the F horizon, factor loadings (absolute values) were higher for stand characteristics (species number and Gleason index) in PC1, which was considered the stand component. The PC2 showed large positive loadings for litter moisture content and was considered the litter component. In the H horizon, litter properties in PC1 (TP, OC:TP and TN:TP) and PC2 (moisture content and TN) with the highest factor loading were considered litter components, and PC3 represented an individual contribution from canopy density and was considered the stand component. Obviously, stand characteristics had a higher contribution rate than litter properties in the F horizon, but an opposite result was found in the H horizon.
The different processes and effects in each horizon suggested that stand characteristics and litter properties mainly act on the F and H horizon, respectively (Andrews, Karlen, & Mitchell, 2002;Askari & Holden, 2015;Govaerts, Sayre, & Deckers, 2006). Before the peak of fall, litter in the F horizon was exposed at the soil surface, which is highly influenced by the light, precipitation, structure, and species composition of forests. On account of the H horizon beneath the F horizon, stand characteristics had relatively weaker effects on litter in this horizon.
According to the fitted spherical or exponential models, the spatial structure variance ratio in the H horizon was smaller than the F horizon, revealing that it had a stronger spatial autocorrelation. The spatial variability of litter amount may be influenced by many factors. Usually, a strong spatial autocorrelation of litter amount can be attributed to intrinsic factors, and a weak spatial autocorrelation can be attributed to extrinsic factors (Chang, Scrimshaw, Emmerson, & Lester, 1998 In the F horizon, litter amount showed several high-value centers and varied along the north-south direction. In the H horizon, lower litter amount was found on the northwestern aspect and higher litter amount was mainly concentrated in the middle and southeastern part. This may also be related to the spatial distribution of tree species and size in addition to the factors previously mentioned in our study. For example, Xia et al. (2015) noticed that different DBH of trees may result in differences in litter input and that the pattern of trees with larger DBH may play a more important role. The spatial patterns of litter amount imply that complex material dynamics exist in forest ecosystems.
In the F and H horizons, the estimated spatial distribution of litter amount using OK, COK P , and COK PCA methods was similar to each other but had different spatial details. Both COK P and COK PCA showed an improvement in the prediction of litter amount in comparison to OK. The COK prediction presented lower bias (ME), smaller prediction errors (MSE and RMSSE), and higher accuracy (ASE and RMSE), simultaneously. Our results agree with Chang et al. (1998) and Yang et al. (2016), who also found that cokriging was superior to kriging. This is mainly because the OK interpolation only considered the spatial information of litter amount, but the COK interpolation improved the estimation accuracy of litter amount by using an auxiliary variable. In previous studies, the spatial factor that had the strongest correlation with the target variable was taken as the auxiliary variable, and to some extent, the influence of other spatial factors on the target variable was ignored (Basaran et al., 2011). In the present study, the litter moisture content, which had the strongest correlation with litter amount, and the comprehensive scores based on PCA in the F and H horizons were used as auxiliary variables to carry out the COK interpolation, respectively. The degree of improvement using COK PCA prediction was obviously higher than that using COK P , revealing that COK PCA can reduce the loss of raw data and information, simplify the data structure and avoid subjective randomness to accurately predict and simulate the spatial pattern of litter amount in the study site. which used litter moisture content and comprehensive scores as auxiliary variables to predict the litter amount, generated a better map than the kriging. The application of COK PCA led to a more unbiased predicted value (ME closer to zero), more accuracy between the predicted value and measured value (smaller ASE and RMSE) and smaller prediction errors (MSE closer to zero and RMSSE closer to one) compared with the OK and COK P datasets. The combination of Pearson correlation analysis especially PCA with geostatistics offered brief and concise information about the spatial variability of litter amount and its relationship with stand characteristics and litter properties at a small scale. The feasibility and accuracy of using cokriging based on principal component analysis not only provided a scientific interpolation method for predicting the spatial distribution of litter amount but also valuable information for better enriching and developing geostatistics in forest ecosystems.

ACK N OWLED G EM ENTS
We thank Fuzeng WANG from Jingouling Forest Farm, Wangqing Forestry Bureau for providing assistance and support in the field work. This work was jointly supported by National Key R&D Program of China (Grant Nos. 2017YFC0504002 and 2017YFC0504101) and Special Fund for Forest Scientific Research in the Public Welfare (Grant No. 201504303).

CO N FLI C T O F I NTE R E S T
We declared that we have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data generated or analyzed during this study are included in this published article.