Spatial Variation and Stock Estimation of Soil Organic Carbon in Cropland in the Black Soil Region of Northeast China

Abstract

Soil organic carbon (SOC) sequestration in cropland is not only instrumental in combating climate change, but it also significantly enhances soil fertility. It is imperative to precisely and accurately quantify the SOC sequestration potential and assess the relative significance of various multiple explanatory factors in a timely manner. We studied 555 soil samples from the cropland topsoil (0–15 cm) across the black soil region in Northeast China between the years 2021 and 2022, and we identified 16 significant impact factors using one-way ANOVA and Pearson correlation coefficient analysis. In addition, the Random Forest (RF) model outperformed the Cubist model in predicting the spatial distribution of SOC contents. The predicted ranges of SOC contents span from 5.24 to 43.93 g/kg, with the average SOC content using the RF model standing at 17.24 g/kg in Northeast China. Stepwise regression and structural equation modeling revealed climate and topography as key factors affecting SOC distribution. The SOC density in the study area varied from 0.51 to 9.11 kg/m², averaging 3.30 kg/m², with a total SOC stock of 1226.64 Tg. The SOC sequestration potential in the study area was estimated at 3057.65 Tg by the categorical maximum method, with a remaining sequestration capacity of 1831.01 Tg. The study area has great potential for SOC sequestration. We hope to transform the theoretical value of SOC sequestration potential into actual SOC sequestration capacity by promoting sustainable agriculture and additional strategies. Our findings provide insights into the global soil conditions, SOC storage capacities, and effective SOC management strategies.

1. Introduction

Cropland is universally acknowledged as a significant carbon reservoir and plays a vital role in maintaining the global carbon equilibrium [1,2]. Soil organic carbon (SOC) in cropland accounts for more than 60% of the terrestrial carbon pool. The content and changes in cropland SOC directly indicate the persistence, stability, or loss of soil carbon [3,4,5]. It is estimated that carbon sequestration in global farmland has the potential to annually offset up to 1.2 Pg of carbon emissions, equivalent to 11.8% of global carbon emissions in 2023 [3,6,7]. Furthermore, SOC sequestration can effectively improve soil quality, thereby buffering the impacts of climate change on crop yields [8,9]. It is imperative to precisely quantify the sequestration potential of SOC and evaluate the comparative significance of various driving factors, particularly crop management practices.

In recent years, there has been a growing trend towards employing machine learning techniques for estimating SOC content or SOC stock, with the results exhibiting promising prospects, as evidenced by Odebiri et al. [10]. Their accuracy depends largely on the parameterization of these models by dedicated experiments. Different machine learning algorithms, especially variations of Random Forest (RF), are investigated for their ability to predict SOC pools. In terms of a machine learning-based RF model, incorporating 15 environmental predictor variables has been adopted for meta-analysis to assess the global response of SOC stocks to global warming at the resolution of 0.0083 °C [11]. Additionally, a myriad of studies has harnessed cutting-edge techniques, including remote sensing, to map SOC content at both landscape and regional scales [12,13]. These state-of-the-art methods have significantly contributed to our understanding of SOC distribution and variability [10].

SOC content levels in cropland is affected by both environmental factors and human activities; its fluctuations influence the trades-off between organic carbon inputs and carbon effluxes [14]. The landscape pattern indeed exhibits a ‘spillover effect’ on the SOC contents, as evidenced by studies conducted by Bowie et al. [15] and Didham et al. [16]. It is crucial to accord special attention to climate change and the management of crop cultivation. Temperature stimulates the SOC turnover by promoting microbial activities and vegetation-derived carbon inputs. Precipitation inhibits SOC decomposition by changing the soil anaerobic environment [14]. However, the escalating global climate change trends, encompassing increasing temperatures and dwindling precipitation, are unfavorable for the accumulation of SOC [11]. SOC dynamics are strongly affected by conservational tillage like straw residues returning and crop rotation [17,18,19,20].

Regrading research of SOC in vast spatial domains, typically, data-driven models, which are trained on publicly available literature and database, are employed to predict the SOC stock maps [6,21]. However, results are not fully representative due to inconsistencies in SOC measurement. Additionally, one significant factor contributing to this is the inherent difficulty associated with large-scale sampling. In this study, we collected a significant number of cropland soil samples from the black soil region in Northeast China in a short period of time, enabling data-driven SOC contents spatial distribution analysis.

To address the inquiries, we undertook a comprehensive evaluation of the factors influencing the spatial variability of SOC in croplands. Drawing from a comprehensive analysis of 555 soil samples, we took into account both agricultural management practices and landscape metrics. Consequently, we devised a predictive SOC model, leading to the formulation of a detailed spatial distribution map. Furthermore, utilizing the SOC stock results, we calculated the SOC reserves within the Northeast Black Soil Region, thereby assessing its potential for SOC sequestration. Ultimately, our objective is to foster responsible utilization and conservation of black soil, offering the insights for strategies that can augment croplands’ contribution towards China’s carbon neutrality goal.

2. Materials and Methods

2.1. Study Area and Sampling

2.1.1. Study Area

The studied area in Northeast China’s black soil region (115°31′ E–135°5′ E, 38°43′ N–53°33′ N), spanning 1.09 million km², sits at 600–1000 m above sea level (Figure 1). It features a basin-like terrain surrounded by mountains and flat ground in the center. The climate is a temperate monsoon climate with an annual temperature of −5 to 10.6 °C and precipitation of 400 to 1000 mm [22]. The cropland area within the black soil region of Northeast China spans a vast 3.59 × 10⁷ ha. Soil types include black soils, chernozems, albic soils, meadow soils, dark-brown earths, brown earths, and paddy soils [23] according to the Chinese Soil Taxonomy [24]. This black soil region in Northeast China holds a pivotal position as a significant commercial grain production base in the country.

2.1.2. Soil Sampling and Preconditioning

