Identification and Evolution of Soil Organic Carbon Density Caused by Coastal Rapid Siltation Based on Imaging Spectroscopy

This article aimed to investigate the feasibility of using imaging spectroscopy (IS) to predict soil organic carbon density (SOCD) either directly or indirectly through soil organic carbon (SOC). Three methods, partial least square regression, support vector machine (SVM), and random forest, were utilized to calibrate the models and map the SOCD. The results showed that direct prediction was better than indirect prediction and the best model SVM had high R2 of 0.94 and 0.93 for calibration and validation, respectively. The measured SOCD was mainly concentrated in the surface soil layer from 0 to 40 cm, and the deeper layer tended to be gentle. The continuous depth variation trend of SOCD in the topsoil (0–40 cm) was relatively close to the measured values, while in the deeper layers (below 40 cm), it was much higher than the measured values. This article also found that the best fitting function of measured SOC stocks over time varied from linear to power and then to logarithmic with increasing depth, indicating less efficient accumulation of SOC in deep soil compared to topsoil, resulting in an overall decrease in SOC accumulation rate across soil depth. The best temporal functions for the predicted values differed from the measured values at each depth, but the changing trends of the three functions were basically consistent. It suggests that IS technology has the potential to quantitatively reveal the process of coastal soil evolution and offers a new insight for the rapid monitoring of soil property changes.


I. INTRODUCTION
S OIL organic carbon (SOC) is a crucial indicator of soil fertility and has a direct impact on crop and vegetation yields in terrestrial landscapes worldwide [1]. Increasing SOC stocks can offset some of the effects of fossil fuel carbon emissions and contribute to mitigating the greenhouse gas effect [2]. Therefore, SOC is often regarded as one of the most essential components to address the issues related to world food security and global climate change [3], [4], [5]. Understanding the dynamic changes in the regional SOC is conducive to regulating farming practices and land management strategies in many countries including China [6].
Until recently, studies on SOC have mainly focused on the surface layer (0-30 cm), with limited understanding of the deep soil profiles [7]. However, more than 50% of the organic carbon is stored at the depth between 0.3 m and 1 m [5], hence, soil sampling from an appropriate depth is essential for accurately calculating the rate, and potential of the soil profile in SOC sequestration [8], [9]. In addition, despite their high accuracy, the determination of SOC density (SOCD) relies mainly on conventional methods, which require considerable manpower and resources, and have the disadvantages such as longer time sampling with, higher sampling cost, as well as pollution [10]. The development of remote sensing technology has addressed these limitations, enabling real-time, rapid quantitative monitoring of soil properties. This capability holds greater importance for quickly obtaining soil information under the ongoing global climate change scenario [11], [12], [13].
Reflectance spectroscopy has been widely used to estimate SOC and other properties worldwide [14], [15], [16], [17], [18], [19], [20]. However, interestingly the reflectance spectroscopy has been rarely used to estimate the SOCD directly; and it cannot map the SOCD for a certain soil profile. Hence, imaging spectroscopy, as a new technique with high spectral and spatial resolution, has become an effective method for mapping soil properties mostly for profiles with fine spatial resolutions [21], [22], [23]. As a promising soil analysis tool [24], imaging spectroscopy can also effectively compensate for the disadvantages arose the fixed-depth sampling and obtain continuous depth functions for soil properties [25]. Lately, imaging spectroscopy has been successfully applied to a soil property analysis [21], [26] and a profile mapping system [23], [27], [28]. For instance, This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Sorenson et al. [29] and Xu et al. [30] mapped SOC and carbon components for soil profile, and extracted their vertical distribution patterns. Their mapping of SOC suggests that imaging spectroscopy helps capture both vertical and lateral variations of soil properties more quickly and efficiently. Today, the advancement of mapping technology has become an essential component of the digital soil morphometrics, and widely believed that it would help estimate, and map SOCD profiles including the SOC stock time-and cost-more efficiently [31].
Apart from mapping, carbon contents and properties in soil spectra are modelled using various calibration methods, which include linear models, such as multiple linear regression, principal components regression, and partial least squares regression (PLSR) [30], [32], [33], [34], [35]. In addition, nonlinear models, such as support vector machine (SVM), random forests (RF), and artificial neural networks, are often used to predict soil properties [32], [36], [37], [38]. Despite having the development of robust modelling, many scholars have explored biases in the accuracy of model predictions in SOC. For example, Hobley et al. [39] used laboratory imaging spectroscopy and various machine learning methods to predict the distribution of arable land SOC, and found that the RF method relatively produce good performance than other methods. Xu et al. [32] compared the accuracy of PLSR and various machine learning methods for predicting SOC and carbon components and found that SVM generated the best prediction result. However, the optimal model suitable for predicting soil properties varies under different soil formation conditions and requires the development of a specific prediction model for each soil property for the selected study area [13], [30], [32], [40].
The coastal soil in the north of Jiangsu Province of China developed from the sediment of the Yellow and Yangtze Rivers over the past 1000 years and has typical soil formation from the parent materials. Many scholars have studied the centennial or millennial scale dynamics of soil organic carbon in this region using traditional methods [9], [41], [42], yet there is limited research evaluating the potential of imaging spectroscopy technology for investigating soil evolution processes. The objectives of this article are as follows: 1) establish prediction models of optimal SOC content and density by using SPXY sample set partitioning methods combined with PLSR, SVM, and RF model; 2) directly and indirectly predict the profile SOCD by using imaging spectroscopy; 3) compare vertical distribution of SOCD and evolution of the SOCs between predicted and measured values at the millennium scale coastline change in the region.

