Soil organic carbon and its relation with the reflectance spectroscopy in soil sampling by layers and compound in the Depositional Piedmont of Villavicencio (Colombia)

doi:10.21203/rs.3.rs-5219195/v1

Download PDF

Research Article

Soil organic carbon and its relation with the reflectance spectroscopy in soil sampling by layers and compound in the Depositional Piedmont of Villavicencio (Colombia)

https://doi.org/10.21203/rs.3.rs-5219195/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Monitoring and quantifying soil organic carbon content in different landscapes at different sampling depths is crucial to assessing the quality of soil, water and the functionality of the ecosystems that are established there. Monitoring results can guide practices that could be adopted to improve soil organic carbon storage and thus enhance and diversify ecosystem services. Inventories of soil organic carbon stocks at local level can be used to assess compliance with international climate change conventions. In the case of the Depositional Piedmont of Villavicencio (Colombia), the results of the integrated principal component variance analysis in spectral reflectance data with subsequent spatial modelling did not show differences in importance either at different depths of soil sampling or in the landscapes themselves, but showed a clear relationship of soil organic carbon content to bulk density and weighted average reflectance of the 1506.93-1846.38 nm band mainly, followed by the 1000.00-1245.02 nm band. This proposed methodology is simple to adopt and can be used involving a larger set of factors that may be associated with reflectance in experimental or observational data structures.

ANOVA-PCA

Bulk Density

Spatial regression model

Weighted average reflectance

Soil organic carbon (SOC) is the property best related to soil health. It could be attributed to the properties of an adhesive hydrogel because it adheres to soil aggregates, adsorbs and “desorbs” nutrients, retains water which provides habitat for soil microorganisms. Crops do not take carbon from the soil (they take it from the atmosphere) and, as such, carbon (C) is not in itself an essential nutrient; however, due to its capacity to retain water in the soil it is key to plant growth and development and it is fundamental to flood control (Ribeiro et al, 2021). The study of SOC is becoming increasingly important because of its role in the fight against climate change, something that can be related to changes in atmospheric properties, but in general the soil, followed by the oceans, is a key element of the climate system (Paterson and Sim, 2013).

Because of its importance in agricultural production, quantification of SOC over large areas is essential to characterize the soil and effects on its management; however, in large areas a greater sampling effort is usually made to achieve greater spatial coverage, raising the costs of field sampling and subsequent laboratory analysis. These costs rise proportionally to the number of layers considered (Liang et al, 2024), which is why efforts are made to involve sensing technologies that simplify the collection of information related to the distribution of its content. These technologies usually generate a large amount of information, as is the case of reflectance spectroscopy, a technique used as a method for chemical and mineral soil analysis (Nocita et al, 2015). NIR (Near-Infrared) reflectance spectroscopy captures spectral signatures or soil spectra that mainly contain information about the concentration of molecules, organic compounds, and inorganic bonds present in the soil sample; all this information is contained in the soil spectrum as a complex combination of absorption bands or peaks (Alomar and Fuchslocher, 1997). With this technique, multiple responses are generated and incorporated into the data matrix, whether observational or experimental, so the data structure becomes complex by having more columns than rows associated with the repetitions made on the sampling units (big-p, little-n) and in many techniques, where it is required to obtain covariance matrices that characterize variability. Some procedures require their inversion; however, it is not possible to obtain the usual inverse and a pseudoinverse may be used (Hastie et al, 2017).

Of the tools usually applied for dimensionality reduction when dealing with spectroscopy data, one can highlight principal component analysis (Yuliara et al., 2020), t-distributed stochastic neighbor embedding (t-SNE) and its variants (Dabrowska et al., 2020), multidimensional scaling (Mcilwaine et al., 2019), uniform manifold approximation and projection (UMAP) (Smets et al., 2020), linear neighbor embedding (LLE) (Shan et al., 2014), isometric mapping (Isomap) (Tang et al., 2014), autoencoders (Tsimpouris et al., 2021), matrix factorization (Pauca et al., 2006), and regression methods based on PLSR, multiple regression (MLR), and Support Vector Machine Regression (Kamboj et al., 2022), et al. While all these techniques reduce dimensionality, the resulting lower-dimensional dataset often requires additional processing. Many of these methods serve more as tools than final solutions, leading to analyses that typically rely on simple correlation matrices of the generated spectral indices and these matrices represent a simplified view of the spectral response (Ramadhan et al., 2021). Other common analyses include cluster analysis for grouping indices (Marino and Alvino, 2021) and univariate or multivariate analysis of variance with p-value adjustments based on the number of responses (reflectance) or latent variables (components) (Riefolo et al., 2021; Tao et al., 2022).

The objective of this research was to relate soil organic carbon (SOC) content to reflectance spectroscopy measurements taken at three soil depths from samples collected in three different landscapes (the first classification factor) in the Depositional Piedmont of Villavicencio (Colombia). Additionally, the study considered samples from three consecutive soil layers as well as a composite sample of these layers. This structure creates two classification factors with 12 scenarios (3 landscapes and 4 samples based on depth), enabling the application of analysis of variance alongside principal components analysis (aov_pca). This combination provides a visual representation to evaluate the interaction effects of the classification and principal factors (Harrington et al., 2005). The results indicate that the first two components explain 98.5% of the variance, showing no interaction between the layered and composite samples. Additionally, there is greater similarity between the mountain and piedmont landscapes compared to the plain landscape in terms of reflectance. This pattern aligns with the SOC content across regions, whether considering layered or composite samples.

