Visible and Near-Infrared Reflectance Spectroscopy Analysis of a Coastal Soil Chronosequence

Abstract

The soil chronosequence is a useful method for investigating pedological theories. Soil chemical, physical and mineralogical properties in chronosequences change over time and exhibit systematic and time-dependent trends, which can be used to analyze the rates and directions of pedogenic changes. The potential of soil spectroscopy as an emerging, rapid and cost-effective technique for predicting soil properties has been widely accepted and has motivated the application of spectroscopic techniques to the analysis of soil chronosequence. We present a soil chronosequence derived from 1000-year-old calcareous marine sediments and examine changes in six soil properties over this period. We evaluated the utility of a soil spectroscopic method to detect soil property changes and to predict the pedogenic properties and soil ages of the chronosequence. The results show that some soil pedogenic processes, such as soil organic matter accumulation, CaCO₃ leaching and clay migration, can be identified in the millennium chronosequence. Power chronofunctions are formulated for soil organic matter (SOM) and Logarithmic chronofunctions are fitted for clay, CaCO₃ and pH. These pedogenic processes are identified in the reflectance intensity and absorption features of soil spectroscopy, and pedogenic properties can be calibrated via soil reflectance spectroscopy. Profile ages can also be predicted via pseudo multi-depth spectra of soil profiles, and soil spectral curves for 0–30 cm generated the best prediction results (RPD = 1.85). We conclude that soil properties, changing due to weathering and soil formation, act as a bridge linking spectroscopy and weathering levels/pedogenic processes. The results imply that applying spectroscopy techniques to chronosequence study and mapping the degree of soil development in certain areas should be possible.

1. Introduction

According to soil genesis theory, soil is a natural body formed through the synthetic actions of parent materials, climate, topography, organisms, and time [1]. These factors continuously determine changes in soil-forming processes, resulting in soil evolution [2]. Every soil profile records its genetic process, developmental history and duration of its exposure to soil-forming factors [3]. A series of soil profiles, developed based on similar soil-forming factors with the exception of time, consist of a soil chronosequence [2] which splices brief segments of the soil evolution history [4]. Consequently, spatial difference between the soil profiles can be translated into temporal differences [5]. Therefore, soil chronosequences can reveal the rates and directions of pedogenic changes and provide invaluable information for examining pedogenic theories and processes [5,6]. Soil properties in a chronosequence exhibit systematic and time-dependent features, i.e., soil organic matter (SOM) [7,8], iron oxide [9,10], clay minerals [11], particle size distribution [12], and CaCO₃ [13,14]. These properties are also called pedogenic properties. Therefore, these time-dependent pedogenic properties can be used for relative age dating [15]. Furthermore, chronofunctions of soil properties were developed to determine the stage and the rates of soil development, and to ascertain the time of reaching steady state [16].

The chemical and physical properties are conventionally measured through laboratory-based “wet chemistry” analyses, which are costly and time-consuming, and not suitable for rapid soil property mapping at multiple spatial scales. Whereas, considerable demand for high-quality and inexpensive soil data for environmental monitoring, modeling and precision agriculture requires more efficient and cost-effective methods of soil analysis nowadays [17]. Since Bowers and Hanks [18] and Condit [19] systematically investigated the relationship between soil properties and reflectance spectroscopy, the potential of soil spectroscopy as an efficient and cost-effective prediction technique has been demonstrated by a large number of subsequent studies [20,21,22,23]. This technique was also used in pedogenic research to explore the relationship between the degree of soil weathering and reflectance spectroscopy [24,25,26,27]. Demattê and Terra [28] investigated the reflectance of weathered soils in a toposequence and noted that soil changes, produced by different weathering levels or pedogenic processes, reflected spectral behaviors along the profile depth. Zheng et al. [29] concluded that the pedogenic properties provided a bridge between spectroscopy and the pedogenic processes, which makes it possible to analyze the chronosequence using spectroscopy. To the best of our knowledge, Awiti et al. [30] first evaluated the capability of spectroscopy to analyze soil property changes for topsoil (0–20 cm) across a forest-cropland chronosequence. However, a chronosequence that utilizes a series of deeper soil profiles with formation age, in addition to topsoils, may provide more comprehensive and accurate soil property changes and pedogenic processes.

The coastal soil is an important source of agricultural lands in China and the understanding of soil dynamics under long-term cultivation is critical to planning the sustainable use of the soil resource. Chen et al. [31] and Cui et al. [32] revealed the dynamic changes in soil properties through the chronosequence derived from the sediments of the Yangtze River. However, few investigations of the soil dynamics of chronosequence formed by the sediments of the Yellow and Yangtze Rivers over the past 1000 years have been published.