A. Field Sample Collection and Processing
The study area was in the coastal plain of Dongtai, Jiangsu Province, which has a subtropical monsoonal climate (mean annual temperature 14.6°C; mean annual precipitation 1042.3 mm). The soil, east of the Fangong Dike, is saline soil, which developed from the sediments transported by the Yellow and Yangtze Rivers over the past 1000 years. Historical changes of the coastline are shown in Fig. 1(a) [9], [43]. Ten soil pits were nearly vertical to the historical coastline from east to west. The age of soil formation in each profile can be estimated according to the history of the coastline. The detailed information on the sampling sites can be found in Supplementary Table I. The intact profiles were sampled and stored in rectangular wooden boxes (100 cm in length, 20 cm in width, and 5 cm in thickness) [see Fig. 1(b)]. In each soil pit, 10 layer samples were collected for the SOC determination at fixed depth interval: 0-5 cm, 5-10 cm, 10-15 cm, 15-20 cm, 20-30 cm, 30-40 cm, 40-50 cm, 50-60 cm, 60-80 cm, and 80-100 cm [see Fig. 1(c)]. A total of 100 soil samples were obtained from 10 sampling sites. The samples for bulk density (BD) determination were collected from the center of each depth interval using a steel cylinder (100 cm 3 volume).

1) Image Spectroscopy Measurement:
The image spectroscopy for the profiles were recorded using an INFINITY V10E imaging spectrometer (Lumenera, Canada) in a dark room. The wavelength ranged between 389 and 1045 nm, and 256 bands were collected. The collection steps were performed in the laboratory: four halogen lamps were used as the light source, the air-dried soil profiles were placed on a mobile platform, and a wooden ruler of 1 m was placed on the side of the wooden box [see Fig. 2(b)]. The profile images were acquired after black and white corrections.
2) Identification and Removal of Invalid Pixels: The images were geometrically corrected to remove wooden borders and to retain the original images with a fixed width (400 columns) [see Fig. 2(c)]. Soil profiles were air-dried because the soil water content may affect the model performance [24], [44], [45]. However, the loss of water in the profile can lead to cracks [see Fig. 3(a)], which will cause errors. The spectroscopy curves of soil and crack are obviously distinct, as shown in Fig. 3(b). The soil reflectance increases significantly between 389-600 nm, and then rises with slight fluctuations. However, the reflectance for crack or shadow pixels is generally low. Therefore, the pixels with reflectance below 0.15 at 600 nm were defined as invalid pixels. This reflectance threshold was fine-tuned according to the interprofile variability to classify the invalid pixels.
3) Computation of Average Spectral Curves for the Layer Samples: The average spectral curves were extracted by averaging the reflectance of valid pixels in the corresponding depth interrval (0-5 cm, 5-10 cm, 10-15 cm, 15-20 cm, 20-30 cm, 30-40 cm, 40-50 cm, 50-60 cm, 60-80 cm, and 80-100 cm). These averaged curves were used as the original spectral data for corresponding layer samples [see Fig. 2 4) Pre-Processing of Spectral Curves: In the process of spectral data collection, certain noise may be generated because of the laboratory environment and the instrument. Then, spectral curve of each pixel in profile image was smoothed and the average spectral curves were smoothed using the Savitzky-Golay method (three orders of 25 points) to remove noise. The bands  at both ends were removed. We retained the data for 230 bands between the range of 410-1000 nm.