Study area

The study was conducted in the municipality of Villavicencio, located in the Depositional Piedmont of the Eastern Cordillera in the department of Meta, southeastern Colombia (Latitude (4° 3’2.67”, 4°14’5.31”, 4° 7’47.95”, 4° 2’35.13” N) Longitude (73°45’40.61, 73°37’29.84, 73°14’20.59, 73°15.39’44.76”)). The total study area covers 673 Km², with altitudes between 232 and 1349 m, representing 51% of the total area of the municipality. According to the Köppen classification system, the regional climate is Am - tropical monsoon; and according to the Holdridge classification system the life zone of the study area is BmhT - very humid tropical forest. This formation corresponds to the transition zone between the mountainous landscape of the Eastern Cordillera and the Depositional Piedmont that occupies the central region of the department (IGAC, 2004).

Geology, geomorphology and soils

With the differential uplift of the Eastern Cordillera during the Quaternary, the Paleozoic formations formed by the alternation of schists with argillites, sandstones, gneiss, and phyllites currently outcrop in most of the Eastern Cordillera and give rise to the mountain soils. At the end of the Pleistocene, the erosive cycle of the Eastern Cordillera and sedimentation in the plain were reactivated due to the last glaciation, since the Cretaceous material of San Juanito, El Calvario and Páramo de Sumapaz (sandstones, argillites and shales) was finally eliminated almost in its entirety and deposited in the partially subsided foothills. These materials of clays and sandy conglomerates were deposited in the foothills, following the course of the main rivers such as the Guayuriba, Guatiquía, et al., as part of their alluvial plains. At the end of the glaciations in the Holocene, the level of regional erosion decreased and created the most recent geoforms that are currently identified in the plains landscape. The rivers also decreased leaving their level below the earlier river level, thus forming a system of terraces. The sedimentation pattern at the foot of the Eastern Cordillera is typically fluvial and comprises a vast plain extending from the Cordillera to the Orinoco River (IGAC, 2004).

The region is dominated by inceptisols, oxisols and, to a much lesser extent, entisols, all of which have developed mainly on a depositional piedmont landscape; about 54% of the total study area corresponds to this type of landscape followed by 26% in an alluvial plain and 13% in a mountain landscape.

Climate

The study area encompasses two major bioclimatic units: the Eastern Cordillera area, a very humid medium climate (m-MH) with precipitation above 4000 mm per year and temperatures that decrease as altitude increases; and the Depositional Piedmont with a warm humid climate (c-MH) where precipitation decreases and temperatures gradually increase as one approaches the alluvial plain from west to east (IGAC, 2004). In general, the Depositional Piedmont climate has an average temperature of 27°C and a bimodal precipitation regime that is largely determined by the Intertropical Convergence Zone and the influence of the mountainous relief of the Eastern Cordillera. The first rainy season occurs between April and June, and a second rainy season peaks in October. The average annual rainfall is 4260 mm and the dry season includes the months of December to February with less than 200 mm in total (Varón and Vargas, 2019).

Soil sampling and laboratory analysis