In this work, we analyze fairly uniform calcareous marine sediments with exact coastline change records in order to evaluate the utility of soil profiles’ spectra for a chronosequence. More specific objectives are (i) to establish a soil chronosequence derived from 1000-year-old calcareous marine sediments and investigate soil property changes over this period, (ii) to evaluate the potential for surface soil spectroscopy to detect soil property changes with time in a chronosequence, and (iii) to assess the feasibility of predicting soil profile age in this area using spectroscopy.

2. Materials and Methods

2.1. Study Area and Soil Formation

The study area is located in the coastal plain of Northern Jiangsu Province, China (Figure 1). The coastal plain east of the Fangong Dike, completed in 1027 AD, is formed of sediments transported mainly by the Yellow and the Yangtze Rivers over the past 1000 years. Historical changes of the coastline in Jiangsu Province are shown in Figure 1. This region is characterized by the northern subtropical monsoon climate. The mean annual precipitation and temperature are 1042.3 mm and 14.6 °C, respectively.

The soil formation and evolution in this area went through three stages. First, the seabeach, which was formed from sediments, evolved into a meadow solonchak (salt marsh) after the salt in the topsoil was leached out and salt-tolerant plants appeared. Secondly, SOM in the topsoil accumulated with the development of salt-tolerant plants, and soils further evolved into saline farming soils with the cultivation of salt-tolerant crops. Finally, dry-farming soils formed with the cultivation of corn– and rice–wheat/barley rotations.

2.2. Soil Sampling and Analysis

Sixteen soil sampling sites are located along three lines (north, middle and south) that generally run perpendicular to the coastlines (Figure 2a), and soil ages were estimated based on the historical coastlines to indicate the development degree of soil profiles [33]. The land use of all the sampling sites was farmland. All the soil profiles belong to the coastal saline alluvial soil, except profile S4, which is paddy soil according to Chinese Soil Genetic Classification System. To specify the soil property changes at various profile depths and compare the profile properties, the following sampling-depth intervals were used for each profile: 0–5, 5–10, 10–20, 20–30, 30–40, 40–60, 60–80 and 80–100 cm and 128 samples were collected. The samples for every specific depth interval were collected from the whole interval (Figure 2b) in June, 2013 and the weight for every sample is about 1–2 kg.

The soil samples were air-dried and gently crushed using a wooden pestle to pass through 2 mm nylon sieves. The subsamples were further ground using agate mortar and passed through 0.250 mm and 0.150 mm sieves to determine their physical and chemical properties. Soil pH was determined using the electrode method with a 1:5 (g:mL) ratio of soil and distilled water. SOM was measured by the potassium dichromate–external heating method, CaCO₃ content was measured through neutralization titration, and the particle size distribution was determined using the pipette method. Total Fe (Fe_t) was determined by the HF-HNO₃-HClO₄ digestion–atomic absorption spectrophotometry method, free Fe oxides (Fe_d) were extracted using the dithionite-citrate-bicarbonate (DCB) method [34]. The descriptive statistics and correlation analysis of soil properties were performed by SPSS 20.0.

2.3. Spectra Determination

The air-dried soil samples passed through a 0.150 mm sieve were used for spectral measurements. Soil reflectance spectra were measured in a dark room with a FieldSpec 3 portable spectrometer (Analytical Spectral Devices Inc., Boulder, Colorado, USA). A 50 W halogen lamp set at a vertical angle of 15° (A) and at a distance 30 cm (L) from the samples was used as a light source. The sensor was vertically positioned 15 cm (H) above the soil samples using a probe set at a 5° viewing angle (half view angle B = 2.5°), producing a radius of 0.655 cm (R2) from the center of the sample dish (radius of sample dish R1 = 2.5 cm) (Figure 3). The dark current was removed prior to collecting the soil spectra, and a 40 cm × 40 cm diffuse reflection standard reference plate was used to obtain the relative reflectance. Throughout the measurement process, the soil sample in the container was rotated on an automatic rotating platform to average out the possible influence of different azimuth angles. Twenty spectral curves were collected for each sample as the sample was turned 360°. The reflectance spectrum for each sample was deemed the average of these 20 spectral curves.

2.4. Statistical Analysis and Transformation of the Spectral Data

Partial least squares regression (PLSR), first proposed by Wold et al. [35], is a commonly used modeling technique for calibrating soil properties to reflectance spectroscopy because it can resolve the issues of high collinearity among predictor variables (reflectance spectroscopy). PLSR integrates compression and regression steps and can effectively extract comprehensive variables (referred to as latent variables, LVs) by excluding unexplainable information, thus ensuring that comprehensive variables are optimally condensed from the raw inputs to explain response variables. Although various data-mining algorithms have been developed to determine soil properties using visible and near-infrared (VNIR) diffuse reflectance spectra and gave better performance than PLSR [36,37], the coefficients of PLSR can give the important wavelength to response variable which is useful to analyze the prediction mechanism. All calibration and prediction were performed by PLSR in Unscrambler X v. 10.0 (CAMO Software AS, Oslo, Norway). The data were centered before calibration and prediction automatically by checking the center data tick-boxes in this software. Full cross-validation was used for the cross-validation method in the calibration process. The optimal number of principal components was determined by the residual variance obtained in the unscrambler.