5) Division for Calibration and Validation Sets:
The SPXY method was used to select the calibration and validation sets based on the independent variable (spectral information) while considering the Euclidean distance from the dependent variable (soil properties) [46]. The distance equations are as follows: where x p , x q represent two different samples; N represents the number of features of the sample. The distance calculated using (6) was used to select calibration and validation samples. In this article, the ratio of the calibration and validation sets was 7:3 [see Fig. 2(e)].

C. Chemical Analysis
After air-drying and ground, the soil samples were sieved by a 0.15 mm mesh sieve. SOC content was determined using the potassium dichromate oxidation-external heating method. The samples for bulk density were weighed after drying in an oven for 48 h at 105°C. The following formula was used to calculate the BD: where m t and m s are the mass of the steel cylinder with the drying soil sample and mass of the empty steel cylinder (g), respectively, and v is the volume of the steel cylinder (100 cm 3 ). The soil organic carbon density per centimeter (SOCD icm , g m −2 ) of layer i refers to the SOC storage of 1 cm interval per unit area. The calculation formula is shown in SOCD i for soil layer i was calculated using The SOC stock (SOCs i ) from surface to certain depth was calculated using where i is the number of soil layer; SOC i , BD i , D i , and C i represent the SOC content (g kg −1 ), soil BD (g cm −3 ), soil depth (cm), and the volume (%) occupied by gravel larger than 2 mm, respectively; A represents 1 cm; C i is negligible because there are no coarse particles in this article area; and 10 is the unit conversion factor.

1) Indirect and Direct Prediction of SOCD:
Indirect prediction of SOCD involves predicting the organic carbon content using imaging spectroscopy, and then combining it with soil bulk density data to calculate the SOCD. Direct prediction of SOCD refers to directly establishing a regression model between imaging spectroscopy and SOCD.D.B. Predictive models and model evaluation.
The calibrations were built with PLSR, SVM, and RF [see Fig. 2(e)]. PLSR was proposed by Wold et al. [47]. This method combines typical correlation analysis, principal component analysis, and multiple linear regression analysis, effectively solving the multicollinearity problem among independent variables. Its discriminant equation is as follows: where i represents the current number of principal components; P is a preset value usually set to 0.05; If the difference between R 2 i+1 and R 2 i is less than P when adding a new principal component to the regression model, only the first i principal components are selected to build the model. SVM can generalize unknown samples for applications such as identification classification, regression modeling, and time series prediction [48]. In this article, the grid search method was used to search the model parameters c and g for the Gaussian radial basis function, and the 10-flod cross-validation was used to determine the best model parameters to minimize the rootmean-squared error (RMSE).
RF can improve prediction accuracy without significantly increasing the amount of computation and is one of the best data mining algorithms available [39]. In this article, the mtry parameter was set to its default value, and the best value of the parameter ntree was selected through 10-fold cross-validation and grid search method.
In this article, three metrics were used to compare the model accuracy: determination coefficient (R 2 ), RMSE, and RPIQ (RPIQ is the ratio of IQ to RMSEp, IQ is the difference between the third quartile (Q3), and the first quartile (Q1) of the sample observations). The higher R 2 and RPIQ and lower RMSE, the better model is [49], [50]. RPIQ is a more objective index than RPD for the model evaluation of non-normally distributed soil data [49], [50]. The subscript C represents the modeling results of the calibration set. The subscript P represents the independent validation results.
The best model was used to map SOC or SOCD for profiles [see

A. Characteristics of Soil Properties 1) Identification of SOC Variation:
The soil properties are shown in Table I. The SOC content of coastal soil in northern Jiangsu was recorded generally low. The mean value, was 4.41 g kg −1 with a large difference in range (1-17.03 g kg −1 ). The coefficient of variation (CV) was as high as 77.21%. The BD, (1.09-1.62 g cm −3 ) had the standard deviation and CV were 0.11 and 7.75%, respectively. The soil organic carbon density (SOCD icm ) ranged from 0.13 to 1.86 g m −2 with a CV of 70.85%. Among 10 soil profiles studied, the highest and the lowest SOCs were 59.47 g m −2 and 23.88 g m −2 . Kurtosis can describe the steepness of the data distribution pattern, with three representing a normal distribution and 1.8 representing a uniform distribution. Skewness measures the degree of skewness of the data distribution, with zero representing a normal distribution. None of the soil properties satisfied the normal distribution.