Soil samples were collected according to the requirements of soil environmental monitoring technical specifications (HJ/T 166-2004). All 555 soil samples, collected from a depth of 0–15 cm from the surface of cropland in the black soil region of Northeast China, were gathered during the months of July and August in 2021 and 2022. In this study, the grid distribution method was adopted in this study, and the factors such as land use type (paddy field, dry land), tillage mode, crop type and area in each grid were fully considered to select the points. The sampling locations were pinpointed using GPS and depicted as red dots on the map (Figure 1). The samples were air-dried in the laboratory and ground using a porcelain mortar, and each sample was ground in a mortar to traverse a 100-mesh polyethylene sieve. The SOC content was then measured according to the determination of soil organic carbon by potassium dichromate oxidation-spectrophotometric method (HJ 615-2011).

2.2. Environmental Variables

The selected environmental variables in this study include natural factors, cropping and management factors, and landscape metrics (Table 1). Then identify significant impact factors using the One-way ANOVA or Pearson correlation coefficient analysis.

2.2.1. Natural Factors

The spatial distribution of ST, TT, MAP, and MAT were downloaded from the Resource and Environmental Science and Data Center of the Chinese Academy of Sciences [27] with 1 km spatial resolution. The DEM data, featuring a spatial resolution of 30 m, were procured from the ASTER satellite elevation image, courtesy of the geospatial data cloud platform [28]. This data was further employed to determine the slope within ArcGIS Pro 3.0.2. The distance factor refers to the distance from each sampling point to the closest water body, as well as the green space (including forest and grassland) within the studied area. The Euclidean distance was used in ArcGIS Pro 3.0.2 to calculate Distowater and Distogreen.

2.2.2. Cropping and Management Factors

According to the statistics of the [22] released by the Chinese Academy of Sciences in 2022, the proportion of grain crops within the total cultivated area was up to 93.3% by the year 2020. The percentage of rice, corn, and soybean in the northeastern region to the sown area of food crops was as high as 98.9%. In addition, the spatial distribution of tillage mode was derived through the integration of crop distribution data spanning from 2017 to 2019, using ArcGIS Pro 3.0.2 [25].

The Modis image was obtained from the [29] platform, and it encompasses the timeframe spanning from September to October 2021. Subsequent to undergoing the preparatory steps of radiometric calibration and atmospheric correction within the ENVI 5.6, the NDVI was derived using the Red and NIR spectral bands. The computation of the NDVI index for each individual month within specified period has been successfully concluded.

(1)$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}={\displaystyle \frac{(\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d})}{(\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{e}\mathrm{d})}}$

The Red and NIR are the reflectance of Modis at 620–670 nm and 841–876 nm, respectively.

2.2.3. Landscape Metrics

Landscape metrics quantitatively reflect the structural composition and spatial configuration characteristics of the landscape [30]. In this context, six distinct types of landscape metrics were evaluated at landscape level and class levels (Table S1). Nine landscape metrics were calculated using Fragstats 4.0 from 30 m land use type raster data, incorporating 300 m landscape metrics into analysis.

2.3. Modelling Methodology

RF regression model relies on regression trees, cultivating an ensemble of such trees by sampling from a random subset of the initial training data [21]. The culmination of these trees’ predictions is achieved by averaging their respective outcomes, thus yielding the final prediction. Notably, the model’s parameters play a pivotal role in determining the accuracy of its predictions. The RF regression model encompasses two crucial parameters: the number of trees (ntree: 1–1000) and the number of nodes (mtry: 1 to the number of independent variables). It is necessary to conduct numerous trials to ascertain the optimal set of parameters. The %InMSE metric represents the extent to which the model’s mean absolute error (MAE) on the test set augments when each independent variable is excluded, thereby reflecting the significance of that particular variable.

The cubist model is a rule-based linear regression tree model that employs an ensemble learning strategy. It first establishes a set of rules by splitting each dependent variable and subsequently divides the data into distinct arranged in a tree-like formation. Subsequently, the model derives the linear regression for each of these subsets. The rules and regressions are formulated with the objective of minimizing the average absolute error for predicting new cases. The cubist’s hierarchical linear regression approach captures the nonlinear relationship between independent and dependent variables, while maintaining interpretability [31]. It is imperative to adjust the committee size (ranging from 1 to 100) and the number of neighbors (between 0 and 9) to enhance the predictive capacities of the model. The attribute usage of the conds and the attribute usage of the model indicates the significance of the variable.

The 555 samples were randomly divided into a train set (a total of 444) and a test set (a total of 111). The following indices were calculated to evaluate model performance: the coefficients of determination (R²), MAE, the root mean square error (RMSE), and Lin’s concordance correlation coefficient (LCCC). Generally, the model is deemed to exhibit superior performance with a high R² and LCCC value, coupled with low RMSE and MAE.

The equations are as follows:

(2)${\mathrm{R}}^2=1-{\displaystyle \frac{\sum _1^n{(O_i-P_i)}^2}{\sum _1^n{(O_i-\overline{O})}^2}}$

(3)$\mathrm{M}\mathrm{A}\mathrm{E}={\displaystyle \frac{1}{n}}{\sum }_1^n|O_i-P_i|$

(4)$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt[2]{{\displaystyle \frac{1}{n}}{\sum }_1^n{(O_i-P_i)}^2}$

(5)$\mathrm{L}\mathrm{C}\mathrm{C}\mathrm{C}={\displaystyle \frac{2r{S}_O{S}_P}{S_O^2+S_p^2+{\left(\overline{O}-\overline{P}\right)}^2}}$

where, n is the number of validation sample points, $O_i$ is the observed value at sample point i, $p_i$ is the predicted value at sample point i, $\overline{O}$ is the average of the observed value,$\overline{P}$ is the average of the predicted value, r is the Pearson correlation coefficient between the observed value and the predicted value, $S_O$ is the standard deviation of the observed value, and $S_P$ is the standard deviation of the predicted value.

2.4. Environmental Variable Analyses for SOC

2.4.1. Stepwise Regression Model

Stepwise regression modeling is a statistical technique that systematically selects independent variables step-by-step to build a regression equation. This process involves gradually incorporating or omitting variables to determine the optimal combination that exerts the most notable impact on the dependent variable, thereby establishing the most effective regression equation. In the modeling process, the stepwise regression model not only considers the correlations among independent variables but also adeptly manages the model’s complexity, therefore preventing overfitting issues. This model is widely employed in exploratory data analysis and predictive modeling, effectively assisting researchers in identifying the most pertinent variables from a vast array of independents, thereby enhancing the explanatory and predictive accuracy of the model. The significant influencing factors determined by One-way ANOVA and Pearson correlation coefficient analysis are used for the stepwise procedure. This process is typically executed using SPSS 27.