A continuum removal (CR) technique can be used to isolate particular absorption features for spectrum analysis and to define the band depth [38,39]. In this work, spectral curves were smoothed following Savitzky and Golay [40], with a second-order polynomial across 21 smoothing points using Unscrambler X v. 10.0 (CAMO Software AS, Oslo, Norway). These curves were then transformed via continuum removal in ENVI 4.7 (ITT Visual Information Solutions Inc., Broomfield, United States). For each soil profile, the smoothed and continuum-removed spectra from the eight depth intervals were joined in sequence to create a pseudo multi-depth soil spectral curve (Figure 4) according to Vasques et al. [41].

2.5. Datasets Partitioning and Accuracy Evaluating

All the samples were sorted according to profile number and depth in a profile. Four out of five samples were used for model development (calibration) and the others (one out of five) for model evaluation (prediction). The number of samples was 102 and 26 for the calibration and prediction (evaluation) sets, respectively.

Two principal components were extracted from the pseudo multi-depth curves of continuum-removed spectra of the sixteen profiles. Six profiles were selected for the prediction of soil age according to a scatter plot of two principal components and others were used for calibration of the profile age.

The model accuracies were evaluated using the coefficient of determination (R², subscripts c for calibration, p for prediction) and residual prediction deviation (RPD). Viscarra Rossel et al. [17] classified RPD values into six classes: RPD < 1.0 indicates very poor results which are not recommended; RPD between 1.0 and 1.4 indicates poor results where only high and low values are distinguishable; RPD between 1.4 and 1.8 indicates fair results which may be used; RPD between 1.8 and 2.0 indicates good results where quantitative predictions are possible; RPD values between 2.0 and 2.5 indicate very good, quantitative model/predictions, and RPD > 2.5 indicates excellent model/predictions.

3. Results

3.1. Basic Statistics of Soil Properties

The descriptive statistics of soil properties are presented in

Table 1. Basic statistics for six soil properties.

Soil Property	Min	Max	Mean	SD	CV (%)	Skew	Kurt
SOM	1.4	33.4	7.1	6.2	86.9	1.5	2.5
Clay	1.9	18.3	7.7	4.0	51.3	0.7	−0.3
CaCO3	0.9	63.8	30.8	16.8	54.4	0.0	−1.0
pH	6.4	9.2	8.5	0.4	5.1	−1.7	5.3
Fed	1.8	6.6	3.3	0.7	20.5	0.9	3.7
Fed/Fet	7.8	26.8	15.0	3.3	21.7	0.8	2.0

Note: Min: minimum; Max: maximum; SD: standard deviation; CV: coefficient of variation; Skew: skewness; Kurt: kurtosis; soil organic matter (SOM): g kg⁻¹; Clay: %; CaCO₃: g kg⁻¹; Fe_d: g kg⁻¹; Fe_d/Fe_t: %; Hereinafter inclusive.

. SOM were mainly distributed between 1.4 g kg⁻¹ and 33.4 g kg⁻¹ with a mean of 7.1 g kg⁻¹, indicating a low SOM level in this area. Clay varied from 1.9% to 18.3% with a mean of 7.7% and a standard deviation of 4.0. CaCO₃ ranged from 0.9 g kg⁻¹ to 63.8 g kg⁻¹ with a mean of 30.8 g kg⁻¹. The value of pH varied from 6.4 to 9.2 with a mean of 8.5, indicating alkalinity because of CaCO₃ in soil. Fe_d and Fe_d/Fe_t ranged from 1.8 g kg⁻¹ to 6.6 g kg⁻¹ and from 7.8% to 26.8%, respectively. SOM was the most variable property among these soil properties (CV = 86.9%), followed by CaCO₃ (CV = 54.4%).

Table 2. Correlation coefficients of six soil properties.

	SOM	Clay	CaCO3	pH	Fed	Fed/Fet
SOM	1.00	0.75 **	−0.64 **	−0.48 **	−0.08	−0.03
Clay		1.00	−0.76 **	−0.41 **	0.20 *	0.21 *
CaCO3			1.00	0.39 **	−0.28 **	−0.33 **
pH				1.00	0.09	0.10
Fed					1.00	0.89 **
Fed/Fet						1.00

Note: ** Significant correlation at the 0.01 level; * Significant correlation at the 0.05 level.