2) Variation of SOC in Soil Profile:
The SOC content in the coastal Jiangsu was found to be varying with depth, but in certainty and in a regular basis [see Fig. 4(a)]. Except for DT09, the SOC content in the profiles decreased with depth, and tended to be gentle in the deep layer. The anomaly in DT09 may be related to the buried peat layer. In addition, the vertical distribution of SOC in the profiles with different soil formation ages differed significantly. The profiles with a short soil formation age exhibited low SOC accumulation and a relatively uniform vertical distribution.
With an increase in soil age, the accumulation of SOC gradually extended to subsurface layers. The SOC accumulation in profile DT10 was more significant than that in the other profiles.
The vertical distribution of SOCD icm showed the same trend as SOC [see Fig. 4(b)]. The SOCD icm decreased along with depth. In general, SOCD icm in older profiles were higher than those in younger profiles for the upper 40 cm. For the entire profile, the average SOCD icm ranged from 0.27 g m −2 (DT02)  to 0.60 g m −2 (DT10) and generally increased with soil age. The power and logarithmic functions showed comparable results, better than linear function, which indicated the decreasing rate of SOC sequestration in this chronosequence [see Fig. 4(c)].

B. Prediction of SOC Density by Spectroscopy 1) Indirect Prediction of SOC Density by Spectroscopy:
Imaging spectroscopy was commonly used for estimating SOC content with high estimation accuracy, and SOC content can be converted to SOCD data by BD [23], [27], [29], [30]. Therefore, we first established a SOC estimation model based on averaged spectra and mapped the SOC content for profiles using imaging spectroscopy, and then obtained the SOCD map combined with the corresponding BD for layer samples.
a) Statistical analysis of SOC content in sample sets: Table II showed the SOC content statistics of the sample set divided by the SPXY method. The mean and standard deviation of the SOC content in the model set were 4.84 g kg −1 and 3.55 g kg −1 , respectively. The kurtosis and skewness were 1.60 and 1.10, which had only had a little difference from that of the total sample set. The content of the SOC in the prediction set varied from 1.09 g kg −1 to 9.20 g kg −1 , and the mean and standard deviation decreased to 3.40 g kg −1 and 2.85 g kg −1 , respectively. The spatial heterogeneity increased, and the CV was 83.92%. The kurtosis was -0.95, indicating that the data distribution was relatively flat compared with the normal distribution.
b) Modeling results for SOC content: Three methods (PLSR, SVM, and RF) were used to establish SOC content prediction model, and the prediction accuracy was compared (see Fig. 5). The results demonstrated that the R 2 c of RF was the highest at 0.93 but that R 2 p was low (0.80). Comparatively, the R 2 c of SVM and PLSR were second to RF, 0.88 and 0.73, respectively. But the R 2 p of SVM was higher than that of RF, at 0.92. The RPIQ of SVM was higher than that of the PLSR and RF, indicating that SVM was more suitable to map SOC content for profiles.