2.4.2. Structure Equation Modeling

Structural Equation Modeling (SEM) is an advanced statistical tool crafted specially to delve into intricate causal dynamics linking diverse variables [32]. This methodology integrates the prowess of factor analysis with multivariate regression analysis, empowering researchers to concurrently evaluate numerous interdependencies and tackle both observable and latent variables [33]. SEM establishes a linkage between latent and observable variable, delving into the interplay between the latent variables and causal pathways via structural modeling. The cornerstone of this methodology lies in its capacity to precisely quantify errors, enhance the model’s explanatory powers, and substantiate the theoretical framework through rigorous statistical analysis.

2.5. Estimation of SOC Stocks and Sequestration Potential

Soil Organic Carbon Density (SOCD) was calculated using the Schwager and Mikhailova formula [34], which was obtained from the SOC content (Equation (6)). The estimation of SOC stock in the upper 15 cm layer of cropland within the northeastern black soil zone is shown in Equation (7). Here, we used the categorical maximum method to estimate the soil organic carbon sequestration potential (SOCsp) in the study area. This method takes the SOC stock when all soils reach the maximum SOC content predicted by the corresponding soil type as the SOCsp. Soil organic carbon potential sequestration capacity (SOCpc) represents the difference between SOCsp and the current SOCstock, calculated using Equation (8):

(6)$\mathrm{S}\mathrm{O}\mathrm{C}\mathrm{D}=\mathrm{C}\times \mathrm{B}\mathrm{D}\times \mathrm{D}\times (1-{\displaystyle \frac{\mathrm{G}}{100}})$

(7)${SOC}_{stock}={\sum }_{i=1}^n{S}_i\times SOC{D}_i/{10}^9$

(8)${SOC}_{pc}={SOC}_{sp}-{SOC}_{stock}$

where C represents SOC content when saturated (g/kg); BD represents bulk density of the soil layer (g/cm³); D represents the thickness (cm); G represents gravel content of the soil sample (expressed as %); $S_i$ denotes the area of soil share/class $i$ (×10⁹ m²); SOC_stock refers to soil organic carbon stock (Tg).

3. Results and Discussion

3.1. Descriptive Statistic

A total of 444 samples were randomly selected to comprise the training set, while the remaining 20% of the sample, total 111, were used as the validation set. The SOC content of the sampling points ranged from 2.55 to 72.39 g/kg, with an average of 19.20 g/kg, which was higher than the median value of 17.34 g/kg, indicating that a substantial presence of sampling points with elevated SOC content (Table 2). In addition, the standard deviation (SD) was 10.67 g/kg, and the coefficient of variation amounted to 55.58%, indicating that a high degree of dispersion in the SOC content across the sampling points.

The spatial distribution of SOC content (Figure S1) exhibits a distinct trend of being lower in the southwest and higher in the northeast, albeit with marked spatial variability. This non-smooth spatial variability potentially stems from a multitude of factors, including disparities in soil properties at minute scales, fluctuations in topographic complexity, and diversity of land utilization patterns. Collectively, the factors influenced organic matter decomposition and SOC accumulation, thereby contributing to the heterogeneity observed in the spatial distribution of SOC content.

Given the necessity of conducting linear regression and ANOVA analyses between the dataset and the environmental variables, both of which relies on the fundamental assumption of the error term’s normality, necessitating a certain degree of data normality, the dataset underwent a square-root transformation, hereafter referred to as SQRT_SOC, measured in units of g/kg. The resulting open-square-root-transformed overall, modeling, and validation sets all showed a morphologically comparable, approximate normal distribution (Figure S2).

3.2. Environment Variables

3.2.1. Natural Factors Analyses

MAP, MAT, TT, ST, DEM, and Distowater have significant impacts on the spatial variability of SOC content (Figure 2 and Figure 3). MAT shows a negative relationship with SOC content at 0–10.65 °C. MAP was positively correlated with SOC when it was lower than 1269.45 mm. With increasing temperatures, soil microbial activity increases, and soil respiration accelerates, resulting in the fast decomposition of SOC. Conversely, when precipitation levels are lower, soil aeration improves, promoting the activity of aerobic microorganisms and thereby accelerating SOC decomposition [21]. Additionally, ST and TT are intimately tied to many soil physical and chemical properties, including soil texture, soil porosity, soil water content and pH, all of which directly impact OC content Chen et al. [21,35] further reveal that as DEM increases, SOC initially rises to a peak at 137 m, subsequently declining. In the rainy season, plains with inadequate drainage conditions are prone to flooding, resulting in waterlogging, which may be the reason why the positive correlation between distance from water sources and SOC content.

3.2.2. Cropping and Management Factors Analyses

Considering the incorporation of agricultural practices into research, it has been demonstrated that CT, TM, NDVI202202, and NDVImax all had notable impacts on the content of SOC (Figure 4 and Figure 5). The ranking of SOC content among various crops was as follows: soybean > other > rice > corn. The heightened SOC content in the soybean-cultivated soil can be attributed to the nitrogen fixation process, the well-developed root system, and the dense stems and leaves during the growth cycle, which collectively contributed organic matter to the soil and facilitate the accumulation of SOC.

The soil cultivated with maize has the lowest SOC content, potentially stemming from the incessant cropping of maize in the agricultural planting system, coupled with substantial water loss and erosion, both of which contributed to a decrease in SOC content. The implementation of crop rotation and fallow practices have important implications for SOC content. By judiciously selecting plant species and their order of cultivation, the decomposition of root residues can furnish a diverse array of organic matter sources, subsequently enhancing SOC content [36]. In addition, fallow tillage fosters the proliferation of soil microbial communities, bolstering soil biological activity and facilitating the gradual accumulation of organic matter, ultimately leading to an increase in SOC content.

NDVI serves as a valuable for monitoring the growth and health of vegetation. High NDVI indicates that robust vegetation growth, which, in turn, indirectly points to sufficient soil nutrients. Nevertheless, this study found that NDVI from December to February in the sampling year was observed to have a significant negatively correlation with SOC. This observation could potential attributed to a combination of factors, including the higher latitude, the lower winter temperature, and the higher snow cover, all of which result in lower NDVI.