A regular sampling grid with 12.5m×12.5 m for every 2 km was established at the study site. This resulted in 152 samples collected at different locations. At each location one sample was compose for three different soil sample for each soil layer 0–10 cm, 10–20 cm and 20–30 cm and a one composite sample for the depth 0–30 cm. The soil samples were air-dried, ground and sieved (2 mm mesh). The total carbon content was determined using dry combustion (Dumas method) in a CN Elemental Analyzer (LECO Inc., St. Joseph, MI, USA). Analysis of soil particle size distribution was performed using the hydrometer method, in which sodium hydroxide (4 g L^-1) and sodium hexametaphosphate (10 g L^-1) were employed as dispersing agents. Elevation data are from the ALOS/PALSAR digital elevation model (DEM) developed by Alaska Satellite Facility (ASF) of spatial resolution 12.5 m × 12.5 m (https://nasadaacs.eos.nasa.gov/).

Acquisition and pre-processing of NIR spectra

To obtain spectral measurements, the samples were oven dried for 24 hours at 45°C and they were ground and sieved (2 mm mesh). Subsequently, a modular FT-NIR NIRFlex N-500 spectrometer (https://www.buchi.com/se-en/node/822) was used to measure the NIR reflectance spectrum for each soil sample, 152 in total. The data were obtained with a resolution of 0.2 nm between 1000 and 2500 nm. Before scanning the soil spectra, the soil samples were placed in petri dishes (9 cm in diameter by 1.5 cm in height) and distributed homogeneously to reduce the influence of the roughness of the material. The instrument was calibrated using a standard Spectralon white plate, and this procedure was performed at first and repeated every 20 minutes during the readings. For each soil sample, the spectrometer automatically calculated an average of 100 scans. In addition, three independent measurements were taken for each soil sample, mixing the soil sample in each measurement to ensure a better representation of the examined surface. The final spectrum was finally calculated as an average of the three spectral scans. NIR reflectance spectra were converted to apparent absorbance (A) using the expression A = log(1/reflectance); the logarithm being natural in all cases. Multiplicative interference and light scattering effects were reduced by normalizing spectra using the standard variable normal algorithm (SNV) (Barnes and Lister, 1989).

Statistical analysis

Once the reflectance data matrix was available for the soil samples of the three landscapes being considered (first factor: mountain, piedmont and plain) and the three soil depths and a composite sample of the three depths (second factor: 0–10 cm, 10–20 cm, 20–30 cm and 0–30 cm), a matrix of dimension 72⋅(2 + 1501) was generated, using two columns for factors and 1501 reflectances, where two columns correspond to the structure of the classification factors of the observational study; this generated 12 combinations of levels. In order to carry out this analysis, two options are proposed: the first one, the evaluation of the whole spectrum using the method combining the AOV variance analysis and the principal component analysis, not for a run as it would be done separately with each technique but generating a series of sub-matrices associated to the two-factor model with interaction given by the decomposition of the response matrix, where it represents the landscapes (l = 1,2,3), d represents the different depths (d = 1,2,3,4) and the number of samples taken (s = 1,2,...,6). The matrix y represents the 1501 soil reflectances for the combined l x d x s = 72 rows of data. The decomposition for a single vector of particular reflectances is written as:

$$\:{y}_{lds}=\stackrel{-}{y}+\left({\stackrel{-}{y}}_{l.}-\stackrel{-}{y}\right)+\left({\stackrel{-}{y}}_{d.}-\stackrel{-}{y}\right)+\left({{\stackrel{-}{y}}_{ld}-\stackrel{-}{y}}_{l.}-{\stackrel{-}{y}}_{d.}+\stackrel{-}{y}\right)+\left({y}_{lds}-{\stackrel{-}{y}}_{ld}\right)$$

where$\:\:\stackrel{-}{\varvec{y}}$ represents the vector of overall reflectance averages, $\:{\stackrel{-}{y}}_{l.}$ represents the vector of means for the first factor with its levels, $\:{\stackrel{-}{y}}_{d.}$ represents the vector of means of the second factor and its levels, and $\:{\stackrel{-}{y}}_{ld}$ represents the vector of means for the combination of levels of both factors. The first operation performed by the algorithm is to standardize the responses, which is achieved with $\:{\varvec{D}}_{1}={\varvec{y}}_{lds}-\stackrel{-}{\varvec{y}}$ (Harrington et al, 2005 notation), thus using the whole spectrum the first working submatrix of the algorithm is achieved. For the first factor analysis associated with the landscape, the effect matrix$\:\:\stackrel{-}{{\varvec{D}}_{\varvec{A}}}=\left({\stackrel{-}{\varvec{y}}}_{l.}-\stackrel{-}{\varvec{y}}\right)$ is calculated for all the reflectances considered in the model, so that using Eq. 1 for a single reflectance would mean that$\:\:{D}_{2}={y}_{lds}-\stackrel{-}{y}-\left({\stackrel{-}{y}}_{l.}-\stackrel{-}{y}\right)={D}_{1}-\stackrel{-}{{D}_{A}}$ (remember that the bold letters refer to all the reflectances, so that the normal text only illustrates a particular reflectance). Similarly, for the second factor associated with depth, we obtain (in a sequential work, i.e., the second factor$\:\:\stackrel{-}{{\varvec{D}}_{\varvec{B}}}=\left({\stackrel{-}{\varvec{y}}}_{d.}-\stackrel{-}{\varvec{y}}\right)$ after considering the first one) first the effects matrix of the second factor and then we obtain $\:{D}_{3}={y}_{lds}-\stackrel{-}{y}-\left({\stackrel{-}{y}}_{l.}-\stackrel{-}{y}\right)-\left({\stackrel{-}{y}}_{d.}-\stackrel{-}{y}\right)={D}_{1}-\stackrel{-}{{D}_{A}}-\stackrel{-}{{D}_{B}}={D}_{2}-\stackrel{-}{{D}_{B}}$. With the calculation of D₃ for a reflectance (or D₃ for the whole spectrum or some subset selected by related theory or by some other applied method). Eq. 1 operated sequentially has on the right-hand side of the equality only the interaction effect and that associated with the error, which in the case of a model with no interaction and a single reflectance observation for each level of the factors (s=1) involved would be written as:

$$\:{y}_{ld}=\stackrel{-}{y}+\left({\stackrel{-}{y}}_{l.}-\stackrel{-}{y}\right)+\left({\stackrel{-}{y}}_{d.}-\stackrel{-}{y}\right)+\left({y}_{ld}-{\stackrel{-}{y}}_{l.}-{\stackrel{-}{y}}_{d.}+\stackrel{-}{y}\right)$$

Thus, if we have the case of Eq. 2, $\:{\varvec{D}}_{3}$, it would already represent the error matrix, however, in the present case, it is still necessary to calculate$\:\stackrel{-}{{\:\varvec{D}}_{\varvec{A}\varvec{B}}}$ the one associated with the interaction effect, so that the error matrix would be obtained from$\:\:{\varvec{D}}_{4}={\varvec{D}}_{E}={\varvec{y}}_{lds}-{\stackrel{-}{\varvec{y}}}_{ld}$.

The evaluation of the first factor with the technique that integrates ANOVA and PCA requires the matrix$\:{\:\varvec{M}}_{A}=\stackrel{-}{{\varvec{D}}_{\varvec{A}}}+{\varvec{D}}_{E}$ to be obtained and it is here where PCA enters on this matrix using the decomposition by singular values of$\:{\:\varvec{M}}_{A}={\varvec{U}}_{A}{\varvec{S}}_{A}{\varvec{V}}_{A}^{T}$, where$\:\:{\varvec{V}}_{\varvec{A}}$ compresses the loadings of the variables, $\:{\varvec{S}}_{A}$ is a diagonal matrix with the singular values and $\:{\varvec{U}}_{A}$ compresses the normalized scores. With the matrix product of $\:{\varvec{U}}_{A}{\varvec{S}}_{A}$, we obtain the values that are finally plotted and that in the technique could reflect the differences attributable to the effect of the first factor. This same procedure is followed for each effect involved in the study to generate the characteristic outputs of the technique, which are used to study the interaction effect, as well as the main effects.

Once the interaction of factors and the main effects of the first database had been studied, we proceeded to work with a second database generated at a single depth (thanks to the results of the first analysis) for all the landscapes. First, a plot was created between the ratio of soil bulk density (BD) and soil organic carbon with each of the reflectances measured at each wavelength. For each iteration of three variables (fixing bulk density and SOC), we obtained the correlation matrix, extracted the first eigenvalue and percentualized with respect to the sum of all other eigenvalues. With the data of the first eigenvalue it was generated a plot for all the lengths, which was used to apply six classification algorithms and thus select three bands from the whole spectrum evaluated intersecting all the ranges obtained from each cluster created; in our case this was three groups (suggested by the algorithm). We proposed a weighted average reflectance ($\:\stackrel{-}{{\rho\:}_{w}}$) for any band with the expression:

$$\:\stackrel{-}{{\rho\:}_{w}}=\frac{\sum\:_{{j=\lambda\:}_{\text{m}}}^{{j=\lambda\:}_{\text{M}}}\sum\:_{i=1}^{n}{\rho\:}_{ij}{\lambda\:}_{j}}{\sum\:_{j={\lambda\:}_{\text{m}}}^{{j=\lambda\:}_{\text{M}}}{\lambda\:}_{j}}$$

where $\:{\lambda\:}_{j}$ represents the j-th wavelength (j = 1,...,_mk), where m_k is an indicator of the number of wavelengths within each band (which can have different sizes), represents the reflectance observed at the i-th observation and j-th wavelength, and $\:{\lambda\:}_{\text{m}}\:$and $\:{\lambda\:}_{\text{M}}$ represent respectively the minimum and maximum wavelength value of a band. With these calculations, we obtained the three weighted mean reflectance vectors, which allowed forming the final data matrix to estimate soil organic carbon using a spatial regression model, specifically the spatial model of the error, which according to Elhorst's (2014) taxonomy is written as:

$$\:\:\:\:\:\:\varvec{Y}=\alpha\:{\varvec{\iota\:}}_{n}+\varvec{X}\varvec{\beta\:}+\varvec{u}\:;\varvec{u}=\lambda\:\varvec{W}\varvec{u}+\varvec{\epsilon\:}$$

where $\:\varvec{Y}$ represents the vector of observed SOC value, $\:\varvec{\alpha\:}{\varvec{\iota\:}}_{\varvec{n}}$ represents the intercept (vectorially), $\:\varvec{X}$ the matrix of explanatory variables, namely the weighted mean reflectances of each of the selected bands and the bulk density, $\:\varvec{\beta\:}$ represents the vector of parameters in the explanatory variables, $\:\varvec{u}$ represents the residuals of the model (with spatial dependence), which is explained by the expression $\:\varvec{u}=\lambda\:\varvec{W}\varvec{u}+\varvec{\epsilon\:}$, where $\:\lambda\:$ represents the spatial autocorrelation coefficient, $\:\varvec{W}$ the matrix of spatial weights, created with the inverses (element to element) in the matrix of Euclidean distances (d) following the pattern d^-1/2 and with a structure of neighborhoods given by the nearest neighbor (with W-type standardization). Finally, $\:\varvec{\epsilon\:}$ represents the vector of normal, independent and identically distributed residuals.

Using the modeling results and metrics for model fit (statistical significance of the spatial autocorrelation coefficient), along with the fulfillment of assumptions (normality and independence of the type $\:\varvec{\epsilon\:}$ residuals, evaluated with the Shapiro-Wilk test and the Moran index), we estimated soil organic carbon values by aggregating data in hexagons with a 3 km diameter in QGIS (Qgis.org, 2024), adjusting the anchor point from the software's default setting. This generated a table of attributes that we append in the supplementary material in the GitHub address (Supplementary Materials). As we built the spatial model based on regular polygons, some landscapes fell into the same polygon for some of the 72 polygons. The estimated value of each of these was assigned equally to each landscape according to its location in the aggregation. With these values, we constructed a scatter plot between the SOC estimated by spatial regression and the mean bulk density of each polygon per landscape. We used the ordinary least squares method to adjust the (non-spatial) regression line to study the bulk density domain and the values where the SOC oscillates in each landscape to characterize the bulk density and the estimated SOC per landscape. We calculated all PCA analyses using the R::ChemoSpec library (Hanson,2024), spatial regression with R::spatialreg (Bivand and Piras, 2015) and Moran index calculations with R::ape (Paradis and Schliep,2019).

Initially, we described the distribution of spectral reflectance in the combination of landscapes and sampling depths (Fig. 1).

Next, we described the results obtained from the application of aov_PCA in the whole evaluated spectrum, ending with the interaction of factors in Fig. 2C, and showing in the upper panels the factors associated to the landscape (Fig. 2A) and the factor associated with the sampling depth (Fig. 2B).

As the previous results showed no differences in reflectance values for the depths, we used the results for a new measurement at points close to the coordinates of the first database and recorded as a second database where SOC, bulk density and spectral curves were obtained again. (Both databases are shown in supplementary material).

In the second database there are 1501 wavelengths with 140 reflectances evaluated for the same three previous landscapes, as well as the SOC and bulk density data. With this information we proceeded to build an algorithm to generate the first eigenvalue of the SOC-BD-reflectance triplet (keeping SOC and BD variables as fixed). We generated this eigenvalue from each correlation matrix (Pearson correlation matrix 3⋅3 dimension) and percentualized with respect to the sum of all eigenvalues. The associated plot is illustrated in Fig. 3.

With the vector of the first eigenvalue, we applied different classification algorithms to group similar bands in terms of observed reflectance. We proposed three groups as the diagrams suggested in all cases. We generated the common band with the intersection of all the ranks generated by each clustering method from the bands provided by each algorithm (Table 1).

Table 1

Selection of a common band by different clustering algorithms.
Method	Cluster 1	Cluster 2	Cluster 3
K means	1000.00- 1401.35	1402.13-1898.25	1899.7–2500.00
PAM	1000.00-1357.22	1357.96-1842.30	1843.66–2500.00
Clara	1000.00-1245.02	1245.64-1896.81	1898.25–2500.00
Hierarchical	1000.00-1506.02	1506.93-1846.38	1847.75–2500.00
Fuzzy	1000.00-1329.08	1329.79-1823.49	1824.82–2500.00
Model based	1000.00-1390.43	1391.21-1875.47	1876.88–2500.00
Banda común	1000.00-1245.02	1506.93-1846.38	1899.7–2500.00

The common band in table two allowed the generation of the weighted mean reflectance vector for each band, which with the bulk density and SOC data allowed the fitting of the spatial model given in (2), from which the results in Table 2 were generated. The image of the aggregated data with the estimated SOC values per landscape for a hexagonal grid of 3 km radius is illustrated in Fig. 4.

Table 2

Spatial regression results based on the spatial error model.
Coefficients	Estimate	Std. Error	Z value	Pr(>\|z\|)
Intercept	9.693	1.792	5.409	6.311e-8
Bulk density	-5.945	1.330	-4.469	7.874e-6
$\:\stackrel{-}{{\varvec{\rho\:}}_{\varvec{w}}}$: 1000.00-1245.02 (nm)	-27.085	6.524	-4.152	3.303e-5
$\:\stackrel{-}{{\varvec{\rho\:}}_{\varvec{w}}}:$ 1506.93-1846.38 (nm)	21.604	5.171	4.178	2.940e-5

Table 3

The error spatial model fits statistics with outputs associated with a revision of assumptions.
	Lambda	LR test	Wald test	Normality	Moran Index
Statistic	0.267	3.723	4.314	W = 0.987	-0.012
P value		0.0536	0.0378	0.693	0.633

Table 2 showed the rejection of the hypotheses associated with the null coefficient of the parameters associated with the bulk density and the 1000.00-1245.02 (nm) and 1506.93-1846.38 (nm) bands. This used the weighted mean reflectance as the vector representing the interval. This fitted model generated the estimated spatial autocorrelation coefficient of 0.267, which, using the likelihood ratio tests (LR test) and the Wald test, allowed the model to be selected to explain the estimated carbon per cell. The latter measures are only useful for the revision of assumptions, especially the p-value of the Moran index where the spatial lag in the residuals is expected to explain the spatial dependence structure found and is finally lost with the incorporation of the term $\:\varvec{u}=\lambda\:\varvec{W}\varvec{u}+\varvec{\epsilon\:}$ of Eq. (2) (Table 3).

Finally, we conducted an optimization to visualize the relationship between BD and SOC estimated by the spatial model, allowing us to validate the typical connection between these two variables using a first-order linear model by landscape. In Fig. 5, “y” in the model should be understood as the expected value of the estimated SOC given the bulk density (x). The image shows the usual relationship between SOC and bulk density (the higher the bulk density the lower the SOC content); however, the domain of the bulk density is different in each landscape, as is the range associated with the model.

Modeling spectral data remains a great challenge due to the high dimensionality involved, since reflectance spectrometry-based technologies provide vectors of values over a wide range of the electromagnetic spectrum. When this data structure is added to the data structure of an experimental or observational study, complication increases, because now when factors are involved, it is understood that each vector of reflectances at each length becomes the matrix of responses to which the factors included are associated through some modeling process. However, the high dimensionality of the responses limits the use of various statistical techniques, so it is usual to find the creation of spectral indices using some wavelengths that can indicate important changes that can be associated with the presence of chemical compounds of interest (Falcioni et al, 2023). Other researchers prefer the selection of subsets of bands of the spectrum using different algorithms, such as genetic algorithms (Weichao et al, 2022).

We chose this second option, using an algorithm that incorporated analysis of variance to account for the study's factorial structure and principal component analysis to reduce dimensionality. This approach allowed us to visually represent the interactions of factors related to landscapes (mountain, piedmont, and plain) from which soil samples were taken and evaluated at different depths (0–10 cm, 10–20 cm, 20–30 cm, and a composite of 0–30 cm) for spectral reflectance. Pinto et al. (2008) used this on MIR (mid-infrared) spectroscopy data to evaluate how concentration and temperature affected the extraction of gels from certain algae; however, there is little research in the literature where this methodology has been used.

Applying this method to study reflectance between 1000 nm and 2500 nm, using a database of 1785⋅(2 + 1501), revealed no interactions among factors (as the concentration ellipses nearly overlap). This indicated a clear lack of differences in reflectance across various depths and landscapes, with mountain and piedmont landscapes being more similar. There was no differentiation in the reflectance patterns between depths, which contrasts with findings from several authors who reported lower spectral reflectance at shallower depths (0-7.5 cm), followed by 7.5–15 cm and 15–30 cm. This trend is typically linked to changes in soil color, as organic matter decreases with increasing depth (Das et al., 2023; Kogge et al., 2021). In general, reflectance gradually decreased with increasing SOC due to absorption by SOC from a descriptive point of view; however, the aov_PCA did not show differences in reflectances of significance enough to state that the depths resulted in contrasting depths.

For the case of individual and interaction factors, the plot was represented with an explanation of the total variability of 98.1% in the space of the first two principal components. Although the diagram associated with the landscapes and interaction is somewhat similar, it is important to note that only a small fraction of the observations of the plains landscape causes the associated concentration ellipse to show a smaller overlap, which could suggest some measure of interaction, which in our case we consider negligible due to the large number of points that overlap in the various ellipses. It is not strange to find in the Colombian case results showed no differences in different landscapes despite the influence that practices associated to soils with different crops have had (Bolívar et al, 2019). However, it is best not to generalize since the physical-chemical and microbiological properties of soils are as diverse as the climatic conditions that make each study a conditional truth inherent to the region where the description and comparison is made.

Based on our first result, we built a second database to evaluate reflectances again, but this time at a single depth, using previous findings. We focused on soil properties like SOC and bulk density, as many studies highlight the link between these properties and different landscapes (Chen et al., 2024; Gebre, 2017). The associated costs in the SOC evaluation only allowed the generation of a database of dimension 140⋅(1501 + 2 + 2) (wavelengths, latitude and longitude coordinates, the two soil properties). This allowed the proposal of an iterative method by generating an eigenvalue on the pair of soil variables and each of the wavelengths.

From the 1501 Pearson correlations, a signature was obtained from the first eigenvalue (normalized by the sum of all the priors in each matrix), to which different classification algorithms were applied yielding the bands 1000.00-1245.02 nm, 1506.93–1846.38 nm and 1899. 7-2500.00 nm, which represents approximately a 25% reduction of the bands Exploratory work was done with the first two eigenvalues and the band that best explains the total variability falls approximately in the 1506.93-1846.38 nm (with approximately 86%). With this first dimensionality reduction, a band-weighted reflectance value was proposed to represent the band reflectance, which generated three reflectance vectors (one per band). This proposal forms a matrix that only involves the SOC, the bulk density and the three weighted reflectances. The georeferencing information was incorporated to this structure and regular hexagonal spatial polygons of 3 km in diameter were constructed, which yielded 140 polygons that were adjusted at their anchor point to minimize the effect of polygons without information in the neighborhood pattern. With this aggregation, the spatial regression model was constructed from which the third weighted mean reflectance was discarded (due to statistical significance) and both the bulk density and the first two reflectance vectors were kept as variables related to the average soil organic carbon per cell. The weighted mean reflectance of the first and second band and bulk density were found to be important variables in explaining the aggregate SOC content. Some studies agree with what was found in the second band of 1506.93-1846.38 nm, which largely overlaps with the 1400–1800 nm band found by several authors (Faria et al, 2023, Moura and Bueno, 2018). This band may be related to the water molecule in clay-textured soils and linked to cellulose in organic residues (Jiang et al, 2017).

From the spatial model, SOC was estimated and assigned to the hexagonal grid cells, and as in this aggregation some cases presented sampling points from two or all three landscapes. The individual measurement was used to show, through a least squares adjustment by landscape, the clear relationship between bulk density and SOC in the three landscapes. The slopes were quite similar, although bulk density had a smaller range in the mountain landscape (1.20–1.30 g/cm³ for bulk density and 1.8–2.3 for SOC), while the piedmont and plain landscapes showed more similarity (1.22–1.40 g/cm³ for bulk density and 1.0-2.4 for SOC). These results align with the PCA analysis, showing no interaction.

The evaluation using aov_PCA as an integrated form of analysis of variance and principal component analysis of soil reflectance between 1000–2150 nm electromagnetic spectrum for the study of SOC content and its relationship with different sampling depths (0–10 cm,10–20 cm, 20–30 cm and 0–30 cm) as well as in different landscapes (mountain, piedmont, and plain) of the Depositional Piedmont of the municipality of Villavicencio (Colombia) showed that in the absence of interaction of factors, there were no differences in the reflectance for each of the depths or for the landscapes, although the descriptive measures showed certain variations, they were not important according to the space of the first two components, which explained more than 98% of the total variability in the data.

The method consisted of extracting the first eigenvalue of each correlation matrix between bulk density, SOC content and reflectance. Using six classification algorithms, common bands were identified and summarized into a single weighted mean reflectance vector, grouped into hexagonal polygons of 3 km diameter. This process estimated SOC values for a single depth (0–30 cm) per landscape, which were statistically related to the first two reflectance vectors and bulk density. This validated the close relationship between bulk density and SOC content across all landscapes, with the best relationship found in the 1506.93-1846.38 nm band.

Supplementary Materials

https://github.com/JorgeJola/Modeling_SoilOrganicCarbon.git

Funding

Private funds were not necessary; it was covered with the participants' employment income.

Alomar D, Fuchslocher R (1997) Fundamentos de la espectroscopia de reflectancia en el infrarojo cercano (NIRS) como método de análisis de forrajes. Agro Sur, 26, 88–104. https://doi.org/10.4206/agrosur.1998.v26n1-11
Barnes RJ, Dhanoa MS, Lister SJ (1989) Standard normal variate transformation and de-trending of near-infrared diffuse reflectance spectra. Applied Spectroscopy, 43(5), 772–777. https://doi.org/10.1366/0003702894202201
Bivand R, Piras G (2015) Comparing Implementations of Estimation Methods for Spatial Econometrics. Journal of Statistical Software,63 (18), 1-36. https://doi.org/10.18637/jss.v063.i18
Bolívar A, Camacho C, Ordoñez N, Gutiérrez J, Álvarez G, Guevara M, Olivera C, Olmedo GF, Bunning S, Vargas R (2021) Estimación de carbono orgánico del suelo en Colombia, una herramienta de gestión del territorio. Ecosistemas 30(1): 2019. https://doi.org/10.7818/ECOS.2019
Chen S, Chen Z, Zhang X, Luo Z, Schillaci C, Arrouays D, Richer-de-Forges AC, Shi Z (2024) European topsoil bulk density and organic carbon stock database (0–20 cm) using machine-learning-based pedotransfer functions. Earth System Science Data. Volume 16, issue 5. ESSD, 16, 2367–2383. https://doi.org/10.5194/essd-16-2367-2024.
Dabrowska A, Kania A,Kopeć D (2020) The t-SNE Algorithm as a Tool to Improve the Quality of Reference Data Used in Accurate Mapping of Heterogeneous Non-Forest Vegetation. Remote Sens. 12(1), 39; https://doi.org/10.3390/rs12010039.
Das B, Chakraborty D, Kumar V, Das D, Narayan R, Aggarwal P, Murgaokar D, Prasad B (2023)Partial least square regression based machine learning models for soil organic carbon prediction using visible–near infrared spectroscopy. Geoderma Regional, Vol. 33,e00628. https://doi.org/10.2139/ssrn.4178232
Elhorst PJ (2014) Spatial econometrics from cross-sectional data to spatial panels. Springer. https://doi.org/10.1007/978-3-642-40340-8
Faria O, Torres G, Di Raimo L, Couto E (2023) Estimate of carbon stock in the soil via diffuse reflectance spectroscopy (vis/nir) air and orbital remote sensing. Rev. Caatinga, Mossoró, Vol. 36, n. 3, p. 675 – 689. http://dx.doi.org/10.1590/1983-21252023v36n320rc
Falcioni R, Antunes WC, Demattê JAM, Nanni MR. Reflectance Spectroscopy for the Classification and Prediction of Pigments in Agronomic Crops. Plants (Basel). (2023);12(12):2347. https://doi.org/10.3390/plants12122347.
Gebre AB (2017) Comparison of bulk density methods in determining soil organic carbon storage under different land use types. Journal of Soil Science and Environmental Management. Vol. 9(1), pp. 13-20. https://doi.org/10.5897/JSSEM2017.0664
Hanson B (2024) _ChemoSpec: Exploratory Chemometrics for Spectroscopy_. R package version 6.1.10, <https://CRAN.R-project.org/package=ChemoSpec>. https://doi.org/10.32614/CRAN.package.ChemoSpec
Harrington P, Vieira N, Espinoza J, Nien J, Romero R, Yergey A (2005) Analysis of variance–principal component analysis: A soft tool for proteomic discovery. Analytica Chimica Acta 544,118–127. https://doi.org/10.1016/j.aca.2005.02.042
Hastie T, Tibshirani R, Friedman J (2017) The Elements of Statistical Learning Data Mining, Inference and Prediction. Second Edition. Springer Series in Statistics. USA.
IGAC. (2004). Estudio General de Suelos y Zonificación de Tierras del Departamento del Meta. Escala 1:100.000. Bogota, Colombia: Instituto Geográfico Agustín Codazzi Aberto, 6(1), 35–44. Retrieved from https://revistas.ufrj.br/index.php/EspacoAberto/article/viewFile/5244/3852
Jiang Q et al. (2017) Estimation of soil organic carbon and total nitrogen in different soil layers using VNIR spectroscopy: Effects of spiking on model applicability. Geoderma, 293: 54 -63. https://doi.org/10.1590/1983-21252023v36n320rc
Kamboj U, Guha P, Mishra S (2022) Comparison of PLSR, MLR, SVM regression methods for determination of crude protein and carbohydrate content in stored wheat using near Infrared spectroscopy. Materials Today: Proceedings,48, 576-582. https://doi.org.10.1016/j.matpr.2021.04.540
Kogge G, Kogge R, Palmer B (2021) Soil Organic Carbon Distribution in a Humid Tropical Plain of Cameroon: Interrelationships with Soil Properties. Applied and Environmental Soil Science. Vol. 2021, ID 6052513, 18 pages. https://doi.org/10.1155/2021/6052513
Liang Y, Rillig MC, Chen H, Shan R, Ma Z (2024) Soil pH drives the relationship between the vertical distribution of soil microbial biomass and soil organic carbon across terrestrial ecosystems: A global synthesis. Catena 238,107873. https://doi.org.10.1016/j.catena.2024.107873
Marino S, Alvino A (2021) Vegetation Indices Data Clustering for Dynamic Monitoring and Classification of Wheat Yield Crop Traits. Remote Sens. 13(4), 541. https://doi.org/10.3390/rs13040541
Mcilwaine B, Rivas M, Leinster P (2019) Using 1st Derivative Reflectance Signatures within a Remote Sensing Framework to Identify Macroalgae in Marine Environments. Remote Sens. 11, 704; https://doi.org.10.3390/rs11060704
Moura-Bueno J M et al. (2018) Erosão em áreas de encosta com solos frágeis e sua relação com a cobertura do solo. Revista Scientia Agrária, 19: 102-112. http://dx.doi.org/10.5380/rsa.v19i1.53738
Nocita M, Stevens A, van Wesemael B, Aitkenhead M, Bachmann M, Barthès B, Wetterlind J (2015) Soil Spectroscopy: An Alternative to Wet Chemistry for Soil Monitoring. Advances in Agronomy, 132(March), 139–159. https://doi.org/10.1016/bs.agron.2015.02.002
Paradis E, Schliep K (2019) ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics,35,526-528. https://doi.org.10.1093/bioinformatics/bty633
Paterson E, Sim A (2013) Soil-specific response functions of organic matter mineralization to the availability of labile carbon. Global Change Biology, 19(5), 1562–1571. https://doi.org/10.1111/gcb.12140
Pauca VP, Piper J, Plemmons RJ (2006) Nonnegative matrix factorization for spectral data analysis. Linear Algebra and its Applications, 416, 29–47. https://doi.org/10.1016/j.laa.2005.06.025
QGIS.org, (2024) QGIS Geographic Information System. QGIS Association. http://www.qgis.org
Ramadhan FD, Hariyanto T, Handayani HH (2021) Analysis of Urban Built-Up Expansion Based on Combination of Spectral Indices in Surabaya City. Geoid Vol. 17, No. 1,21-37. http://dx.doi.org/10.12962/j24423998.v17i1.8676
Ribeiro SG, dos Santos A, Rabelo M, Gomes M, da Silva I, Jario L, Bezerra F (2021) Soil Organic Carbon Content Prediction Using Soil-Reflected Spectra: A Comparison of Two Regression Methods. Remote Sens., 13(23), 4752; https://doi.org/10.3390/rs13234752
Riefolo C, Antelmi I, Castrignanò A, Ruggieri S, Galeone C, Belmonte A, Muolo MR, Ranieri NA, Labarile R, Gadaleta G, Nigro F (2021) Assessment of the Hyperspectral Data Analysis as a Tool to Diagnose Xylella fastidiosa in the Asymptomatic Leaves of Olive Plants. Plants (Basel). 10(4):683. https://doi.org/10.3390/plants10040683
Shan R, Wensheng C, Shao X (2014) Variable selection based on locally linear embedding mapping for near-infrared. Chemometrics and Intelligent Laboratory Systems. Vol. 131, pp. 31-36. https://doi.org/ 10.1016/S1872-2040(17)61065-X
Smets T, Waelskens E, De Moor B (2020) Uniform Manifold Approximation and Projection. : Anal. Chem. 92, 5240−5248. https://doi.org/ 10.1021/acs.analchem.8b05827
Tang X, Gao K, Cheng H, Ni G (2014) Hyperspectral unmixing based on ISOMAP and spatial information. Optik, 125,4283-4287. https://doi.org/ 10.14358/PERS.87.7.445
Tao M, He Y, Bai X, Chen X, Wei Y, Peng C, Feng X (2022) Combination of spectral index and transfer learning strategy for glyphosate-resistant cultivar identification. Front. Plant Sci., Sec. Technical Advances in Plant Science. Vol. 13, https://doi.org/10.3389/fpls.2022.973745
Tsimpouris E, Tsakiridis NL, Theocharis JB (2021) Using autoencoders to compress soil VNIR–SWIR spectra for more robust prediction of soil properties. Geoderma, Vol. 393, 114967. https://doi.org/ 10.1016/j.geoderma.2021.114967
Varón SD, Vargas G (2019) Análisis de la susceptibilidad por inundaciones asociadas a la dinámica fluvial del río Guatiquía en la ciudad de Villavicencio, Colombia. Cuadernos de Geografía: Revista Colombiana de Geografía, 28(1), 152–174. https://doi.org/10.15446/rcdg.v28n1.70856
Weichao Sun, Shuo Liu, Xia Zhang, Yao Li (2022) Estimation of soil organic matter content using selected spectral subset of hyperspectral data. Geoderma, Vol. 409, https://doi.org/10.1016/j.geoderma.2021.115653.
Yuliara IM, Ratini NN, Windarjoto, Tri Suandayanid NK (2020) Spectral Reflectance and Principal Component Analysis on the Distribution of Clove Vegetation Using Landsat 8. International Journal of Physical Sciences and Engineering. Vol. 4 No. 3, pages: 27-37. https://doi.org/10.29332/ijpse.v4n3.611

No competing interests reported.

Download PDF

Submission checks completed at journal
17 Oct, 2024
First submitted to journal
07 Oct, 2024

You are reading this latest preprint version

Soil organic carbon and its relation with the reflectance spectroscopy in soil sampling by layers and compound in the Depositional Piedmont of Villavicencio (Colombia)

Status:

Version 1

Abstract

Figures

Introduction

Materials and Methods

Study area

Geology, geomorphology and soils

Climate

Soil sampling and laboratory analysis

Acquisition and pre-processing of NIR spectra

Statistical analysis

Results

Discussion

Conclusions

Declarations

References

Additional Declarations

Status:

Version 1