The study was conducted in Longxi County (34°50′-35°23′N, 104°18′-104°54′E), which is in the middle of the Loess Plateau, in northwest China ( Figure 1 ). There is a typical continental climate, with four distinct seasons, abundant sunshine and a mild climate. The mean annual precipitation is 445.8 mm, while the mean annual average temperature is 6.54^°C. Influenced by the general circulation of East Asia and the special topography around the Qinghai Tibet Plateau. The area is located between the edge of loess plateau in northwest China and the foothills of the Qinling mountains. The terrain is high in the northwest and low in the southeast and the altitude ranges from 1612 to 2762m. Longxi County can be divided into three main climatic regions, namely, the temperate-warm semi-arid region which has a relatively low altitude, mild climate, more accumulated temperature, longer frost-free period and rich light resources, the temperate-cool semi-arid region where the climate is warm and cool with moderate heat but is colder in winter, and the temperate-cold semi-humid area, compared with the other areas, in which precipitation is abundant, but heat is insufficient and accumulated temperature is the least. The loess hilly, mountain and river plain are the most typical landforms. The agrotype in Longxi Country can be divided into five soil types (1): Ustochrept, formed by direct cultivation and ripening on the parent material of loess, without obvious profile development level, having good water storage, but poor soil nutrients and anti-scour ability (2); Chernozem, as a special soil type developed in the semi-arid climate of warm temperate zone with aeolian loess as parent material and grassland xerophytic vegetation as main vegetation forms, it has a uniform texture but poor permeability (3); Krasnozem, the parent material is tertiary red layer, being coarse in texture and poor in water and fertility preservation (4); Aquoll, as a kind of micro-regional soil formed by river sediments, having poor cultivability and low fertility (5); Alfisols, as a kind of forest soil in the vertical zone of mountain, with complex parent materials including Loess, Tertiary sediments, Cretaceous sediments and residual slope deposits of Paleozoic rocks, having high humus content and sticky texture. As the major producing areas of Chinese medicinal materials in Gansu Province, and even in China, more than 70 Chinese medicinal herbs are cultivated there, including 50 wild species and 20 cultivated species, such as Codonopsis pilosula (Franch.) Nannf, Astragalus propinquus Schischkin, Scutellaria baicalensis Georgi, Glycyrrhiza uralensis Fisch, and Bupleurum spp., among others (39).

Longxi County mainly contains planted traditional Chinese medicinal herbs, and the planting area accounts for about 70% of the cultivated land. In 2018, soil sampling was conducted at 200 locations in 134 villages of 17 towns in which traditional Chinese medicinal herbs were cultivated ( Figure 1 ). The geographic co-ordinates of the sampling locations were recorded using a hand-held GPS set. Soil sampling points were randomly set in each location, and the distance among pairwise points was more than 2m. Three mixed soil samples at 0–10 cm, 10–30 cm and 30–50 cm soil depths were collected in each sampling location. Soil bulk densities were determined by cutting ring method in site. In total, 600 soil samples were collected. The weeds, roots and gravel in soil were eliminated from the air-dried soil samples, and then sieved using a 2 mm sieve. The SOC was determined via potassium dichromate titration and calculated using the following formula (40):

Where SOC is the soil organic carbon content, g kg^-1; c is the standard solution concentration of 0.8000 mol L^-1 (1/6 K₂Cr₂O₇); 5 is the volume of the potassium dichromate standard solution added, mL; V₀ is the blank titration volume of FeSO₄, mL; 3.0 is the molar mass of 1/4 carbon atom, g mol^-1; 1.1 is the oxidation correction factor; m represents the quantity of the air-dried soil samples, g; and k is the conversion coefficient of the dried soil.

Where SOCD_i is the soil organic carbon density, 10³ kg ha^-1; SOC_i is the soil organic carbon content, g kg^-1; D_i is the soil bulk density; and E_i is the thickness of the soil layer, cm.

where SOCR_i is the organic carbon storage on each soil layer, kg; SOCD_i represents the organic carbon density of each soil layer, kg m^-2; and S_i is the pixel area, m².