3.2.3. Landscape Metrics Analyses

SHDI, PLAND4, PLAND5, PLAND6, and PLAND7 have a negative correlation with SOC content, while both of COHENSION and PLAND7 demonstrate positive correlation (Figure 6). The intricate landscape patterns hinder the migration of species across diverse patches, particularly those fragmented landscape by water bodies, roads and residential construction land, this fragmentation reduces the genetic diversity of soil organisms across regions and impedes the stabilization of SOC [30]. In addition, the complexity of the landscape is often due to the conversion of natural vegetation into farmland and the long-term intensive farming practices, which disrupt the structure of large soil aggregates and accelerates the loss of soil nutrients. Consequently, a highly fragmented landscape is detrimental to the accumulation of SOC.

In summary, the significant natural factors influencing the SOC contents in the studied area encompass soil type, soil texture, MAT, MAP, DEM, and distance to water. Among the farmland cultivation and management factors, Land use, crop type, tillage mode, NDVImax and NDVI202202 stand out. Additionally, the significant land-use landscape pattern indices include COHENSION, SHDI, PLAND4, PLAND5, PLAND6, and PLAND7. In conclusion, a total of fifteen significant influences on SOC within the studied area have been identified.

3.3. Spatial Distribution of Predicted SOC

3.3.1. Model Performance

RF regression model and Cubist model were constructed using factors that demonstrated a statistically significant correlation at the p < 0.05 level, as analyzed through Pearson’s correlation coefficient and ANOVA in the previous section of this study. Upon assessing the statistical metrics, we discerned that RF model, characterized by a higher mean value of R² and LCCC values alongside a lower RMSE, exhibited superior prowess in estimating observed SOC content compared to the Cubist model (Figure 7). Therefore, the RF model was deemed as the most suitable for the SOC content spatial distribution prediction.

The %InMSE represents the extent to which the model’s the MAE on the test set escalates upon the exclusion of each individual variable, thereby underlining the variable’s significance (Figure S3). The descending order of importance, the model’s sensitivity to various environmental variables influencing SOC was as follows: MAT > DEM > NDVI202202 > SHDI > COHENSION > MAP > PLAND6 > PLAND7 > ST > NDVImax > Distowater > CT > PLAND4 > TT > PLAND5 > TM.

The cubist model possesses the capability to delve into the intricate nonlinear inertia existing between independent and dependent variables. By strategically segmenting the dataset based on predefined rules, it ensures the key SOC regulating factors within distinct strata are duly considered. Specifically, the conditional influence of MAT stands at 70%, whereas DEM contributed 18%, indicating a disparity in the interplay between cultivated land SOC and environmental variables across varying MAT and DEM ranges. Among the four most influential factors contributing to the model’s accuracy, MAT, DEM, MAP, and PLAND7 emerge as prominent (Table S2), reinforcing the notion that the cubist models primary controlling factor remains rooted in natural element.

3.3.2. Spatial Distribution of SOC

According to the spatial distribution map of SOC content predicted by the RF model, the SOC content in cropland within the Northeast Black Soil Region ranged from 5.24 to 43.93 g/kg (Figure 8). The SOC content in cropland located in the Songnen and Sanjiang plains was generally higher, whereas it was relatively lower in the western section of the Liaohe Plain. This pattern aligns with the overall distribution observed at the sampling sites. The overarching distribution trend of SOC content depicted in the spatial distribution map mirrored the spatial pattern of mean annual air temperature and soil type, thereby reinforcing the notion that the annual air temperature serves as the primary determinant of SOC content in cropland within the Northeast Black Soil Region. The study exhibited a discernible gradient SOC content, which lower levels in the southwest and progressively higher levels towards the northeast. Consequently, the average SOC content across the region amounted to 17.24 g/kg.

3.4. Environmental Variable Analyses for SOC

3.4.1. Stepwise Regression Model Analyses

In order to further understanding of the underlying mechanism of each significant SOC influencing factor operates, as well as the collective interplay of various influencing factor on the SOC, we conducted separate stepwise regression analyses for natural elements, anthropogenic activity elements, and all the environmental variables, respectively (Table 3). The magnitude of the absolute value of the standardized coefficients in the stepwise regression model serves as an indicator of the degree to which a given environmental variable contributes to explaining the spatial variability of SOC (Figure S4).

In the SQRT_SOC stepwise regression model, using exclusively natural elements, we observed that 48.8% of the spatial variation in SOC is accounted (Adj. R² = 0.488). Among these factors, MAP emerges as the most dominant factor affecting the spatial variation in SOC, with a standardized coefficient of −0.587. Alternatively, approximately 43.0% of the spatial variation in SOC can be elucidated considering solely anthropogenic elements. In this regression framework, NDVI202202 stands out as the most important environmental variable exhibiting a standardized coefficient of −0.384. When all environmental variables are considered, the explanatory power increases to 56.6% of the spatial variance in SOC. MAT emerges as the strongest single variable in the explaining SOC variance. Therefore, in the Northeast Black Soil Region, natural factors slightly outweigh anthropogenic activities in elucidating spatial variation in SOC content of cropland. Nevertheless, both factors are equally crucial contribute to explaining the spatial variability of SOC content in the region’s cropland.

3.4.2. Structure Equation Modeling Analyses

In order to further dissect the collective influence of various environmental variables, both individually and collectively, on SOC, as well as the driving relationship among these, the stepwise regression model outcomes guided the selection of latent variables. These encompass, climate factors, geographic factors, soil attributes, planting and management factors, and landscape pattern indices. Specifically, MAP and MAT serve as the observational variables for climate factors, whereas DEM and the distance to the water bodies represent the geographic factors’ observable variables. Soil type and texture are utilized as observable indicators of soil attributes. The latent variables within planting and management factors encompass crop type, tillage mode, NDVI202202 and NDVImax. For landscape pattern indices, PLAND4, PLAND5, PLAND6, PLAND7, and LSI_C are the observed variables. Based on these variables, a SEM model of SOC’s primary governing factors was formulate, incorporating random effects (Figure 9).