lists the correlation coefficients of six soil properties. SOM had significant positive correlation with clay, while negative correlations with CaCO₃ and pH. There were significant negative correlations between CaCO₃ and clay, Fe_d and Fe_t, while there was a positive correlation with pH.

3.2. Soil Property Changes in Chronosequence and Pedogenic Interpretations

SOM, clay, pH, and CaCO₃ exhibited direct dependences on age, and the vertical distribution of these properties changed regularly with depth (Figure 5). Trends over time and the vertical distribution of Fe_d were not evident. The changes in the properties in the top surface (0–5 cm) were fitted to analyze the chronofunctions. The four functions, including linear, power, exponential and logarithmic models, were used to fit the scatter, and the best fitted results were determined by the coefficients of determination and shown in Figure 6.

SOM increased with soil age between profiles and decreased with depth in a profile (Figure 5a). The vertical distribution of SOM in younger profiles (M1 and N1) is relatively uniform throughout the 0–100 cm soil depth. However, the differentiation in the vertical distribution of SOM in the profiles increased with soil age, i.e., older soils presented greater SOM differentiation between the topsoil and subsoil. In other words, the shapes of depth function of SOM changed gradually from uniform to exponential with soil age [42]. In addition, SOM was similar in the subsoil under a depth of 40 cm, except in the M3 and N4 profiles. The results suggested that SOM accumulation mainly occurred in the top 40 cm due to cultivation and biological activity of salt-tolerant plants. Brauer et al. [43] showed that paddy rice cultivation will lead to a dense plough pan and seriously reduces, but not totally prevents, downward transport of organic matter.

A scatter was plotted with SOM in the top soil sample (0–5 cm) as dependent and soil age as independent (Figure 6a). The functions, fitting to the scatter, were named chronofunctions, which were formulated relating the pedogenic properties to the ages of the soils [44]. A logarithmic model for SOM is common because it suggests that the SOM will approach a steady state [16]. However, power function got the best result (largest R²) among linear, logarithmic and power functions in this study (Figure 6a), which suggests that SOM in this area has not reached a steady state. The soil in this area can continue to play the role of carbon sequestration. This curve still displayed the general pattern of rapid increase in the first period and then a lower rate of accumulation [32,44,45]. The rate of SOM accumulation is different in many studies. Chen et al. [31] analyzed some studies about the time taken to reach a steady state and concluded that different environmental conditions may produce different response times. However, many studies proved that positive anthropogenic activities, especially rice–barley rotation and fertilization [6,7], will enhance the accumulation of SOC. Furthermore, the mass fraction of SOC was widely used for modeling chronofunction, but we think SOC density or stock is better for SOC accumulation analysis.

The changes in clay content were similar to those in SOM (Figure 5b). Clay content in the youngest profile was uniform, and differentiation between the topsoil and subsoil increased with age because larger particles physically and chemically weathered to smaller particle sizes over time [46,47]. The vertical distribution of clay content within a depth of 0–40 cm or less in every profile was more uniform than that of SOM, suggesting that clay migration and illuviation mainly occurred at this depth during the past 1000 years. Artificial irrigation and drainage during cultivation can favor clay accumulation and differentiation [48], which may be the cause of the difference of vertical distribution between clay and SOM. The total mass of accumulated clay is an excellent indicator of relative soil age during all stages of soil development [15]. However, the relationship between clay and age is different. Bockheim [49] concluded that the logarithmic model resulted in the best fit of clay in B horizon and Huang et al. [10] got the same result. However, VandenBygaart and Protz [44] indicated that a linear chronofunction was the best model to fit the increase in clay with soil age. In this study, logarithmic chronofunction obtained the best result (Figure 6b).

The CaCO₃ content showed an inverse trend relative to SOM and clay. The CaCO₃ contents were high and remained uniform in profile M1 and N1 and decreased with soil age due to leaching (Figure 5c), consistent with the results of previous studies [13,44]. A logarithmic curve fit to the scatter plot yielded a high determination coefficient (Figure 6c). In this study, high CaCO₃ contents and rapid leaching rates appeared in the profiles of less than 600 years; low CaCO₃ contents and lower leaching rates were found in the >600-year profiles in this study area. However, Chen et al. [31] found that CaCO₃ was completely leached from paddy soil within a 300-year period. Wissing et al. [50] indicated that the upper 20 cm of soil was free of carbonates in 100-year-old paddy soils. Vidic and Lobnlk [13] discovered that carbonate clasts vanished from the upper horizons after 62 ka. Precipitation and the farming system of the study areas caused this difference, as precipitation (and resulting soil wetness) and vegetation type were the two main factors that affected the rate of carbonate weathering and the depth of leaching in carbonate-rich soils [44]. To the contrary, CaCO₃ accumulation is the most characteristic age-related pedofeature in the semiarid regions where precipitation is insufficient [51]. The content of CaCO₃ can be seen as an intrinsic threshold [52] because the changes in CaCO₃ may be responsible for many successive processes. Firstly, the leaching of carbonates in soils is commonly a prerequisite to clay migration by eluviation. Secondly, Fe_d illuviation occurred only when carbonates were eventually leached because CaCO₃ will retard the transformation of silicate Fe to Fe oxide [31]. There are significant negative correlations between CaCO₃ and clay and Fe_d in this study (Table 2), which suggests that the leaching of CaCO₃ is very important to clay and Fe_d.