c) Mapping results of SOC content in profiles:
The SVM model was used to map the SOC content for the profiles, and the results were shown in Fig. 6. There were significant differences in the vertical distribution of SOC in profiles of different age. The SOC contents in DT01 and DT02 were relatively low with mean value 4 g kg −1 . The low carbon content could be due to the shorter time of soil formation, especially in DT02 profile. In DT03, DT04, DT05, DT06, DT07, and DT08, the SOC content in the surface layer (30 cm) ranged from 4-8 g kg −1 . The SOC content in DT09 at 0-40 cm layer ranged from 4-8 g kg −1 while in DT10 at 0-15 cm layer was measured above 8 g kg −1 and at 15-40 cm layer, ranged from 4-8 g kg −1 . Interestingly, the SOC content in the 40-100 cm layer was less than 4 g kg −1 in almost all profiles. The SOC contents in DT01, DT03, and DT04 increased slightly at 90-100 cm depth, because of the organic residues buried during deposition. The mapping results showed similar vertical distribution and temporal trend of the SOC content as the measured values presented in Fig. 4(a). It is s clear that the SOC was found to be more concentrated in the surface layer, decreased along with depth and eventually stabilized for ageing. d) Indirect mapping of SOCD icm for profiles: The SOCD icm distribution was calculated by combining the SOC map with the corresponding BD (see Fig. 7). As shown in Figs. 6 and 7, the vertical distribution of the SOC content and density in the profiles were roughly similar. The SOCD icm in the younger  profiles, DT01 and DT02, was low and relatively consistent with, less than 0.50 g m −2 . The SOCD icm content at 0-25 cm depth of the older (formed during 100 to 500 years ago) profiles such as DT03, DT04, DT05, DT06, DT07, and DT08, the was relatively higher than the younger profiles. The DT09 and DT10 soil profiles of 600 and 940 years old, where the SOCD icm content was measured highest at 0-35 cm depth. Except for DT01 and DT02, the vertical distribution of SOCD icm varied significantly with depth, and most of the SOCD icm at 35-100 cm depth was less than 0.5 g m −2 .
2) Direct Prediction of Soil Organic Carbon Density by Spectroscopy: a) Correlation between SOC content and SOC density: Fig. 8 shows the correlation between SOC content and SOCD icm . The fitting results for the linear and power functions were R 2 values of 0.98 and 0.99, respectively. It indicated that there was a strong correlation between SOC content and SOCD icm . In theory, SOCD icm could be predicted directly using imaging spectral data. b) Statistical analysis of SOCD icm in sample sets: The same sample sets as before were used for calibration and prediction sets (see Table III). The mean and standard deviation of the modeling set were 0.64 g m −2 and 0.40 g m −2 , respectively, and the kurtosis and skewness were -0.58 and 0.42, which had a small difference with the total set. The CV was 62.55%, which was smaller than the total set. The SOCD icm of the prediction set varied from 0.13 to 1.80 g m −2 , with an average value of 0.50 g m −2 and a CV as high as 93.27%. The kurtosis and skewness were 0.49 and 1.23, indicating that the data was skewed to the right compared to the mean value.
c) Modeling results of SOC density: The smoothed spectral data and SOCD icm were used to build PLSR, SVM, and RF models and the prediction performance was compared (see Fig. 9). The modeling and prediction performance of SVM was higher than those of PLSR and RF, up to 0.94 and 0.93, respectively. The RPIQ for SVM was 5.96, which was much higher than those of the other two models. The SVM got the best model for this article, which was suitable for mapping the SOCD icm for the profiles. d) Mapping results of SOCD icm in profiles: The SVM model was used to map the SOCD icm for the profiles directly (see Fig. 10). The SOCD icm of the invalid pixel was expressed as the average of the nearby values. The results (see Fig. 10) were similar with the maps predicted indirectly by SOC and BD (see Fig. 7). The SOCD icm was low and there was no significant difference in the vertical distribution in younger profiles (DT01 and DT02). In other profiles, the SOCD icm of the surface and subsurface soil were significantly different. These SOCD icm maps also showed that the vertical distribution of SOCD icm in profiles with different ages differed significantly. And the amount and accumulation depth of SOC increased with age.

3) Evaluation of Mapping Results:
To evaluate the mapping accuracy for SOC and SOCD icm , the maps for SOC and SOCD icm were averaged according to the sampling depths. Then, the averaged mapping values for the 100 layer samples were compared with the measured values. The results showed that the R 2 and RMSE of the predicted SOC content reached 0.77 and 2.07 g kg −1 , and the RPIQ was 2.79, indicating that the model produced a good mapping result for SOC [see Fig. 11(a)]. The R 2 for indirect prediction of SOCD icm was 0.68, RMSE was 0.27 g m −2 , and RPIQ was 2.85, respectively, indicating that the indirect mapping of SOCD icm by the model was generally effective [see Fig. 11(b)]. The mapping result for direct prediction of SOCD icm (R 2 = 0.78, RPIQ = 2.83) was slightly better compared to the indirect prediction [see Fig. 11(c)]. These three scatter plots also showed a phenomenon: the high value tend  to be underestimated and low value tend to be overestimated. As a result, the vertical variation for SOC and SOCD icm in the profiles was underestimated.
The SOCs i is the accumulated SOCD i from surface to certain depth (ith layer), representing the amount of SOC stock within certain depth in the profile. The R 2 was 0.84 and the RPIQ was 2.45, indicating that the map for SOCD icm can be used to estimate SOC stock for the profiles with better accuracy [see Fig. 11(d)]. The SOCs i were underestimated within 40 cm and then the values were overestimated beyond 40 cm due to the phenomenon of overestimation for low values and the underestimation for high values. That is, the underestimation for high SOCD icm led to the underestimation within 40 cm soil. The overestimation for low SOCD icm and the greater thickness of subsoil led to the overestimation for SOCS i for deeper soil and the whole profile.