The results demonstrated that the established SEM model aligns well with Fisher.C = 0.735, df = 2, and p = 0.692 (p > 0.05) (Figure 9). In certain scenarios, neglecting random effects can potentially skew the significance of a parameter or overlook critical variations within the data, ultimately leading to biased or distorted results. This consideration becomes paramount when dealing with experimental design incorporating repeated measures, multilevel data, or randomized groupings. The R²_Marginal metric indicates the proportion of total variance explained by fixed effects, with the he R²_Marginal of 0.52, while incorporating random effects elevates this explanatory power, as evidenced by an R²_Conditional of 0.94, significantly improving the model’s ability to account for the total variance. Additionally, aligning with the outcomes of the stepwise regression model, MAT was the primary driver influencing the individual variables under consideration.

Figure 10 illustrates the direct and indirect effects of the five latent variables within the SEM, where the climate factor t exerts the greatest direct effect, with a path coefficient of 0.546. The indirect effect of the climate factor is 0.051, which suggests that not only directly impacts the sequestration of SOC by soil respiration but also indirectly modulates the spatial variability of SOC by altering soil properties and influencing crop planting and management, which in turn affect the spatial variability of SOC. Despite the modest indirect effect of 0.051, the climate factor remains the dominant latent variable, contributing a total effect of 0.597, outranking the geographic factor with a total effect of 0.444. Notably, the geographic factor exerts is the largest indirect influence, affecting SOC indirectly through its effect on SOC. Beyond its direct impact, the geographic factor further modulates SOC by mediating influences of climate, soil properties, and planting and management factors.

3.5. Stock and Sequestration Potential of SOC

3.5.1. Spatial Distribution of Predicted SOCD

The SOCD of 0–15 cm of cropland in the Northeast black soil zone ranged from 0.51 to 9.11 kg/m² with an average value of 3.30 kg/m² (Figure 11). Accordingly, the current SOC storage capacity within this depth range in the region was estimated to be 1226.64 Tg. The SOCD was derived from SOC content according to the formula of Schwager and Mikhailova (2002), which comprehensively took into account the soil thickness, soil bulk weight and gravel content [34]. The overall spatial distribution maps of SOC content and SOCD exhibited remarkable similarity, indicating a negligible disparity in their overall patterns. This congruency extends to the observation that regions with abundant SOC content also display elevated SOCD values. Specifically, lower SOCD levels were recorded in the southern Songnen Plain and Liaohe Plain, while higher SOCD values were prevalent across Heilongjiang Province and the eastern segment of Jilin Province, mirroring the spatial distribution of SOC content.

3.5.2. Assessment of SOCsp of Cropland

In this study, we opted to employ the maximum value methodology for assessing SOCsp. SOCsp and SOCpc were calculated for each soil type according to Equations (6)–(8) (Table 4). The maximum values of SOCD in each soil type were, in order, Luvisols, Semi-hydromorphic soils, Hydromorphic soils, Semi-Luvisols, Caliche Soils, Skeletol primitive soils, Anthrosols and Saline soils. The SOCsp in the topsoil (0–15 cm) of the Northeast Black Soil Region, within agricultural lands, was projected to reach 3057.65 Tg, with a potential sequestration capacity of 1831.01 Tg determined by the maximum SOCD method for each soil type method. This SOCsp was predominantly concentrated in four major types: Hydromorphic soils, Luvisols, Caliche Soils, and Semi-Luvisols, which accounted for 32.22%, 25.28%, 15.71%, and 15.66% of the total sequestration potential, respectively.

4. Conclusions

The SOC content of the 555 sampling sites varied from 2.55 to 72.39 g/kg, with an average value of 19.20 g/kg. Our study identified a total of 16 significant influencing factors, including crucial natural factors (soil type, soil texture, MAT, MAP, DEM, and Distowater), significant farmland cultivation and management practices (land use, crop type, tillage mode, NDVI202202, and NDVImax), and notable land use landscape pattern indices (COHENSION, SHDI, PLAND4, PLAND5, PLAND6, and PLAND7). Among the three spatial distribution models constructed for SOC content, the RF model displayed as the most effective, exhibiting an R² of 0.70 and an RMSE of 0.62 g/kg. Consequently, it was chosen for predicting the spatial distribution of SOC content in this study. The predicted SOC content ranged from 5.24 to 43.93 g/kg, with a mean value of 17.24 g/kg, mirroring the overall distribution pattern observed at the sampling sites, where the overall trend was predominantly shaped by the climate and topographical factors. Additionally, the SOCD in the studied area ranged from 0.51 to 9.11 kg/m², and the average value was 3.30 kg/m², translating to a current SOC reserve of 1226.64 Tg and a SOCsp of 3057.65 Tg. This analysis underscores the carbon sequestration and sink enhancement potential of the northeastern black soil region, a vital terrestrial carbon reservoir aimed at achieving carbon neutrality. We hope to transform the theoretical value of SOC sequestration potential into actual SOC sequestration capacity by promoting sustainable agriculture and additional strategies.

Author Contributions

W.L. contributed to the conceptualization, data curation, methodology, validation, visualization, writing—original draft, and writing—review and editing. Z.Y. contributed to the writing—review and editing. J.J. and G.S. contributed to the conceptualization, funding acquisition, supervision, and writing—reviewing and editing. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are available on request from the authors.

Conflicts of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Locations of 555 soil sample points used for model development in the black soil region in Northeast China.

Figure 2. Correlation matrix of Pearson correlation coefficient (upper triangle) and scatter plot matrix among SOC_SQRT (N = 555) and nature factors at 300 m scale (lower triangle) in the study area. * indicates p [less than] 0.05, ** p [less than] 0.01, and *** p [less than] 0.001.

Figure 3. Box plot and one-way ANOVA of different nature factor with SQRT_SOC in study area (Soil type including: 1. Luvisols; 2. Semi-Luvisols; 3. Caliche Soils; 4. Skeletol primitive soils; 5. Dark Semi-hydromorphic soils 6. Hydromorphic soils; 7. Saline soils; 8. Anthrosols). If there is no same letter in the lower-case letters indicated for both factors, it means that they are significantly different from each other at the 0.05 level.