Generally, the soil pH value decreased with soil age, following changes in the CaCO₃ content (Figure 5d and Figure 6d). After about 1000 years of development, the soil pH decreased from >8.5 to nearly 7.5. The relationship between pH and CaCO₃ is significant and positive (Table 2). Therefore, the best-fit curve of pH scatter was logarithmic, the same as that of CaCO₃ (Figure 6d). There were significant negative relationships between pH and SOM and Clay (Table 2). Vidic and Lobnlk [13] demonstrated that carbonate leaching and the corresponding decrease in pH created favorable conditions for clay translocation.

In contrast with the above-listed properties, Fe_d and Fe_d/Fe_t did not exhibit evident trends correlated with soil age. The vertical distributions of these two properties were relatively uniform in each profile, compared with the above four properties (Figure 5e,f). The fitted chronofunctions of these two properties still showed positive linear relationships with soil age (Figure 6e,f). In general, Fe_d and Fe_d/Fe_t increased with age due to the irreversible transformation degree from primary iron silicate minerals to secondary iron oxide [32,53].

3.3. Spectroscopy Analysis of Top Soil

Soil spectroscopy is a cumulative property derived from the intrinsic spectral behavior of a heterogeneous combination of soil constituents [54]. SOM, clay minerals, carbonate, iron oxides, moisture, and particle size are important factors for reflectance [55]. These soil properties simultaneously determine the reflected energy [24] and form spectral curves. Therefore, changes in the soil properties can be studied based on diagnostic absorption features [56].

In order to describe the changes in the soil spectroscopy clearly, these profiles were divided into three groups, i.e., north line, middle line and south line, and each group is perpendicular to historical coastlines (Figure 2a). The original and continuum-removed spectral curves of topsoil were plotted in Figure 7. In the range of 350–1300 nm, the original reflectance density increased rapidly and the curves were convex with some broad absorption features by electronic processes of iron (Figure 7a,c,e). In the range, 1300–2500 nm, the original reflectance curves were almost parallel with some sharp absorption features by vibrational processes. The absorption features were more prominent after transformation by the continuum removal algorithm [38] (Figure 7b,d,f).

SOM, an important soil component, has a significant effect on soil reflectance. Generally, reflectance intensity decreases due to strong absorption of SOM and increases considerably following SOM removal [24]. In this study, the original reflectance intensity in the whole wavelength decreased with soil age (Figure 7a,c,e). This difference in intensity was mainly attributable to SOM changes, i.e., overall changes in reflectance intensity demonstrated the process of SOM accumulation with age. The reflectance intensities of N6, M4 and S4 were lowest in every group because their contents of SOM are highest in their respective group. Some exceptions existed in every group because the reflectance intensity was affected by other properties such as clay minerals, iron oxide and carbonate. This is also the reason that some spectral curves were similar, although SOM content is different (for example, S1 and S2, M3 and M4).

The iron oxides are the most active factor in the visible and near-infrared (VNIR) range because nearly all electric features in this range are produced from various forms of iron [57]. Crystalline iron was responsible for concavities in the 450, 850 and 1100 nm, but amorphous iron only reduced the intensity of reflectance and did not alter the concavities [58]. Goethite and hematite are the most well-known iron oxides and have respective characteristic absorption features (Figure 8). Ben-Dor et al. [25] indicated that the position of the absorption feature of goethite and hematite mixture shifted between the original absorption of the pure mineral absorption and this shift in the Fe absorption features enabled the recognition of hematite and goethite occurrences. In this study, the presence of distinct absorption features near 420, 485–490, 660–670 and 900–910 nm in the spectral curves following continuum-removed spectra indicated the presence of more goethite and little hematite in these soil samples. However, these absorption features of Fe did not present evident changes with soil age.