C. Temporal Function of SOCs Accumulation
According to the relationship between soil formation age and SOCs at different depths, the evolution of SOCD in this region can be inferred on a millennium scale (7). The scatter diagrams between the measured and predicted SOCs at different depths and soil ages are shown in Fig. 12. The R 2 values for the linear, logarithmic, and power functions were used to determine the best fitting results for SOCs and soil age. In general, the measured SOCs were greater than the predicted values for surface soil within 20 cm, and the predicted SOCs were greater than the measured values for deeper soil (>50 cm).
The temporal function of measured SOCs varies continuously with the sampling depth. The accumulation of surface SOCs continued at a stable rate and was roughly linear on a millennial time scale. With increasing soil depth, the fitted function of SOCs with time was closed to a power function and then changed to a logarithmic function, indicating that the deep soil could not reach the organic carbon stock accumulation rate of the surface soil, thus slowing the rate of SOCs accumulation in the whole soil profile. For the predicted values of SOCs, the linear function was the best fit for all soil depths, which was a result of the overestimation of low values and underestimation of high values of SOCD from the spectral data.

A. Calibration Methods and Their Performance
Many scholars have used imaging spectroscopy to predict SOC/SOM (see Table IV). PLSR was used in earlier studies while SVM or RF has been used in recent studies. Our research indicates that the models established by SVM and RF methods are superior to PLSR, which is consistent with previous research results (see Table IV). This is mainly because machine-learning techniques can effectively handle the nonlinear relationship between spectral features and soil properties [34]. However, although these nonlinear models have higher predictive accuracy, there is a "black box" problem where the model structure is complex and the transparency is generally low [51], making it difficult to reveal the functional relationship between spectra and soil properties. It is also difficult to understand the contribution of each band to property prediction in model construction. Linear models have lower predictive accuracy but simpler structures and typically have lower computational costs [13], [30], which are advantageous for identifying response characteristics between spectra and properties from large amounts of data. Therefore, the predictive performance of different models is related to their underlying mechanisms, and each method has advantages and disadvantages that must be considered based on research purposes and data availability [32].
Previous studies showed that the imaging spectroscopy could effectively map the SOC content for profiles [22], [27], [29], [30], [33], [35], [39]. However, SOCD icm needs to be calculated based on the BD [52], [53], [54], which is complicated and laborious. Therefore, a question arises whether SOCD icm can be predicted by spectral data directly with improved performance. In this article, we compared the consistency between direct and indirect predictions of SOCD icm with the measured values, and the results showed that direct prediction accuracy was better than that of indirect prediction [see Fig. 11(b) and (c)]. Three modeling methods achieved good performance in predicting SOCD icm (see Fig. 9), with SVM method obtaining the best model, with R 2 values of 0.94 and 0.93 for the calibration and prediction sets, respectively. This means that it is possible to predict SOCD icm using spectral data directly and accurately, thus filling the gap in BD data in SOC inventory calculations, which is of great significance for estimating global soil carbon stocks. The 1:1 relationship indicated that all three models exhibited both over-and underestimation of SOCD icm predictions (see Fig. 9). The reason for overestimation of high values may have been due to the limited number of soil samples with high SOCD icm concentrations, which were only found in the top 0-10 cm layer of the DT10 profile with contents exceeding 1.5 g m −2 [see Fig. 10(j)]. Underestimation of low values may have been related to the weakening of the spectral prediction ability of soil organic matter with decreasing content [55]. Zheng et al. [13] used imaging spectroscopy combined with multivariate regression to predict the content of five soil fertility properties in Cixi area, and the results showed that all six models exhibited over-and underestimation of the five soil properties, which was consistent with this article. Since the prediction of soil properties by spectroscopy is influenced by multiple factors such as spectral bands, soil characteristics, and preprocessing procedures [13], [30], [56], the underlying mechanisms of these prediction errors were not fully understood, and further research was needed for more accurate prediction. In addition, the variable soil parent material and carbon source can also affect the model's predictive performance [57], [58], [59]. The soil samples selected in this article developed in the typical parent material area of the Yellow River and Yangtze River accumulated under the interaction of land and sea in the past 1000 years, which has reference significance for studying soil chronosequence in other coastal areas. However, soils developed on different parent materials from other regions may have different applicable models and results [60], [61]. For instance, Liu et al. [60] investigated the predictive performance of three regression models for soil organic carbon content developed on different parent materials (trachyte and basalt) and found that the best model and spectral response characteristics differed. Therefore, the model we used was applicable to this region and was not considered globally applicable. It is, however, necessary to collect soil samples from different regions for SOC/SOM and type evaluations and to extend the spectral prediction model to a larger range.