Figure 4. Correlation matrix of Pearson correlation coefficient (upper triangle) and scatter plot matrix among SOC_SQRT (N = 555) and the different cropping and management factors at 300 m scale (lower triangle) in the study area. * indicates p [less than] 0.05, ** p [less than] 0.01, and *** p [less than] 0.001.

Figure 5. Box plot and one-way ANOVA of different cropping and management factors with SQRT_SOC in study area (tillage mode including: 1. Maize; 2. Rice; 3. Soybean; 4. Maize-Soybean; 5. Maize-Rice; 6. Maize-Rice-Soybean; 7. Rice-Soybean; 8. Fallow; 9. Other crop). If there is no same letter in the lower-case letters indicated for both factors, it means that they are significantly different from each other at the 0.05 level.

Figure 6. Correlation matrix of Pearson correlation coefficient (upper triangle) and scatter plot matrix among SOC_SQRT (N = 555) and the different landscape metric at 300 m scale (lower triangle) in the study area (PLAND1 to PLAND7, respectively, represent the proportional areas of paddy fields, dryland, forests, grasslands, water, urban and rural land, and unused land). * indicates p [less than] 0.05, ** p [less than] 0.01, and *** p [less than] 0.001.

Figure 7. RF (a) and Cubist (b) model fitting effect and prediction effect evaluation.

Figure 8. Spatial distribution of SOC content in cropland in the study area.

Figure 9. The responses of SOC stock to environmental variables as derived from the SEM. Note: The thickness of the arrows indicates the magnitude of the effect, with red representing a positive effect, blue a negative effect, and dashed lines indicating not significant; * indicates p [less than] 0.05, ** p [less than] 0.01, and *** p [less than] 0.001.

Figure 10. Standardized effect of each latent variable in SEM on SOC content.

Figure 11. Spatial distribution of cropland SOCD in the northeast black soil region.

Table 1

The environmental variables of SOC.

Drivers	Environmental Variables	Resolution	Abbreviation	Unit	Data Origin and Reference
Natural factors	Soil type	1 km	ST	-	https://www.resdc.cn/ (accessed on 20 March 2023)
	Soil texture	1 km	TT	%
	Mean annual temperature	1 km	MAT	°C
	Mean annual precipitation	1 km	MAP	mm
	Elevation	30 m	DEM	m	https://www.gscloud.cn/ (accessed on 17 September 2023)
	Slope	30 m	Slope	°
	Distance to water	300 m	Distowater	m
	Distance to green space	300 m	Distogreen	m
Cropping and management factors	Crop type	10 m	CT	-	[25]
	Cropland land-use	30 m	CL	-	https://www.resdc.cn/ (accessed on 20 March 2023)
	Normalized difference vegetation index	250 m	NDVI	-
	Irrigation	30 m	IR	-	[26]
	Tillage modes	10 m	TM	-
Landscape metrics	Percentage of landscape	300 m	PLAND	%
	Patch cohesion index	300 m	COHESION	-
	Shannon’s diversity index	300 m	SHDI	-

Table 2

Descriptive statistical results of soil organic carbon content (g/kg).

Data	Number	Mean	SD	Skew	CV (%)	Maximum	Minimum	Median
Dataset	555	19.20	10.67	1.27	55.58	72.39	2.55	17.34
Trainingset	444	19.13	10.71	1.01	56.01	72.39	2.54	17.35
Validation set	111	19.48	10.54	1.07	54.11	52.57	5.42	17.05

Table 3

Stepwise regressions of different types of factors with SQRT_SOC.

Variable Types	Variable	Regression			Parameters
		F	Sig.	Adj. R2	Ratio	Error	t	Sig.
Naturefactors	Constant	122.828	<0.001	0.488	4.819	0.176	27.434	<0.001
	MAT				−0.251	−0.263	−21.074	<0.001
	MAP				0.01	0.000	5.147	<0.001
	Distowater				0.00002	0.000	5.356	<0.001
	Soiltype4				0.743	0.225	−3.034	<0.001
	Soiltype3				−0.445	0.085	−5.242	<0.001
	Soiltype7				−1.342	0.591	−2.271	<0.05
	Soiltexture1				−0.429	0.193	−2.227	<0.05
Humanactivitiesfactors	Constant	69.573	<0.001	0.430	1.814	0.456	3.982	<0.001
	NDVI202202				−5.046	0.486	−10.376	<0.001
	Croptype3				0.787	0.114	6.879	<0.001
	NDVImax				3.293	0.523	6.300	<0.001
	Tillagemode8				0.457	00.207	2.205	<0.05
	Tillagemode9				−0.913	0.345	−2.650	<0.05
	Tillagemode3				−0.231	0.099	−2.346	<0.05
Allfactors	Constant	80.232	<0.001	0.566	3.891	0.446	8.728	<0.001
	MAT				−0.255	0.017	−15.233	<0.001
	Soiltype3				−0.450	0.101	−4.468	<0.001
	Soiltype4				−0.888	0.203	−4.370	<0.001
	Distowater				0.00003	0.000	4.791	<0.001
	NDVImax				1.828	0.106	3.807	<0.001
	Soiltype1				0.285	0.091	3.133	<0.05
	Croptype1				−0.247	0.075	−3.283	<0.05
	Croptype4				−0.450	0.191	−2.352	<0.05
	PLAND6				0.013	0.006	2.161	<0.05

Note: Dummy variables were set for the category variables, where Soiltype1–7 indicated Luvisols, Semi-Luvisols, Caliche Soils, Skeletol primitive soils, Dark Semi-hydromorphic soils, Hydromorphic soils, Saline soils, Anthrosols, respectively; Soiltexture1–3 indicated Sandy, Clay Loam, and Clay; Croptype1–4 indicated Maize, Rice, Soybean, and Others, respectively; Tillagemode1–9 indicated Maize, Rice, Soybean, Maize-Soybean, Maize-Rice, Maize-Rice-Soybean, Rice-Soybean, Fallow and Other crop, respectively; and Distowater refers to the distance to the water body.

Table 4

SOC sequestration potential of 0–15 cm soil of different soil types.