Hunt and Salisbury [59] concluded that five characteristic bands appear near 1900, 2000, 2160, 2350 and 2550 nm for carbonates in the NIR range. The absorption feature of CaCO₃ near 2341 nm can be used to estimate carbonate content via continuum-removed reflectance spectroscopy [60,61,62]. In this study, the absorption depth at 2345 nm decreased from younger soil to older soil with a decrease in CaCO₃ in every group (Figure 7b,d,f). This absorption depth at 2345 nm was positively correlated with CaCO₃, i.e., the absorption depth increased linearly with increasing CaCO₃ content (Figure 9a). In addition, the change in the absorption depth at 2345 nm with soil age is distinct and a logarithm function fitted the relationship between the absorption depth and soil age (Adj. R² = 0.7103) (Figure 9b). Therefore, we can conclude that the spectral feature near 2345 nm can be used to estimate the CaCO₃ content and to reveal the soil CaCO₃ leaching process in this area.

The obvious absorption features near 1412, 1918, 2210, 2345 and 2445 nm and the two shoulders near 1465 and 2095 nm in these CR curves (Figure 7b,d,f) suggested that the clay mineral in these soil samples is mainly illite. Montmorillonite has absorptions near 1412, 1918, 2210 and the two shoulders, but the absorption depth at 2210 nm is near half of that at 1918 nm (Figure 10). Therefore, montmorillonite is also present in these samples.

3.4. Soil Property Predictions

The PLSR results for the soil properties of the calibration and prediction with smoothed spectra are shown in

Table 3. Calibration and prediction of the soil properties.

Soil Property	LVs Number	Calibration (n = 102)		Cross Validation		Prediction (n = 26)
		RMSEC	${\mathbf{R}}_{\mathbf{c}}^2$	RMSECV	${\mathbf{R}}_{\mathbf{c}\mathbf{v}}^2$	RMSEP	${\mathbf{R}}_{\mathbf{p}}^2$	RPD
SOM	9	1.76	0.92	2.22	0.88	1.98	0.89	3.04
Clay	9	1.35	0.89	1.74	0.82	1.62	0.81	2.33
CaCO3	10	4.40	0.93	5.40	0.90	7.13	0.83	2.42
pH	3	0.31	0.38	0.34	0.31	0.42	0.37	1.28
Fed	3	0.49	0.48	0.52	0.43	0.52	0.34	1.26
Fed/Fet	3	2.61	0.37	2.75	0.31	2.56	0.31	1.23

Note: LVs is latent variables, RMSEC is root mean square error for calibration, RMSECV is root mean square error for cross validation, RMSEP is root mean square error for prediction, ${\mathrm{R}}_{\mathrm{c}}^2$ is coefficient of determination for calibration, ${\mathrm{R}}_{\mathrm{cv}}^2$ is coefficient of determination for cross validation, ${\mathrm{R}}_{\mathrm{p}}^2$ is coefficient of determination for prediction, and RPD is residual prediction deviation. Hereinafter inclusive.

. Good and excellent calibrations and prediction (${\mathrm{R}}_{\mathrm{c}}^2$ ≥ 0.80, ${\mathrm{R}}_{\mathrm{p}}^2$ ≥ 0.80 and RPD ≥ 2.30) were obtained for SOM, CaCO₃, and clay. The pH, crystalline Fe and Fe_d/Fe_t produced an acceptable result (${\mathrm{R}}_{\mathrm{c}}^2$ ≥ 0.37, ${\mathrm{R}}_{\mathrm{p}}^2$ ≥ 0.31 and RPD ≥ 1.23).

Complex SOM components complicate the assignment of absorption features to specific functional groups [63,64], although Viscarra Rossel and Behrens [35] concluded that SOC is related to wavelengths that represent absorptions due to organic molecules and proteins. SOM is the most frequently studied soil property via reflectance spectroscopy and the prediction results were superior [20], even at a large scale [65]. In this study, the SOM model provided the best results among six properties (RPD = 3.04).

The CaCO₃-leaching process exhibited an inverse trend with SOM accumulation (Figure 5c and Figure 6c), producing a negative correlation between CaCO₃ and SOM (Table 2). This correlation and the absorption features of carbonate (2345 nm) may render it possible to predict CaCO₃ levels in the soils studied.

The fact that clay is primarily composed of phyllosilicate minerals with special features and the correlation between clay and iron oxides (Table 2) might have contributed to the clay prediction via reflectance spectroscopy [35,63,66]. The pH value, a chemical property, is spectrally featureless, thus its correlation with clay, CaCO₃ and SOM (Table 2) was an important prediction mechanism [67]. The diagnostic absorption features can be used to calibrate the content of Fe_d in soil [68,69,70], which makes it possible to predict Fe_d.

3.5. Soil Age Calibration Using Spectra From Multiple Depths

In this study area, the absolute soil age is indicative of the development degree of the soil profile. Eight partial least squares regressions were constructed with soil age data as the response variable and spectral data as the predictor variable. The predictor variables were multi-depth soil spectral curves after continuum-removal transformation from the topsoil (one layer) to the entire profile (eight layers), respectively.