B. Variation of SOCD in Soil Profile
Many scholars have already conducted research on the vertical variation of SOCD, and their findings indicated that the surface soil had a higher SOCD compared to the lower layers [9], [62], [63], [64]. In this article, we also observed a similar distribution pattern of SOCD [see Fig. 4(a)]. The reason may be that the decomposition of plant leaves and roots formed organic matter that accumulated rapidly in the surface layer, leading to a significant increase in the surface SOCD [5], [65]. Imaging spectroscopy is an effective new method for mapping soil properties, which allows us to observe the vertical distribution characteristics of soil properties in more detail compared to traditional reflectance spectroscopy [24], [30]. To evaluate and validate the prediction results, we extracted the continuous depth variation curve of SOCD from the soil profile mapping results, and displayed the vertical distribution of measured values on the right side of the profile map (see Fig. 10). The continuous depth variation trend of SOCD in the topsoil (0-40 cm) was relatively close to the measured values, while in the deeper layers (below 40 cm), it was much higher than the measured values. The main reason for this was that the model had poor accuracy in estimating the SOCD for lower contents (less than 0.5 g m −2 ). The mapping results of soil organic matter content in the soil profile by Zheng et al. [13] showed that the best model overestimated the organic matter content in soil samples with values lower than 5 g kg −1 to varying degrees, consistent with this article. These results indicated that when mapping for SOC/SOCD, the mapping results in the low-content areas may have been unreliable. Although many studies have conducted high-resolution mapping of various soil properties in both vertical and horizontal directions, all the validation results of their established models have exhibited overestimation of high values and underestimation of low values [13], [30], [32], [40]. Further research is needed to determine the accuracy in estimating the content of attribute in high and low value regions. In addition, the generalized ability of the model may also affect the mapping effect. Steffens and Buddenbaum [23] pointed out that a high-precision regression model does not always represent the best mapping results. In their study, the accuracy of PLSR was the lowest, but its generalized potential for reasonable mapping of profile properties was the highest, while SVM had the highest accuracy but the poorest robustness and lower generalization potential.
This article filled the gap in the application of imaging spectroscopy technology in soil chronosequence research by mapping SOCD in soil profiles at different ages. The results suggest that imaging spectroscopy technology could be used to directly map the topsoil (0-40 cm) SOCD of older profiles (DT03-DT10), but the predicted accuracy of younger profiles (DT01, DT02) and lower SOCD areas (below 40 cm) needed to be improved. Nonetheless, imaging spectroscopy technology provide a fast and non-destructive method for observing the spatial structure of soil properties in profiles [30]. Compared to the limitations of traditional fixed-depth sampling intervals, this technology could better serve the study of profile horizon division, dynamic changes in soil properties, and achieve better soil survey tasks [13], [23], [30].

C. Temporal Variation of Organic Carbon Density
Soil chronosequence is of greater significance to explore the direction and rate of overall soil evolution, providing essential information for soil genetic theory [66], [67], [68]. In this article, power functions and logarithmic functions were found to be better fit the changes in 1 cm average SOCD over time compared to linear functions [see Fig. 4(c)], indicating that the soil carbon sequestration rate gradually slowed down over time for the entire profile. This was mainly because the deep soil could not accumulate organic carbon as efficiently as the topsoil, leading to a decrease in the accumulation rate of organic carbon in the entire depth of the soil [see Fig. 4(b)]. The best-fit function for the topsoil (0-15 cm) was linear, but with increasing soil depth, it became power functions and logarithmic functions (see Fig. 12), which also supported this view. However, our results were different from that of other soil evolution studies in which the accumulation rate of SOM with time at different depths all showed a linear growth trend [13]. The main reason may be that soil with different parent materials had different carbon sequestration capacity [69].
For the predicted values of SOCs, the linear function was the best fit for the 0-15 cm topsoil, which was consistent with the measured results (see Fig. 12). However, the best fitting function for predicted values at other depths was different from the measured results, with the best fitting function for predicted values being linear and for measured values being power function. Traditional sampling methods were commonly used in previous studies to analyze the direction and rate of soil evolution [69], [70], [71], [72], making it difficult to directly compare the results of the best fitting functions obtained from imaging spectroscopy with those of other studies. Nevertheless, the fitting trends of SOCs over time showed minimal differences among the three functions for both measured and predicted values. This suggests that imaging spectroscopy technology holds great potential for quantitatively revealing the process of soil evolution. Despite the current limitations in accuracy, this article provided a new insight for the rapid monitoring of soil property changes under global climate change scenarior. Future research should focus on exploring more intelligent prediction models [73], [74], [75], improve the generalized ability of the models, and reducing the extent of overestimation and underestimation of high and low values, respectively. This will enable more accurate soil property mapping to support strategic planning for regional and even national-level sustainable development goals related to soil conservation and management.