Soil Type	SOCDmean(kg/m2)	SOCDmax (kg/m2)	Area(109 m2)	SOCs(Tg)	SOCsp(Tg)	SOCpc(Tg)
Luvisols	3.88	9.11	84.86	329.26	773.07	443.81
Semi-Luvisols	3.72	7.90	60.61	225.47	478.82	253.35
Caliche Soils	2.65	7.81	61.52	163.03	480.47	317.44
Skeletol primitive soils	2.31	7.44	16.13	37.26	120.01	82.75
Semi-hydromorphic soils	3.14	8.30	118.71	372.75	985.29	612.54
Hydromorphic soils	3.91	8.29	16.25	63.54	134.71	71.17
Saline soils	1.97	5.30	1.68	3.31	8.90	5.59
Anthrosols	2.78	6.63	11.52	32.03	76.38	44.35
Sum				1226.64	3057.65	1831.01

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/agronomy14112744/s1, Figure S1: Spatial distribution trends in SOC content (a) Latitude-longitude-SOC three-dimensional map; (b) Latitude-SOC scatter; (c) SOC contour map; (d) Longitude-SOC scatter); Figure S2: Dataset, training set and test set SQRT_SOC normal distribution plots; Figure S3: The contribution rate of each variable of the random forest model; Figure S4: Ranking of standardized coefficients for each environmental variable in stepwise regression models for each type and all environments; Table S1: The landscape metrics and meaning description used in this study; Table S2: Conditional contribution rate and model contribution rate of each factor in the Cubist.

References