The meteorological data (1958–2018) was provided by China Meteorological Science data sharing Service Network (http://cdc.nmic.cn) and the local meteorological bureau. In order to fully exhibit the spatial distribution of temperature (AT) and precipitation (PR), the multiple regression equations (formula 4 and 5) among the meteorological longitude (LON), latitude (LAT), altitude (ALT), AT and PR were established based on the data of 14 meteorological stations, which assisted us in deducing the meteorological data of each sampling plot.

The digital elevation model (DEM) data were derived from the geospatial cloud platform (www.gscloud.cn/). ArcGIS 10.6 software was used to obtain the basic topographic data of Longxi County, including the elevation, slope, aspect, surface roughness, relief amplitude, slope of aspect (SOA) and composite terrain variable parameters, which were represented by the stream power index (SPI), sediment transport index (STI) and topographic wetness index (TWI). The SPI denoted the moving capacity of the runoff to surface materials of the slope; STI mirrored the transport status of the surface sediments; TWI is an index that describes the soil moisture distribution ^[7]. The formula used was as follows:

A_S represents the specific catchment area of any point on the surface; β represents the slope of the point.

Data analyses were carried out using R 4.0.4 with its add-in packages: “car”, “Box-Cox”, “MASS” and “multcomp”. One-way ANOVA was adopted to determine the differences in SOCD in the different townships and soil layers. The normal distribution of SOCD was analyzed using the Kormolov-Sminov test. Box-Cox function package was applied to calculate the normal conversion coefficient, and the bcPowerBox-Cox function package was adopted to carry out the conversion. The conversion formula used was as follows:

where NSOCD_i represents the soil organic carbon density of each soil layer; SOCD_i is the density of soil organic carbon in each soil layer; and λ_i stands for the conversion coefficient of each soil layer.

To eliminate the influence of dimensionality, we standardized the environmental variables, and then analyzed the correlation between the SOCD and climate and terrain factors using the Pearson correlation coefficient. The regression model between SOCD and the environmental factors was built using a multiple linear regression. Principal component analysis was used to transform the data. The principal components were extracted using the Kaiser-Harris criterion, and the eigenvalues greater than 1 were retained. According to the score coefficient of the environmental variables on each principal component, both the regression expression of each principal component variable and environmental variable could be obtained. The calculated principal component variables of F1 _i and F2 _i were taken as new independent variables, which were then incorporated into the regression modelling process of principal components and NSOCD; thus, the principal component regression fitting models were established through synthesizing the information of the main environmental variables. The residual values of principal component regression analysis were interpolated using ordinary Kriging and Arcgis 10.6. The formula used is the following:

where λ represents the normal conversion coefficient; principal component regression Kriging (PCRK) represents the SOCD based on the principal component regression Kriging model prediction; PCA stands for the SOCD values predicted using principal component regression models; REKR represents the SOCD values predicted via Kriging interpolation for principal component regression residuals; a represents the regression coefficient; ∑ 1 n b i × F i represents the weighted sum of n principal components; F_i is the I principal component, b_i represents the principal component regression coefficient; and n is the principal component score.

Random sample data collected from 30 sampling points were selected as validation datasets and the remaining 170 samples served as training datasets. The mean error (ME), mean absolute error (MAE), root mean square error (RMSE) and relative accuracy (RA) were selected to estimate the prediction accuracy of the model. When the ME value is greater than 0, the prediction value is lower than the actual measured value; when the ME value is lower than 0, the prediction value is higher than the measured value. The values of MAE, RMSE, or RA were used to evaluate the prediction accuracy of the model; the smaller the value, the higher the accuracy.

The SOCD in the 0–50 cm soil layer in Longxi County ranged from 15.68 Mg ha^-1 to 38.20 Mg ha^-1 ( Table 1 ). The average SOCD was 6.60, 11.52 and 10.18 Mg ha^-1 in the 0–10 cm 10–30 cm and 30–50 cm soil layer, respectively. Except for Fuxing Town, Shuangquan Town and Tong’anyi Town, the SOCD in the 10–30 cm soil layer was the highest, which was 1.07 to 1.44 times that of the 30–50 cm soil layer, and 1.28 to 2.03 times that of the 0–10 cm soil layer. The SOCD in the 30–50 cm soil layer took the second place and was 1.15–1.58 times that of the 0–10 cm soil layer, and the lowest level of SOCD was observed in the 0–10 cm soil layer. The mean standard deviation (σ) increased gradually with the soil depth, indicating that the variations in SOCD became more and more significant with the increase in soil depth.

The SOCD of each soil layer showed non-normality, and the normality conversion coefficient λ of the 0–10, 10–30 and 30–50 soil layers were 0.7, 0.25 and 0.3, respectively. The results of correlation analysis showed that the temperature, precipitation, elevation and TWI were significantly correlated with the NSOCD (P< 0.05), while the NSOCD_10cm was significantly (P< 0.05) correlated with the slope, roughness and relief amplitude. The NSOCD_30cm was also significantly correlated with the slope, and the NSOCD_50cm was significantly correlated with the TWI (P< 0.05, Table 2 ).

The correlation analysis results showed that a total of eight environmental variables were related to the NSOCD, of which the NSOCD_10cm, NSOCD_30cm and NSOCD_50cm was related to seven, five, and five environmental variables, respectively ( Table 3 ). The results of principal component analysis showed that the first two principal components related to NSOCD_10cm, NSOCD_30cm, and NSOCD_50cm contained 89.41%, 91.38% and 89.41% of the information. F1 _i and F2 _i were used to represent the two principal components of each soil layer. (Equations 12, 13, and 14).

The dimensionality reduction of the environmental variables was performed through correlation analysis and principal component analysis, and the relationship between the environmental variables and the SOCD was fully excavated by using a principal component regression model. However, the regression model only considers the structure of the data, which has a low degree of explanation for the SOCD, and the residual part cannot be fully explained. Therefore, a Kriging interpolation was carried out for the residual regression model, and then the residual results and the fitting results of the regression model were added together to obtain the PCRK to predict the spatial distribution characteristics of the SOCD ( Table 4 ).

The results showed that the total SOCD (0-50cm) was about 26.7 ± 10.2 Mg ha^-1 and the estimated stock of soil organic carbon reached 6.29 ± 2.41 Tg. The SOCD and stocks were different at different depths; in the first layer (0–10 cm), second layer (10–30 cm), and in the last layer (30–50 cm), the average SOCD reached 6.3 ± 1.7, 11.0 ± 3.8 and 9.3 ± 4.8 Mg ha^-1 with a stock of 1.49 ± 0.41, 2.58 ± 0.90 and 2.21 ± 1.11 Tg, respectively ( Figure 2 ).

Longxi County, which is located in the Loess Plateau, is the most important agricultural and traditional Chinese medicine planting area in Gansu Province. Its SOCD (0-20cm) level was 10.9 Mg ha^-1, which was lower than the cropland of the Loess Plateau and the average levels of the Chinese croplands (20.5 to 21.7 Mg ha^-1). The main soil types of Longxi County include Ustochrept, Chernozem, Krasnozem, Aquoll, and Alfisols, accounting for 82.99% of the county area (42, Figure 1 ). Among them, Ustochrept, chernozem and krasnozem show a coarse texture and serious leakage of water and fertilizer, which is not conducive to root development, resulting in a low content of soil organic matter (41). The SOCD of above-mentioned three soil types were 11.4, 18.6 and 16.1Mg ha^-1 (42), respectively. The area of these three soil types in Longxi County accounts for 75.16% of the total area, which is an important cause of the low SOCD levels in the tillage layer of croplands.

Climatic factors (temperature, precipitation) and topographic factors (elevation, slope, surface roughness, among others) are important factors affecting the SOC content. These factors have direct/indirect effects on the mineralization rate of SOC, microbial activity, aboveground biomass, soil erosion and human activities, resulting in spatial heterogeneity of organic carbon (5, 20, 43). Results showed that there was a significant negative correlation between the SOCD and temperature. With the increase in temperature, both heterotrophic respiration and CO₂ emissions increased, which led to a decrease in the SOCD levels (44–46). The rainfall is one of the main sources of soil water, affecting the soil moisture, permeability and the soil redox reactions, including mineralization, synthesis and decomposition of soil organic matter (47, 48). Studies had shown that the SOC content and precipitation were related, but the correlations varied from region to region (49). In this study, the SOCD in each soil layer was positively correlated with the annual rainfall. The reason is that the research area belongs to a semi-arid climate, and the increase in rainfall promotes the growth of vegetation and the decomposition of litter, which is conducive to the input and accumulation of organic carbon (50).

Topographic factors affected the redistribution of the soil hydrothermal resources and the process and intensity of material circulation in the soil ecosystem, thus affecting the organic carbon contents in soil (51, 52). This study showed that elevation and the SPI had a significant correlation with the SOCD in each soil layer, with the former being a significant positive correlation. With an increase in soil depth, the correlation between elevation and the SOCD was enhanced, and the correlation with the SPI was significantly negative. Climate change induced by altitude is a major determining factor; generally, a higher altitude leads to a lower temperature, less activity of the soil microbes, a lower decomposition rate of carbon, but a higher content of soil carbon (53, 54). The SPI is used to describe the erosivity of the surface flow. Higher levels of SPI indicate a greater extent of soil erosion, which ultimately leads to lower SOC levels (55). Additionally, the slope, relief amplitude and surface roughness can characterize the erosion of surface soil, the disturbance of the soil surface and the irregular changes in the micro-geomorphology (56). The surface soil was greatly affected by human activities, litters, crop residues, regular application of organic fertilizer and regular tillage (57). The behaviors mentioned above mostly affected the erosion status and micro topography of the soil surface, which resulted in a correlation between the slope, topographic fluctuation, surface roughness and the SOCD, which gradually decreased with an increase in the soil depth (8). Therefore, no obvious correlation was found in the deeper soil layer. In general, the TWI characterizes the dry and wet state of the soil and measures the water content of the soil. The results showed a negative correlation between the TWI and the SOCD with the increase in soil depth. Deep soil, which has a lower oxygen content and worse permeability, generally shows lower levels of SOCD and lower decomposition rates. For the agricultural land, human production activities have a great impact on soil, especially in traditional Chinese medicine planting areas; for example, different planting densities, technologies, fertilizer use, and continuous cropping or rotation tillage methods can affect the soil nutrients (6).