Soil chronosequence reflectance spectroscopy can be used to predict soil ages (

Table 4. Calibration and predication of the profile age using multi-depth spectroscopy.

Layer Number	Depth	LVs Number	Calibration (n = 10)		Cross Validation		Prediction (n = 6)
			RMSEC	${\mathbf{R}}_{\mathbf{c}}^2$	RMSECV	${\mathbf{R}}_{\mathbf{c}\mathbf{v}}^2$	RMSEP	${\mathbf{R}}_{\mathbf{p}}^2$	RPD
1	0–5 cm	2	131.11	0.80	188.88	0.66	138.72	0.62	1.79
2	0–10 cm	2	137.59	0.77	198.75	0.62	157.58	0.51	1.57
3	0–20 cm	2	141.37	0.76	206.70	0.59	155.19	0.53	1.60
4	0–30 cm	2	133.98	0.79	201.96	0.61	133.63	0.65	1.85
5	0–40 cm	3	99.76	0.88	167.47	0.73	140.78	0.61	1.76
6	0–60 cm	4	80.94	0.92	177.19	0.70	145.15	0.59	1.71
7	0–80 cm	4	78.10	0.93	187.05	0.66	148.85	0.57	1.66
8	0–100 cm	4	72.38	0.94	201.58	0.61	154.29	0.53	1.61

Note: Conventions are shown in Table 3.

). In general, all the models achieved fair or good results (RPD > 1.57). The number of principal components increased with the soil layer number. The calibration results (${\mathrm{R}}_{\mathrm{c}}^2$) decreased slightly and increased remarkably with the soil layer number. Otherwise, the prediction results began to decrease after they reached the best result (RPD = 1.85). That is to say that more principal components extracted from more soil layers did not produce better prediction results. The topsoil regression results (0–5 cm) were acceptable because topsoil properties, such as SOM, clay, and CaCO₃, changed with soil age (Figure 6a–c). Obtaining the chemical and physical properties of more layers provided more information on pedogenic processes than those of the topsoil (0–5 cm). Therefore, better results should be obtained when additional layer spectral curves are used to calibrate the soil age. Although the two- or three-layer spectra did not produce better calibration and prediction results than topsoil as we supposed, the prediction results for four layers (0–30 cm) were best. The possible reason is that four-layer soil contains more pedogenic information than topsoil and the pseudo multi-depth spectra, joined by four-layer spectra, includes enough information about the difference between soil profiles of different degrees of development.

The regression coefficients can show the importance of the specific wavelength. These important wavelengths (Figure 11) included all the absorption features in the continuum-removed spectra in Figure 7. These absorption features are attributed to Fe (372, 444, 535, 666, 908 and 1131 nm), clay minerals (1412, 1914 and 2217 nm) and CaCO₃ (2347 nm). This suggested that these soil properties, producing these absorption features, were very important in the calibration and prediction of profile age.

The correlation coefficients were calculated to analyze which soil properties were correlated with the two latent variables in the topsoil model (

Table 5. Correlation coefficients of latent variables and soil properties.

	SOM	Clay	CaCO3	pH	Fed	Fed/Fet
LV-1 (68%)	0.78 **	0.52 *	−0.62 *	−0.41	−0.14	0.33
LV-2 (12%)	0.46	0.66 **	−0.71 **	−0.37	0.39	0.35

Note: ** significant at the 0.01 level (two-tailed); * significant at the 0.05 level (two-tailed); LV-1: the first latent variable; LV-2: the second latent variable.

). These two latent variables explained 68% and 12% of the spectroscopy information, respectively. The 1st latent variable is positively correlated with SOM and clay and negatively correlated with CaCO₃. The correlation coefficients between the second latent variable and clay and CaCO₃ were significant. Although it is very difficult to assign the absorption features to SOM, SOM is one of the most important properties in the calibration and prediction of profile age. We can conclude that these three soil properties, SOM, clay and CaCO₃, changed with soil age and primarily resulted in the observed changes in the spectra. Thus, the changing soil properties with time act as a bridge and make it possible to analyze soil age or pedogenic processes with spectroscopy.

4. Discussion

4.1. Reflectance Spectroscopy and Soil Genesis