1. Schlesinger, W.H. On fertilizer-induced soil carbon sequestration in China’s croplands. Glob. Change Biol.; 2010; 16, pp. 849-850. [DOI: https://dx.doi.org/10.1111/j.1365-2486.2009.01958.x]

2. Qin, Z.C.; Huang, Y.; Zhuang, Q.L. Soil organic carbon sequestration potential of cropland in China. Glob. Biogeochem. Cycles; 2013; 27, pp. 711-722. [DOI: https://dx.doi.org/10.1002/gbc.20068]

3. Lal, R. Soil carbon sequestration impacts on global climate change and food security. Science; 2004; 304, pp. 1623-1627. [DOI: https://dx.doi.org/10.1126/science.1097396]

4. Bangroo, S.A.; Najar, G.R.; Rasool, A. Effect of altitude and aspect on soil organic carbon and nitrogen stocks in the Himalayan Mawer Forest Range. Catena; 2017; 158, pp. 63-68. [DOI: https://dx.doi.org/10.1016/j.catena.2017.06.017]

5. Teixeira, R.F.M.; Morais, T.G.; Domingos, T. Global process-based characterization factors of soil carbon depletion for life cycle impact assessment. Sci. Data; 2021; 8, 237. [DOI: https://dx.doi.org/10.1038/s41597-021-01018-2]

6. Zuo, W.A.; Gu, B.X.; Zou, X.W.; Peng, K.; Shan, Y.L.; Yi, S.Q.; Shan, Y.H.; Gu, C.H.; Bai, Y.C. Soil organic carbon sequestration in croplands can make remarkable contributions to China’s carbon neutrality. J. Clean. Prod.; 2023; 382, 135268. [DOI: https://dx.doi.org/10.1016/j.jclepro.2022.135268]

7. IEA. CO₂ Emissions in 2023: A New Record High, But Is There Light at the End of the Tunnel?; International Energy Agency: Paris, French, 2024; Available online: https://www.iea.org/reports/co2-emissions-in-2023 (accessed on 30 March 2024).

8. Qiao, L.; Wang, X.H.; Smith, P.; Fan, J.L.; Lu, Y.L.; Emmett, B.; Li, R.; Dorling, S.; Chen, H.Q.; Liu, S.G. et al. Soil quality both increases crop production and improves resilience to climate change. Nat. Clim. Chang.; 2022; 12, pp. 574-580. [DOI: https://dx.doi.org/10.1038/s41558-022-01376-8]

9. Ma, Y.Q.; Woolf, D.; Fan, M.S.; Qiao, L.; Li, R.; Lehmann, J. Global crop production increase by soil organic carbon. Nat. Geosci.; 2023; 16, pp. 1159-1165. [DOI: https://dx.doi.org/10.1038/s41561-023-01302-3]

10. Odebiri, O.; Mutanga, O.; Odindi, J.; Naicker, R.; Masemola, C.; Sibanda, M. Deep learning approaches in remote sensing of soil organic carbon: A review of utility, challenges, and prospects. Environ. Monit. Assess.; 2021; 193, 802. [DOI: https://dx.doi.org/10.1007/s10661-021-09561-6]

11. Wang, M.M.; Guo, X.W.; Zhang, S.; Xiao, L.J.; Mishra, U.; Yang, Y.H.; Zhu, B.A.; Wang, G.C.; Mao, X.L.; Qian, T. et al. Global soil profiles indicate depth-dependent soil carbon losses under a warmer climate. Nat. Commun.; 2022; 13, 5514. [DOI: https://dx.doi.org/10.1038/s41467-022-33278-w]

12. Luo, C.; Zhang, W.; Zhang, X.; Liu, H. Mapping of soil organic matter in a typical black soil area using Landsat-8 synthetic images at different time periods. Catena; 2023; 231, 107336. [DOI: https://dx.doi.org/10.1016/j.catena.2023.107336]

13. Guo, L.; Fu, G.; Shi, T.; Chen, Y.; Zeng, C.; Zhang, H.; Wang, S. Exploring influence factors in mapping soil organic carbon on low-relief agricultural lands using time series of remote sensing data. Soil Tillage Res.; 2021; 210, 104982. [DOI: https://dx.doi.org/10.1016/j.still.2021.104982]

14. Lal, R. World cropland soils as a source or sink for atmospheric carbon. Adv. Agron.; 2001; 71, pp. 145-191. [DOI: https://dx.doi.org/10.1016/S0065-2113(01)71014-0]

15. Bowie, M.H.; Black, L.; Boyer, S.; Dickinson, N.M.; Hodge, S. Persistence of biodiversity in a dryland remnant within an intensified dairy farm landscape. N. Z. J. Ecol.; 2016; 40, pp. 121-130. [DOI: https://dx.doi.org/10.20417/nzjecol.40.14]

16. Didham, R.K.; Barker, G.M.; Bartlam, S.; Deakin, E.L.; Denmead, L.H.; Fisk, L.M.; Peters, J.M.; Tylianakis, J.M.; Wright, H.R.; Schipper, L.A. Agricultural intensification exacerbates spillover effects on soil biogeochemistry in adjacent forest remnants. PLoS ONE; 2015; 10, e0116474. [DOI: https://dx.doi.org/10.1371/journal.pone.0116474]

17. Amelung, W.; Bossio, D.; de Vries, W.; Kogel-Knabner, I.; Lehmann, J.; Amundson, R.; Bol, R.; Collins, C.; Lal, R.; Leifeld, J. et al. Towards a global-scale soil climate mitigation strategy. Nat. Commun.; 2020; 11, 5427. [DOI: https://dx.doi.org/10.1038/s41467-020-18887-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33110065]

18. Zhu, K.; Ran, H.Y.; Wang, F.F.; Ye, X.; Niu, L.A.; Schulin, R.; Wang, G. Conservation tillage facilitated soil carbon sequestration through diversified carbon conversions. Agric. Ecosyst. Environ.; 2022; 337, 108080. [DOI: https://dx.doi.org/10.1016/j.agee.2022.108080]

19. Yang, X.L.; Xiong, J.R.; Du, T.S.; Ju, X.T.; Gan, Y.T.; Li, S.; Xia, L.L.; Shen, Y.J.; Pacenka, S.; Steenhuis, T.S. et al. Diversifying crop rotation increases food production, reduces net greenhouse gas emissions and improves soil health. Nat. Commun.; 2024; 15, 198. [DOI: https://dx.doi.org/10.1038/s41467-023-44464-9]

20. Gorooei, A.; Aynehband, A.; Rahnama, A.; Gaiser, T.; Kamali, B. Cropping systems and agricultural management strategies affect soil organic carbon dynamics in semi-arid regions. Front. Sustain. Food Syst.; 2023; 6, 1016000. [DOI: https://dx.doi.org/10.3389/fsufs.2022.1016000]

21. Chen, F.Z.; Feng, P.Y.; Harrison, M.T.; Wang, B.; Liu, K.; Zhang, C.X.; Hu, K.L. Cropland carbon stocks driven by soil characteristics, rainfall and elevation. Sci. Total Environ.; 2023; 862, 160602. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2022.160602]

22. CAS. Northeast Black Land Protection and Utilization Report (2021); China academic of science: Beijing, China, 2022; Available online: https://www.cas.cn/sygz/202210/t20221002_4849787.shtml (accessed on 30 March 2023).

23. Shi, J.P.; Song, G. Native Database in China--Based on the Second National Soil Soil Survey Data Sets [DS/OL]. V1. Science Data Bank, 2016. Available online: https://www.scidb.cn/en/detail?dataSetId=633694460877733888 (accessed on 4 April 2024).

24. Gong, Z. Chinese Soil Taxonomy; China Science Press: Beijing, China, 2001.

25. You, N.; Dong, J.; Huang, J.; Du, G.; Zhang, G.; He, Y.; Yang, T.; Di, Y.; Xiao, X. The 10-m crop type maps in Northeast China during 2017–2019. Sci. Data; 2021; 8, 41. [DOI: https://dx.doi.org/10.1038/s41597-021-00827-9] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33531510]

26. Zhang, C.; Dong, J.; Ge, Q. Mapping 20 years of irrigated croplands in China using MODIS and statistics and existing irrigation products. Sci. Data; 2022; 9, 407. [DOI: https://dx.doi.org/10.1038/s41597-022-01522-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35840621]

27. Xu, X.L. Annual spatially interpolated dataset of meteorological elements in China. Resour. Environ. Sci. Data Regist. Publ. Syst.; 2022; [DOI: https://dx.doi.org/10.12078/2022082501]

28. Gscloud. ASTER GDEM 30M Resolution Digital Elevation Data. Geospatial Data Cloud. 2023; Available online: https://www.gscloud.cn/#page1/2 (accessed on 30 October 2022).

29. LAADS&DAAC. MODIS—Level-0 and Level-1. The Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC). 2024; Available online: https://ladsweb.modaps.eosdis.nasa.gov/ (accessed on 1 March 2023).

30. Liu, X.C.; Bian, Z.X.; Sun, Z.T.; Wang, C.Q.; Sun, Z.Q.; Wang, S.; Wang, G.L. Integrating Landscape Pattern Metrics to Map Spatial Distribution of Farmland Soil Organic Carbon on Lower Liaohe Plain of Northeast China. Land; 2023; 12, 1344. [DOI: https://dx.doi.org/10.3390/land12071344]

31. Saidi, S.; Ayoubi, S.; Shirvani, M.; Azizi, K.; Zhao, S. Digital mapping of soil phosphorous sorption parameters (PSPs) using environmental variables and machine learning algorithms. Int. J. Digit. Earth; 2023; 16, pp. 1752-1769. [DOI: https://dx.doi.org/10.1080/17538947.2023.2210314]

32. Lefcheck, J.S. PIECEWISESEM: Piecewise structural equation modelling in R for ecology, evolution, and systematics. Methods Ecol. Evol.; 2016; 7, pp. 573-579. [DOI: https://dx.doi.org/10.1111/2041-210X.12512]

33. Tarka, P. An overview of structural equation modeling: Its beginnings, historical development, usefulness and controversies in the social sciences. Qual. Quant.; 2018; 52, pp. 313-354. [DOI: https://dx.doi.org/10.1007/s11135-017-0469-8]

34. Schwager, S.J.; Mikhailova, E.A.J.S.S. Estimating Variability in Soil Organic Carbon Storage Using the Method Of Statistical Differentials. Soil Sci.; 2002; 167, pp. 194-200. [DOI: https://dx.doi.org/10.1097/00010694-200203000-00004]

35. Wu, Z. Research on Spatial Variation and Mechanism of Farmland Soil Organic Carbon in Plains, Wuhan, China. Ph.D. Thesis; Wuhan University: Wuhan, China, 2021.

36. Hu, N.J.; Shi, H.; Wang, B.J.; Gu, Z.H.; Zhu, L.Q. Effects of different wheat straw returning modes on soil organic carbon sequestration in a rice-wheat rotation. Can. J. Soil Sci.; 2019; 99, pp. 25-35. [DOI: https://dx.doi.org/10.1139/cjss-2018-0094]