D. Limitations and Perspectives of This Article
The results of this article indicated that utilizing imaging spectroscopy data for direct or indirect estimation of SOCD distribution in soil profiles was a feasible option for studying soil evolution processes. However, there were still some limitations to this approach. First, the soil samples used for chemical analysis were grounded, which may have reduced the influence of environmental factors. However, the surface of intact soil profiles was much rougher and may decrease the accuracy of estimating SOC/SOCD. Therefore, it is necessary to develop algorithms to eliminate the light scattering effect caused by soil surface roughness. Second, the developed model for estimating SOCD and the fitted soil chronosequence function in this article are only applicable to the studied region. Under different environmental conditions, the accumulation rate and vertical distribution of SOC may present varying results. Therefore, it is necessary to compare the results obtained from applying imaging spectroscopy technology to study soil evolution to the other region.
In the future research, it is desired to collect soil profile samples under differential parent material conditions to construct soil chronosequence and use more advanced imaging spectrometers to detect more wavelengths to enrich soil information. Additionally, it may be possible to extend the developed model in the article to larger spatiotemporal scales. These can provide a theoretical basis for the soil evolution and pedogenesis.

V. CONCLUSION
This article utilized PLSR, SVM, and RF methods to predict the SOCD directly or indirectly in the entire depth (0-100 cm) of soil profiles at different ages. The results showed that the accuracy of direct prediction was better than that of indirect prediction, and the SVM method produced the best model, with R 2 values of 0.94 and 0.93 for the calibration set and validation set, respectively. The variation of measured SOCD in soil profile showed that SOC mainly accumulates in the topsoil (0-40 cm), with little accumulation in the deeper soil (below 40 cm). The mapping of SOCD achieved good results in the 0-40 cm layer of older profiles (DT03-DT10), but the prediction accuracy in younger profiles (DT01, DT02) and deep soil layers needed to be further improved. Power functions and logarithmic functions were found to be better fit the changes in 1 cm average SOCD over time compared to linear functions, indicating that the soil carbon sequestration rate in this area gradually slowed down over time for the entire profile. For the predicted values of SOCs, the best temporal functions were different from that of the measured values at most depths, but the fitting trends of changes with time by three functions were basically consistent. It suggests that the IS technology had the potential to quantitatively reveal the process of coastal soil evolution. Future research should focus on exploring more intelligent prediction models, and improve the existing generalization in the ability and prediction accuracy of the models.
Xianli Xie received the Ph.D. degree in geography from the Nanjing Normal University, Nanjing, China, in 2004.
She is currently an Associate Researcher with the Nanjing Institute of Soil Science, Chinese Academy of Sciences, Beijing, China. Her main research interest includes soil remote sensing.
Rong Zeng received the Ph.D. degree in resource environment and remote sensing information from the Nanjing Institute of Soil Science, Chinese Academy of Sciences, Beijing, China, in 2018.
She is currently a Lecturer with the School of Geographic Sciences, Nanjing University of Information Science and Technology, Nanjing, China. Her main research interests include soil hyperspectral remote sensing and digital soil mapping.
Chengyi Zhao received the Ph.D. degree in soil science from the China Agricultural University, Beijing, China, in 2002.
He is currently a Professor with the School of Geographic Sciences, Nanjing University of Information Science and Technology, Nanjing, China. His main research interests include intelligent remote sensing and artificial intelligence technology research and development, geographical Big Data and cloud platform technology application, and surface process simulation.
Giri Kattel received the Ph.D. degree in geography from the University College London, London, U.K, in 2004.
He is currently a Professor with the School of Geographic Sciences, Nanjing University of Information Science and Technology, Nanjing, China. His main research interest includes global environmental change.
Ying Xiong is currently working toward the bachelor's degree in maths and applied mathematics with the Reading Academy, Nanjing University of Information Science and Technology, Nanjing, China.
Her research interests include global environmental change and soil remote sensing.
Dian Zhou is currently working toward the bachelor's degree in human geography and urban and rural planning from Nanjing University of Information Science and Technology, Nanjing, China.
Her research interests include costal area's humannatural relationship, remote sensing, and human geography.