Soil visible and near-infrared (VNIR) diffuse reflectance can provide rich information on soil physical and chemical properties, which implies the possibility of using soil spectroscopy to aid in the analysis of soil genesis. Our results showed that pedogenic processes, such as SOM accumulation and CaCO₃ leaching, can be evaluated, and iron oxides and clay minerals can be identified by soil spectroscopy. The application of soil spectroscopy in soil genesis is mainly divided into the following aspects: firstly, the soil weathering degree or weathering indices can be estimated using spectral reflectance [24,71]. Furthermore, Crouvi et al. [72] found that it is possible to map alluvial ages by measuring spectral parameters of three independent chromophores, namely, iron, clay, and carbonates. Secondly, pedogenetic processes, such as organic carbon accumulation, oxide iron variation, and soil rubification process, were identified with soil spectral reflectance [26,28]. Thirdly, soil genetic horizons were distinguished or recognized, and horizon boundaries were defined using soil spectroscopy [73,74,75]. Fajardo et al. [74] indicated that soil spectroscopy offers more information than the traditional description and avoids observer bias when describing a soil core. In addition, soil classes were identified or allocated with surface or profile spectral reflectance [76,77,78,79,80,81]. Since the use of soil spectral reflectance applied to soil science studies has seen exponential growth in the past decades, a new term “spectral pedology” was defined by Demattê and Terra [28].

4.2. Combination of Horizontal Spectroscopy

The study object of soil genesis should be the entire soil profile, not only topsoil or a certain horizon. Therefore, a critical issue when using this new methodology to characterize the entire soil profile is how to combine the spectroscopy of soil horizons for individual profiles [78]. At present, there are four methods to integrate spectral curves for the soil profile. Viscarra Rossel and Webster [73] averaged the spectra for topsoil and subsoil to represent the spectra of a profile. This simple average method may weaken some spectral features of a certain genetic horizon. Vasques et al. [41] constructed a pseudo multi-depth spectrum for profiles by combining the spectra for three depth layers sequentially. This method effectively retains the spectral information for every horizon, but it can be used only when the number of spectral curves in each profile is the same. Thirdly, Ogen et al. [75] used hyper-depth spectral data to classify soils and distinguish soil horizons, where the X-axis presents the wavelength in nm, the y-axis is the depth in cm, and the color scale displays the reflectance value. This method translated the spectral curves in a profile into a spectral surface and retained all the information contained in soil spectroscopy. However, Fajardo et al. [74] found that the depth interval should be ≤2 cm in order to efficiently capture the soil depth changes in a soil core. Therefore, this spectral surface can be used to classify soil classes, but a smaller sampling interval is needed to distinguish the genetic horizon. The last method is to integrate the spectra of a genetic horizon by depth-weighted averaging [78], similar to the method of Viscarra Rossel and Webster [73]. The above methods serve as efficient tools for integrating spectra of legacy soil samples. For the future profile study, Ben-Dor et al. [82] noted that imaging spectroscopy (IS) technology can provide high spectral and spatial resolution data and has the ability to track chemical and physical properties that can serve as tracers to monitor soil genesis and formation.

4.3. Limitation of Research

This work attempted to estimate the profile age in a chronosequence with combined spectral data and obtained a good performance. In this study, the pedogenic properties, such as SOM, clay, and CaCO₃, which change with time and are spectrally active, act as the bridge between the profile age or pedogenic processes and spectra. However, this result only suggested the feasibility of this study area and more study in this direction is still required for other chronosequences. It is possible that other properties, such as clay minerals and iron oxides, would be important for spectral analysis and age prediction in a chronosequence of tens of thousands or millions of years.

5. Conclusions

Changing soil properties in a chronosequence that are produced through different weathering levels and pedogenic processes can be reflected through spectral behaviors. These properties serve as a bridge linking soil spectroscopy and weathering levels/pedogenic processes. This study, conducted using a coastal soil chronosequence developed on fairly uniform calcareous marine sediments, demonstrated that profile spectroscopy can detect property changes in a chronosequence and be further used to calibrate soil age.

In this millennium chronosequence, soil pedogenic processes, such as SOM accumulation, CaCO₃ leaching and clay migration were detectible. Fe_d and Fe_d/Fe_t did not show remarkable sequential changes with soil age. These three soil properties changing with soil age, SOM, clay and CaCO₃, primarily resulted in the observed changes in the spectra. Therefore, these pedogenic processes can be reflected through the reflectance intensity and absorption features of soil spectroscopy. These soil properties, some of which are featureless, can be calibrated through soil reflectance spectroscopy. The pseudo multi-depth spectral curve of a soil profile, a method used to integrate reflectance spectroscopy data for various depths, provided more information on pedogenic processes and produced better soil age calibration and prediction than a spectral curve of topsoil. However, more profiles are needed in the future to further prove this conclusion, and this idea and method can be used in other regions and calibration models according to local soil samples.

Therefore, soil spectroscopy for chronosequence study not only serves as a new technique for pedogenesis research but also represents a new application of spectroscopy techniques. Further studies on how spectroscopy can be applied to chronosequences of other time scales are needed, e.g., a chronosequence of millions of years. The other properties, such as clay minerals and iron oxides, may be important properties for spectral analysis and age prediction in the thousands or millions of years chronosequence.