Can a history of crop rotations improve the prediction of soil organic carbon in the Andes? integrating machine learning multi-annual crop classification as a proxy of soil management

Figures
Abstract
Soil organic carbon (SOC) is a crucial component related to various processes that ensure soil health and function. Its modeling is vital for assessing and monitoring soil degradation caused by the potential impact of agricultural activities.This study aimed to model SOC in the Northern highlands of Peru, characterized by a high amount of SOC, which is being affected by crop expansion. Crop rotation (CR) history was linked to ground-truthed soil data via a multi-year crop classification model trained on data from 534 fields across 2022–2024. Each cropland field was represented as a polygon delineating its boundaries and indicating its dominant crop cover. Time series of multispectral Sentinel-2 Level-2A Top of Canopy imagery were used to derive phenological features—such as the timing of maximum canopy cover and the length of the growing period—based on Normalized Difference Vegetation Index (NDVI) time series. A Random Forest classifier was used as the baseline model. The cropland classification model demonstrated a strong overall performance, with F1 scores ranging from 0.81 to 0.98 across the different classes. The model performed well for lupin and pasture but scored lower for beans and potatoes. Predictions of cropland classes from 2019 to 2022 were created, resulting in frequency layers that represent crop rotations. Four feature configurations were evaluated: (i) including all features as a benchmark, (ii) excluding climatology, (iii) excluding crop rotation history, and (iv) excluding soil properties. Configurations including all features and excluding crop rotation history showed the highest performance (R2 = 0.63), while those excluding climatology or soil properties performed worse (). Although soil features were the most important, fallow frequency emerged as the most critical predictor of SOC in crop rotations. When soil data were excluded, fallow frequency, combined with climatic features, explained over half of the SOC variability. The findings emphasize the importance of incorporating remote sensing-derived CR into SOC mapping efforts.
Citation: Bueno M, Loayza H, Ninanya J, Rinza J, Briceño P, Silva L, et al. (2026) Can a history of crop rotations improve the prediction of soil organic carbon in the Andes? integrating machine learning multi-annual crop classification as a proxy of soil management. PLoS One 21(7): e0353966. https://doi.org/10.1371/journal.pone.0353966
Editor: RunGuo Zang, Chinese Academy of Forestry, CHINA
Received: October 15, 2025; Accepted: July 1, 2026; Published: July 16, 2026
Copyright: © 2026 Bueno et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All data required to replicate the findings of this study are publicly available. The Sentinel-2 time-series data cube (2019–2024) and associated crop labels are available through Zenodo at: https://zenodo.org/records/19949051 (DOI: 10.5281/zenodo.19949051). The soil dataset supporting this study, including soil physical and chemical properties and soil organic carbon (SOC) measurements for the 189 sampled fields in the Northern Highlands of Peru, is publicly available through the International Potato Center (CIP) Dataverse at: https://doi.org/10.21223/P3/RHE9QI. The source code and example data required to reproduce the workflow are available on GitHub at: https://github.com/kundun14/multi-annual-crop-classification-with-phenological-features.
Funding: This work was supported, in whole or in part, by the Bill and Melinda Gates Foundation (BMGF; grant number OPP1070785), the CGIAR Sustainable Farming Science Program, and the Deutsche Gesellschaft für Internationale Zusammenarbeit (GIZ) through the fund for the promotion of innovation in Agriculture (i4Ag; Agreement No. 81275993). CGIAR is a global research partnership for a food-secure future dedicated to transforming food, land, and water systems in a climate crisis. •The CGIAR Sustainable Farming Science Program provided support for laboratory analyses, as well as personnel salaries and time. •The Deutsche Gesellschaft für Internationale Zusammenarbeit (GIZ), through the i4Ag fund (Agreement No. 81275993), supported field travel and contributed to personnel time. •The Bill and Melinda Gates Foundation (BMGF; grant number OPP1070785) provided support through equipment acquired during the implementation of the project, which was subsequently used in this study, including GPS devices, RGB cameras, and field sampling tools.
Competing interests: The authors have declared that no competing interests exist.
Introduction
To meet the growing global demand for food, agricultural crop production has intensified in recent decades. This intensification threatens soil quality and increases atmospheric greenhouse gas levels [1]. Soil organic carbon (SOC) is widely recognized as a crucial component of soil quality in agricultural systems, playing a significant role in sustaining crop production and ensuring food security. Furthermore, sequestering SOC is considered a promising method to offset greenhouse gas emissions and facilitate carbon neutrality [2]. However, intensive conventional agricultural practices have damaged soil quality, particularly by reducing SOC. To address these issues, it is crucial to enhance SOC accumulation in agroecosystems through the adoption of sustainable agricultural management practices. Crop rotation (CR) is one such practice in regenerative agriculture [3], where a series of crops are grown sequentially over time on the same croplands [4,5]. The adoption of CR can enhance agricultural production and soil fertility [6]. Additionally, CR significantly impacts SOC dynamics, as crop species diversity in CR systems directly affects the quantity and quality of residue-derived carbon input into the soil [7]. Andean farming systems are typically mixed, combining livestock with a high level of inter- and intra-specific crop diversity. These systems primarily rely on native crops such as potato (Solanum spp.), oca (Oxalis tuberosa), olluco (Ullucus tuberosus), mashua (Tropaeolum tuberosum), maca (Lepidium meyenii), and quinoa (Chenopodium quinoa), as well as non-native crops like barley (Hordeum vulgare) and oats (Avena sativa) [8]. Different cropping rotations in the Andes are often accompanied by sectoral fallows, known locally as “laymis”, within communal territories [9]. These fallow periods help restore soil nutrients [10]. However, due to socioeconomic factors and population growth, these regions’ fallow periods have been reduced [11]. The encroachment of crops into highland grasslands has increased in the Andes, primarily driven by global warming [12,13]. This trend poses a significant threat to the SOC stocks preserved by these grasslands [14]; thus, urgent restoration activities are needed to recover these ecosystems [15]. To enhance the quantification of SOC and analyze the key activities that help preserve soil carbon, it is crucial to focus on the Andes, which holds great potential for carbon sequestration [12,16].
However, the degree and variability of SOC accumulation in response to CR, along with the underlying factors driving this response, remain uncertain [17–19]. This uncertainty stems from the complex interactions among influential factors—such as climate, soil properties, and agronomic practices—that often make predicting SOC responses to individual variables difficult [17]. In specific contexts, inherent environmental variables and specific soil covariates may exert a significantly stronger influence on SOC stock variations than crop management history [17,20,21]. Nevertheless, the reported importance of soil properties varies across scales [17]. While climate covariates often outweigh soil predictors in broad regional or global studies like SoilGrids [22–24], local-scale studies highlight the dominant explanatory role of soil texture. Because carbon accumulation is ultimately modulated by the interplay between clay concentration, land use, and climate [25,26], site-specific analyses are essential to determine whether management-driven signals can be identified above the dominant influence of inherent soil characteristics.
Information about CR can be gathered from multi-annual crop inventories, government maps, surveys, and, more recently, earth observation (EO) technologies. Using satellite image time series has emerged as a promising approach for crop type classification in remote sensing. It captures spectral-temporal profiles from multispectral images taken throughout the growing season, reflecting crop phenology and seasonal variation. This approach enhances classification accuracy by extracting unique spectral-temporal features for various crops [27–29]. However, a notable gap remains in knowledge regarding the integration of soil management into agricultural remote sensing frameworks.
This study aimed to assess the importance of crop rotation features in predicting SOC. Multi-year crop rotation covariates were developed to evaluate their contribution to SOC variability and to provide insights into cropping practices that support sustainable soil management. The specific objectives of the study were to: i) develop a phenological-based multi-year (2019–2023) cropland classification model integrating Sentinel-2 time series; ii) evaluate the predictive power of remote sensing-derived CR features for SOC estimation against traditional climatology and soil properties and assess the relative importance of these feature groups. The working hypothesis was that spatially explicit CR features, captured through multi-temporal classification would significantly enhance SOC prediction in heterogeneous Andean landscapes—offering a viable modeling alternative even in the absence of ground-based soil features.
Materials and methods
Study area
The study was conducted in the highlands of the Chugay district (), La Libertad region of Peru, focusing on a study area of approximately 255 km2 (Fig 1A). This area is dominated by grasslands, where smallholder agriculture has been established, with potatoes as the main crop [30]. Potatoes are typically grown under two rotation schemes: (i) fallow – potato – cereal (oat, barley, wheat) – legume (lupin, bean), and (ii) fallow – potato – cereal – fallow.
The study area (indicated by a rectangle over a Sentinel-2 scene from 17/08/2019) is situated in the Chugay District, La Libertad Department, Peru (A). Red triangles denote soil organic carbon sampling points, while green patches depict cropland fields in the study area, shown over a Digital Elevation Model from ALOS PALSAR data distributed by Alaska Satellite Facility (B). Monthly mean minimum (blue) and maximum (red) temperatures, and monthly accumulated precipitation (sky-blue bars) for the 1981–2016 period, derived from PISCO datasets (C).
Fallow refers to arable land left uncultivated to recover from intensive cropping, typically covered by low-density Rumex acetosella weeds. Potatoes are primarily cultivated during the rainy season, from October to April of the following year. To better understand the cropping dynamics in this region, it is important to distinguish between the main (MS) and secondary (SS) growing seasons. MS aligns with the onset of rains from late September to early October, during which potatoes are the predominant crop. SS may occur between May and September, although its timing is highly variable due to fluctuations in soil moisture availability. As a result, the MS spans two calendar years, while the SS falls within a single year, reflecting both inter- and intra-annual cropping phenology.
The climate is rainy and temperate with dry autumns and winters [31]. Monthly minimum and maximum temperatures range from 5.4 to 8.0 ◦C and from 18.5 to 21.5 ◦C, respectively, with an average annual rainfall of 757.7 mm, approximately 90% of which occurs during the potato growing season (Fig 1C), based on the 1981–2016 period from PISCO datasets [32].
According to the national soil map of Peru [33], the study area is predominantly characterized by Eutric Regosols–Eutric Cambisols (RGe–CMe code; i.e., young to moderately developed soils with slightly acidic to neutral pH and high surface organic matter), with smaller areas of Eutric Leptosols–Lithic outcrops (LPe–R code; i.e., shallow, somewhat acidic soils with coarse fragments near the surface) in the north and south. These soils are suitable for grazing or pastures, but present limitations for other crops due to edaphic constraints [34], as they are naturally acidic and lack carbonates despite their high organic matter content [35].
Cropland labels and digitization of prediction masks
Cropland cover information of 189 fields was collected in 2022 and 2024 through field visits [36] and farmer surveys with prior authorization from the District Municipality of Chugay, facilitated by Asociación Pataz https://asociacionpataz.org.pe/, an organization whose mission is to promote sustainable development opportunities in vulnerable communities.
Planting and harvesting dates for most fields in 2022 and 2024 were validated directly in the field through direct observation and surveys conducted with farmers. The coordinates of the central point of each field were obtained using a handheld GPS (Fig 1B). This was essential for spatially aligning the crop monitoring points with farmers’ fields for the soil sampling and the subsequent remote sensing analysis.
The field boundaries were determined using central point GPS coordinates and Planet imagery with a 1.5 m ground sample distance (GSD), acquired in August 2022 [37], ensuring synchronized field visits for crop monitoring. Each cropland was represented as a polygon delineating the field boundaries and its dominant crop type coverage (fallow lands were treated as an additional modeling category) from 2022 to 2024, for both the MS and SS. To increase the number of samples for cropland classification, a dataset of 345 potato (Solanum tuberosum L.) fields containing only sowing and harvest information from the 2021 MS was included from the e-Agrology platform database [36]. All non-agricultural land covers (e.g., water bodies, forests, and urban areas) were omitted because SOC sampling was conducted exclusively in agricultural areas.
In addition to these 534 monitored cropland fields used in the subsequent cropland type classification modeling, most of the remaining croplands in the study area were identified and digitized, resulting in an additional 1,843 cropland fields. The digitization process used high-resolution imagery [37] for precise delineation. Finally, this dataset was used to create a comprehensive cropland and non-cropland mask to enable spatially continuous predictions across the study area.
Cropland label processing
The cropping calendar in the study area is highly variable due to the region’s diverse climate conditions. Therefore, to support season-specific modeling, the Sentinel-2 (S2) imagery time series was grouped by season type into MS and SS, enabling separate analysis of interannual and intra-annual cropping periods. To address class imbalance, only the most frequent crop classes, representing 97% of the dataset, were kept while less frequent crops (i.e., quinoa—Chenopodium quinoa—and corn—Zea mays—) were dropped.
Additionally, fields that were too small () or excessively narrow, as identified in the original database, were removed to minimize mixed pixels and facilitate classification following [38] recommendations. Finally, the croplands considered were: barley (Hordeum vulgare L.), oat (Avena sativa L.), pasture, lupin (Lupinus mutabilis), beans (Vicia faba), potato (Solanum tuberosum L.), and fallow. Thus, fallow was included as a core component of the CR definition, as higher frequencies of covered fallow periods are generally associated with increased SOC contents [11,39].
These cropland labels generally span several seasons within the same field, thus providing a reliable multi-season temporal record.
In this study, an instance is defined as a uniquely labeled cropland field for a specific season (therefore, a single field may have multiple instances—for example, potato MS 2023, barley SS 2024, etc.). These instances served as the basis for the subsequent S2 time series acquisition and phenological analysis, which provided the input data for the cropland type classification modeling.
Sentinel-2 time series acquisition
A custom algorithm was developed to download Harmonized Level-2A surface reflectance time series for each cropland using the Google Earth Engine (GEE) Python API and geemap [40]. The data were retrieved from the COPERNICUS/S2_HARMONIZED ee.ImageCollection, which provides 10 m GSD visible and near-infrared (NIR) bands and 20 m GSD shortwave infrared (SWIR) bands. Scenes with >70% cloud cover were excluded using the Scene Classification Layer metadata, while residual cloudy pixels were subsequently masked using a cloud probability layer [41]. Data were downloaded from 2019 to 2024; however, for 2024, data were only collected up to September, before the start of a new seasonal cycle, which was not considered in this study. The S2 SWIR bands were resampled to 10 m GSD using the Python library Xarray version 2023.3.0 [42]. Subsequently, reflectance time series were smoothed using the Savitzky–Golay filter implemented in the signal R package [43]. The resulting processed data, which provides the time-series data cube of Sentinel-2 from 2019 to 2024 imagery and crop labels, is openly accessible at https://zenodo.org/records/19949051.
Since planting and harvesting dates identified through field surveys were not always reliable or available, crop growth windows were estimated using a consistent criterion. The zonal medians of NDVI [44], derived from red and NIR reflectance (see S1 Table), and their temporal patterns were analyzed for each cropland field, using medians instead of means to reduce the influence of outliers, given that the median is more robust [45]. Because NDVI reflects vegetation vigor and biomass, higher values are typically observed during peak growth stages (e.g., flowering), while values near zero are indicative of planting and senescence dates [46]. Finally, plausible planting and harvesting dates were estimated through preliminary surveys and then refined using local knowledge of cropping calendars, rainfall onset (e.g., in Andean rainfed agriculture, planting dates are determined by rainfall patterns), and visual inspection of NDVI time series.
Sentinel-2 time series cropland phenological analysis
Achieving crop separability during each growing season involved combining three types of features that were used for multi-year cropland type classification modeling (S1 Table): (i) phenological characteristics derived from NDVI time series, (ii) vegetation indices, and (iii) preliminary clustering of NDVI time series of cropland types using Dynamic Time Warping (DTW).
First, a Gaussian function was used to model crop seasonal phenology, fitting NDVI time series in a way that preserves the maximum amount of information while simultaneously reducing the dimensionality of the data [47], which was parameterized using the following equation:
(1)here is the NDVI value at peak canopy cover time, is the time to reach maximum canopy cover, and is the full width at half maximum of the fitted curve, related to the growing period length. Months After Planting (MAP) tracks crop development time. Only cropland fields with at least eight valid Sentinel-2 scenes were retained. The Gaussian function was fitted using the gslnls library in R [48] by non-linear least squares using the Marquardt method [49]. The resulting parameters—, and —were used as features for cropland classification.
To test for significant differences among cropland covers for the Gaussian parameters, the Kruskal–Wallis rank sum test was used, followed by Dunn’s post hoc test with Holm’s correction for multiple comparisons. The analyses were performed using the R packages stats [48] and DescTools [50]. Additionally, box plots were used to present the distribution of fitted Gaussian parameters by cropland type.
Second, the zonal medians from vegetation indices—bare soil index (BSI), specific leaf area vegetation index (SLAVI), normalized burn ratio 2 (NBR2), and normalized difference moisture index (NDMI) (see S1 Table)—were evaluated at the time of maximum canopy cover, i.e., closest to for each cropland type. Vegetation index calculations and zonal medians were performed in Python using Xarray version 2023.3.0.
Third, Dynamic Time Warping (DTW) clustering [51] was used to group similar NDVI growth patterns. DTW uses the Euclidean distance for similarity measurement between two temporal sequences. Given two time series sequences and , the base distance matrix (where each element stores the Euclidean distances between points in the two sequences) is defined as:
(2)The DTW distance matrix is obtained recursively as:
(3)The lower the distance between two time series, the more similar they are. These distances were used to perform unsupervised clustering of cropland type temporal patterns using the dtwclust R library [52]. The distance parameter was set to dtw_basic, and partitional clustering was applied with k = 7. Fuzzy clustering was also performed using the TSclust library [53], with clusters set to k = 7 in both cases to match the number of cropland types under study.
Multi-year cropland type classification modeling
Given the sample size and label imbalance, a single multiannual model was trained using instances from 2021–2024, following the approach of [54] and [55]. Multiannual training helps capture time-independent features, averages out meteorological variations, and mitigates the impact of class imbalances by increasing the dataset size. To predict the 2019–2022 CR using a model trained on the 2022–2024 field monitoring period without introducing temporal bias, we incorporated seasonal (MS, SS) indicator features into the model, as described in the labeling processing section. These features used windows defined by WorldClim precipitation averages and observed seasonal sowing and harvesting patterns. By establishing these consistent seasonal boundaries, the model could reliably project seasonal crop classifications back to the 2019–2022 period.
The Random Forest (RF) algorithm [56] was used as the modeling framework due to its proven effectiveness for cropland classification with limited labeled data [29,57].
To further address class imbalance in the dataset, data augmentation was applied using the themis library in R with the SMOTE–ENC algorithm [58], setting k = 5 and over_ratio = 1. Before augmentation, categorical features such as the DTW cluster groups were transformed into dummy variables using the fastDummies package.
Hyperparameter tuning was conducted using the tuneRanger package [59], which performs automatic model-based optimization [60] through the ranger package [61]. The tuned parameters included mtry (number of variables randomly sampled at each split), min.node.size (minimum number of samples in a leaf node), and sample.fraction (fraction of observations used for each tree), while the number of trees (num.trees) was fixed at 1000. The tuning process relied on out-of-bag accuracy as the performance criterion, allowing for faster yet robust estimation compared to cross-validation. Finally, with the optimal model configuration, the model was refitted using the optimized parameters on a 70% split of the complete dataset, following the recommendations of [62].
A confusion matrix [63] was used to summarize the classification results on the remaining 30% test set. Because the dataset was highly imbalanced, overall accuracy could be biased toward dominant cropland types (e.g., potato). To address this, the F1 score [64] was calculated for each class. Overall accuracy was also reported for comparison. Predictions were then performed for all cropland fields without an observed crop history.
A working example of the procedure suggested in this publication, including source code and example data, is provided on GitHub (https://github.com/kundun14/multi-annual-crop-classification-with-phenological-features).
Multi-year cropland frequency as a rotation feature
To integrate the spatially explicit, seasonally predicted crop type into an extended CR (2019–2022) dataset that could be meaningfully linked to SOC, we derived a ”cropping frequency” feature for each crop field. Integrating crop sequence information facilitated the incorporation of CR history into our SOC modeling framework. Consequently, cropping frequency was derived from the extended CR to quantify the prevalence of each crop type throughout the 2019–2022 campaigns. This approach enabled the integration of multi-year remote sensing classifications as a predictive feature for measured SOC during the 2022 sampling campaign.
The cropping frequency feature was derived according to the following Eq 4:
(4)where is the cropping frequency (%) for crop i, is the number of seasons in which crop i was detected in a given field, and is the total number of seasons analyzed.
These frequency features reflect the dominance of various crop types over time and may help explain spatio-temporal trends in SOC [65]. This methodology resulted in seven frequency raster layers, one dedicated to each crop type.
Soil sampling
In August 2022, 189 locations (Fig 1B) were sampled using composite sampling of five points centered in the middle of the farmers’ fields, following the procedure of [15]. Soil samples were taken at 0–24 cm depth for physical-chemical analysis and at 15 cm for dry bulk density (Bd) analysis with the core method [66], which involved extracting two samples from the center of the croplands (0.15 and 0.30 m deep) in two volumetric cylinders with diameters of 50 mm and heights of 50.8 mm.
Soil samples were analyzed in 2023 at the Soil Laboratory of the La Molina National Agrarian University (Lima, Peru), following protocols described by [67]. The measured variables included pH (1:1 soil-water suspension), electrical conductivity (EC) (1:1 extract), available phosphorus (P) (Olsen method), available potassium (K) (buffered acetate extraction), soil organic carbon (SOC) using the Walkley & Black method, cation exchange capacity (CEC) determined with buffered ammonium acetate at pH 7.0, exchangeable cations–calcium (eCa), magnesium (eMg), potassium (eK), and sodium (eNa)—and exchangeable acidity (eAC) analyzed by atomic absorption spectrophotometry.
Particle size distribution (sand, silt, and clay) was determined using the Bouyoucos hydrometer method, and exchangeable acidity was also assessed. Additionally, the nitrogen isotopic composition () of dried soil samples was analyzed at the Stable Isotope Facility, University of California, Davis [68]. The soil dataset supporting the findings of this study, which includes the physical-chemical properties and soil organic carbon (SOC) measurements for the 189 sampled fields in the Northern Highlands of Peru, is publicly available through the International Potato Center (CIP) Dataverse at: https://doi.org/10.21223/P3/RHE9QI.
Features for soil organic carbon content prediction
The features used to predict SOC were classified into three groups (S2 Table): (i) soil physical-chemical properties, (ii) gridded historical climate, and (iii) multi-year cropland frequencies.
For the second set of features, monthly temperature, precipitation accumulation, and solar radiation spanning 1961–2000 were acquired from the WorldClim dataset at 1 km2 (30 arc sec) spatial resolution (https://www.worldclim.org, accessed January 20, 2025).
In the case of the cropland frequency features, we followed the approach introduced by [69] and described in the previous sections.
Soil organic carbon modeling
An RF model was trained using 72 soil, climatic, and agronomic features (see S2 Table). A variable-importance analysis was used to identify the most influential features impacting SOC. Following [70], recursive feature elimination [71] was used to select essential predictive features of SOC. This is an iterative process that ranks significant predictors and removes nonexplanatory ones. The process continues until a convergence criterion is met, using the coefficient of determination (R2) as the evaluation metric.
Additionally, a sensitivity analysis was conducted by grouping rotation features, climatology-based features, and soil properties separately to assess their influence on SOC prediction. Four scenarios were explored: (i) including all 72 features as a benchmark (Case 1); (ii) removing climatology, yielding 24 features (Case 2); (iii) omitting crop rotation history, leaving 64 features (Case 3); and (iv) excluding soil properties, resulting in 57 features (Case 4).
Model performance was assessed using k-fold cross-validation, where the dataset was split into 5 folds, each fold used once for validation while the remaining folds were used for training. The R2 and root mean square error (RMSE) were used to evaluate the ability of the Machine Learning model to predict observed SOC. R2 represents the proportion of variance explained by the model, and RMSE indicates the accuracy of the predicted values [72]. R2 and RMSE were calculated as follows:
(5)(6)where , , and N are the predicted SOC, observed SOC, and total number of SOC sample values, respectively. A higher R2 (close to 1) and lower RMSE (close to 0) indicate better model performance.
Partial Dependence Plots (PDP) were used to visualize the relationship between features and SOC predictions, showing whether the relationship is linear, monotonic, or complex. For this purpose, the implementation available in the pdp library [73] was used within the R environment [48].
Results
Phenological-based multi-year cropland classification model
A total of 788 observations were analyzed, with 49.7%, 13.8%, 10.1%, 7.9%, 6.9%, 6.3%, and 5.2% representing potato, fallow, barley, oat, beans, pasture, and lupin, respectively (see Fig 2A). Most fields are small, clustered below 0.5 ha, with a strongly right-skewed distribution, indicating that while few larger fields exist, they are rare, reaching up to about 8 ha (Fig 2B).
(A) Number of fields for each cropland type per crop season recorded via field farmers’ visits and surveys—main season (MS) and secondary season (SS). (B) Field size distribution of farmers’ fields where SOC was measured.
The adjustment of the Gaussian function showed an average R2 of . The best average performance showed the following pattern by cropland: beans ()> oats ()> barley ()> potato ()> fallow ()> tarwi ()> pasture () (Fig 3). Models fitted better in the secondary campaign for most crops (except for lupin) than the main campaign, with a global R2 average of and , respectively (Fig 3).
Based on the coefficient of determination (R2) by crop in the main (MS) and secondary (SS) seasons, 2021–2023. Significance: ns (); * (); ** (); *** (p < 0.001). Symbols indicate statistical significance levels.
values ranged between 0.55 and 0.80, beans () and fallow () showed the highest and lowest global average value in this parameter (S1 Fig), with significant differences among beans vs. fallow–potato, and fallow vs. barley (S3 Table). values ranged between 4 and 6.5; beans () and barley () showed the highest and lowest global average value (S1 Fig), with significant differences among barley vs. all croplands (S3 Table). values ranged between 3.2 and 4.6; pasture () and barley () showed the highest and lowest global average value (S1 Fig), with significant differences among barley vs. all croplands and pasture vs. potato (S3 Table). Crops showed different patterns of clustering using DTW on S2 NDVI-based time series (S2 Fig), with clusters 2, 3, 6, and 7 being the most important in crop differentiation.
The optimal RF architecture for the multi-year cropland classification model was found when , , and . Validation on the test sample set (30% of the full dataset, n = 788 instances) shows strong overall model performance, with F1 scores ranging from 0.81 to 0.98 for the cropland classes (S3 Fig). The model performance ranking was: lupin (0.98)> pasture (0.96)> oat (0.93)> barley (0.92)> fallow (0.91)> beans (0.90)> potato (0.81). The overall accuracy and Kappa Score were 0.91 and 0.90, respectively (S3 Fig). The resulting classification outputs enabled the identification of spatial and temporal crop rotation patterns across MS and SS during the 2019–2023 study period (Fig 4).
Spatial distribution of cropland classification outputs across main and secondary growing seasons in a subarea of the study region, overlaid on Sentinel-2 RGB imagery from the Copernicus Programme. Colors represent the classified crop types and illustrate field-level crop rotation dynamics between 2019 and 2023.
Soil, climate, and crop rotation features for SOC prediction
SOC values ranged from 2.4 to 22.6 %, with a mean value of 10.61 ± 2.95 % (Table 1). Soils were predominantly loam (34.9 ± 6.88%, 49.4 ± 6.43%, and 24.16 ± 7.70% of Sa, Si, and Cy, respectively), with a moderately acidic pH (4.68 ± 0.85), relatively low EC (0.20 ± 0.24 dS m-1), and moderate to high CEC (15.54 ± 6.81 kg-1) (Table 1). Available P and K showed high variability, ranging from 3.1 to 237.4 mg kg-1 and from 19.0 to 950.0 mg kg-1, respectively (Table 1). Mean values of exchangeable bases were 5.27 ± 5.77, 0.60 ± 0.40, 0.42 ± 0.32, and 0.03 ± 0.04 kg-1 for eCa, eMg, eK, and eNa, respectively, while eAC averaged 1.68 ± 1.40 kg-1 (Table 1). Mean Bd and N were 1.09 ± 0.16 g cm-3 (range: 0.79–1.57 g cm-3) and 5.95 ± 0.77 ‰ (range: 3.74–7.97 ‰), respectively (Table 1).
Based on historical climate data, monthly averages varied as follows: air temperature ranging from 3.5–16.9 ° C, solar radiation from 13,400–15,580 kJ m-2 day-1, and precipitation from 22.9–154 mm. Regarding CR features, the cropland cover frequencies from 2019 to 2022 showed the following distribution: barley (38.7%)> potato (21.8%)> beans (12.8%)> pasture (8.6%)> fallow (6.6%)> oats (5.9%)> lupin (5.7%) (Fig 5). Barley was the most dominant crop, covering over 50% of seasons in 73 fields, while potato and beans only dominated in 5 and 2 fields, respectively, and the remaining cropland types rarely dominated any parcel. The least persistent crops, appearing in fewer than 10% of seasons in most fields, were sorted in descending order as follows: oats, fallow, lupin, pasture, and beans.
Spatial distribution of the features generated by the multi-year cropland classification model from 2019 to 2022, including both surveyed and digitized fields for the study area croplands, overlaid on Sentinel-2 imagery from the European Union Copernicus Programme.
Among the four feature-set configurations, the SOC predictive model performed best in Cases 1 and 3, as evidenced by higher R2 values (0.63 in both cases) and lower RMSE values (1.99% and 1.98%, respectively). Conversely, Cases 2 and 4 exhibited comparatively poorer performance, with lower R2 values (0.53 and 0.52, respectively) and higher RMSE values (2.24% and 2.28%, respectively) (Fig 6).
SOC values for four modeling cases: (1) all features included; (2) without climatic features; (3) without rotation features; and (4) without soil features. Solid red lines represent linear regressions between observed and predicted SOC. The dashed line indicates the 1:1 relationship. The regression equation, adjusted determination coefficient (R2), and Root Mean Square Error (RMSE) are shown for each case.
In general, the SOC predictive model showed a tendency toward overestimation for values below ∼10%. In contrast, for values above it, the model showed a tendency toward underestimation (Fig 6). In Cases 1–3, which included soil features, these variables were among the most important predictors of SOC, with Cy being the most relevant and, together with CEC, EC, eAC, eCa, and eK, consistently ranking among the top 5–7 features (Fig 7). In contrast, soil pH and the other textural fractions (Sa and Si) showed much lower importance, falling below the top 15 in Case 2 and below the top 20 in Cases 1 and 3 (Fig 7). In general, CR features ranked higher than climatic features, with fallow frequency being the most critical CR feature (Fig 7). It ranked first in Case 4 (absence of soil features), while in Case 1 (together with bean frequency) and Case 2, it ranked within the top 9 and 6 respectively (Fig 7). Potato frequency ranked within the top 20 only in Case 2 (Fig 7). Among the climatic features, precipitation variables consistently ranked higher than solar radiation. Notably, the average precipitation in January was the most significant, ranking within the top 11, 9, and 2 in Cases 1, 3, and 4, respectively. In contrast, temperature variables (both mean and maximum) did not appear among the top 20 in any case (Fig 7). In Case 4, 19 of the top 20 features corresponded to climatic features (Fig 7).
Comparison across four modeling scenarios: Case 1 (all features), Case 2 (without climatic features), Case 3 (without rotation features), and Case 4 (without soil features). Variables are color-coded as green (climate features), purple (rotation features), and orange (soil features). The bars represent the normalized importance of each variable in the model. The soil properties feature abbreviations are described in S2 Table.
The PDP analysis revealed an inverse relationship between Cy and predicted SOC, with sensitivity to SOC variations declining above 20% Cy in Cases 1 and 3, and above 40% Cy in Case 2 (Fig 8). In contrast, CEC, EC, and eAC exhibited a direct proportional relationship with predicted SOC (Fig 8). SOC changes in response to CEC were notably faster below 12 kg-1 than at higher values in Case 1 and Case 2. In contrast, SOC responded more slowly to eAC in Case 3 compared to Case 2 (Fig 8). Fallow frequency exhibited a “V”–shaped pattern‘with predicted SOC, showing an inverse relationship below 25% and a direct relationship above 25% in Case 4 (Fig 8). Monthly precipitation (January and May) showed a direct relationship with predicted SOC across specific ranges, with minimal sensitivity observed below 115 mm for January and below 40 mm or above 49 mm for May in Case 4 (Fig 8).
SOC as a function of top three most important variables for each case: clay content (%), cation exchange capacity ( kg-1), electrical conductivity (dS m-1), exchangeable acidity ( kg-1), fallow frequency, and monthly precipitation in January and May (mm). Shaded areas show confidence intervals. Line and ribbon colors indicate the feature category: orange for soil variables, purple for rotation variables, and green for climatic variables. The dashed line represents the mean of observed SOC values.
Discussion
4.1. Phenology-informed features enable accurate crop classification in heterogeneous Andean smallholder system
In this study, phenological trajectories were modeled with a simple symmetric Gaussian function, parameterized by planting and harvesting dates obtained from field surveys. Despite its parsimony, it requires only three coefficients compared to more complex functions such as the asymmetric Gaussian (four parameters; [74]) and the Double Logistic (five parameters; [75]). The range of R2 values between observed and fitted NDVI across crops and seasons obtained in this study (0.22–0.89, Fig 3) reached levels comparable to Gaussian-based analyses of MODIS-EVI time series in boreal forests (0.89–0.94; [76]) and with barley classification in Australia (0.91–0.96; [47]). Performance was particularly enhanced during the secondary growing season, when reduced rainfall and cloud cover improved the quality of Sentinel-2 imagery (Fig 1C). Among the extracted parameters, the Gaussian mean () proved especially informative, as it reflects the duration of the crop growth cycle. Significant differences in between barley and other crops (S3 Table) underscore its value for distinguishing species with contrasting phenological patterns. NDVI evaluated at also emerged as a critical indicator of biomass and vigor, consistent with prior research demonstrating its utility for yield prediction in wheat, rice, and cotton [77–79].
The RF classification model achieved strong overall performance, with an accuracy of 0.91 and a Kappa coefficient of 0.90 across seven crop categories. Similarly, RF has proved highly competitive against more complex deep learning models (i.e., CLSTM, 3D U-Net) simulating cloud-prone croplands belonging to smallholders in Africa, particularly when labeled data were limited [29]. This resonates with the study area context, where frequent cloud cover poses similar challenges; yet, RF demonstrated robustness and adaptability across diverse crop categories, even with scarce labeled field data (F1 of 0.71 to 0.93), underscoring the adaptability of time series approaches in smallholder contexts. Nonetheless, classification accuracy varied by crop, reflecting differences in spectral distinctiveness and exposure to seasonal constraints. Lupin and pasture (F1 = 0.98 and 0.96, respectively) yielded the highest accuracies, confirming the ability of Sentinel-2 time series to capture strongly contrasting growth patterns. The performance achieved for lupin is particularly noteworthy, as this study, to our knowledge, represents one of the first documented cases of its reliable discrimination from spaceborne optical imagery. Similarly, the accuracy of pasture classification aligns with studies that successfully integrated Sentinel-1 and −2 data [80] or multi-level neural network approaches [81], reinforcing the conclusion that temporal optical data are sufficient for identifying this crop type.
Other categories also displayed high performance: oat (0.93), barley (0.92), and fallow (0.91). These findings echo recent evidence demonstrating the separability of these classes using temporal metrics. For instance, oats have been reliably identified at the parcel level using Sentinel-2 imagery [82,83], while barley has been differentiated from wheat using spectral features and classical classifiers [84]. Likewise, fallow land has been effectively mapped using multi-sensor time series and phenological curve fitting, with F1 scores above 0.95 [85]. That such accuracies were obtained in the present study, using a considerably smaller training dataset (224 ha versus more than 2000 ha in comparable work), underscores the robustness of the approach. By contrast, beans (0.90) and potatoes (0.81) exhibited lower classification accuracies. This decline is consistent with previous findings that highlight challenges in separating beans from other legumes and cereals [86]. Potato discrimination is minimal when cultivation coincides with the rainy season (Fig 1C), as cloud cover and incomplete time series reduce the availability of reliable spectral data [87]. These limitations highlight that classification accuracy can be affected when crop cycles overlap with challenging atmospheric conditions, even with robust phenology-based approaches, and to some extent because of survey data uncertainty.
4.2. In the absence of comprehensive soil measurements, crop rotations provide an effective alternative for modeling SOC
Clay content exhibited a robust and inverse relationship with SOC. Partial dependence analysis revealed that SOC declined sharply in soils with clay contents exceeding 20%, while consistently lower SOC levels were also observed when the CEC fell below 10 cmolc kg-1. These findings align with global reviews showing that clay is often the most influential soil predictor of SOC [17] and digital soil mapping studies where gridded clay, sand, and silt covariates explained significant variation in SOC [20,21]. Nevertheless, the reported importance of clay varies across contexts. For instance, [22] found climate covariates outweighed soil predictors in Argentina, while [23] identified soil depth as the dominant feature in Canada. Such divergences reflect scale, data resolution, and environmental conditions. Since global products like SoilGrids interpolate across broad regions [24,88], climatic factors often overshadow soil variables [89]. By contrast, local studies, as in the present work, highlight the strong explanatory role of clay. Mechanistically, clay influences SOC through surface adsorption and aggregate occlusion, protecting organic matter from microbial decomposition [90,91]. However, accumulation also depends on land use and climate [25,26,92]. Thus, while clay concentration is consistently important, its role is modulated by context, reinforcing the need for site-specific analyses.
Although soil features were dominant, fallow frequency emerged as the most critical predictor of SOC among the crop rotation features. When soil data were excluded, fallow frequency, combined with climatic features, explained over half of the SOC variability. Moderate frequencies of fallow were associated with higher SOC, consistent with studies showing that rotation history enhances soil health [65,69,93,94]. However, the statistical dominance of soil properties over CR features in this model may be attributed to two factors. First, soil properties such as clay content are directly related to the physical stabilization of SOC [19], consistently ranking as primary predictors in high-resolution models. Second, while CR impacts SOC dynamics, measurable changes in carbon stocks typically require rotation cycles of 8–10 years to manifest [7,17], a longer horizon than the modeling period. Nevertheless, the positive impact of reduced soil disturbance through permanent cover rather than annual cropping [89] remains evident. These findings emphasize that CR data, readily accessible through remote sensing, can serve as a practical surrogate for costly soil sampling, enabling moderate prediction performance without the need for extensive laboratory analyses.Potato frequency was not correlated with SOC, contrary to the expectation that intensive tillage associated with potatoes would result in reduced carbon. Tillage has been shown to accelerate soil carbon losses in conventional systems [95]. However, Andean potato-based systems often involve less soil disturbance [96], which may mitigate these effects. This highlights the need for context-specific assessments of management practices. Traditional Highland Andean agricultural systems incorporate multiyear fallow periods to ensure soil fertility and the recovery of soil functionality, as a key component of complex CR [11]. However, the increasing food demand from cities and markets has driven a reduction in the length of fallow periods [97]. More broadly, the role of fallow in Andean agriculture remains a complex issue [98].
Thus, in the highlands of the Central Peruvian Andes, where Asian markets promoted the cultivation of maca (Lepidium meyenii), an energetic root, [15] did not find significant differences in SOC among fallows of varying ages after maca cultivation. In this ecosystem, [13] reported higher SOC in maca than in fallow plots caused by the incorporation of manure during the management of this crop. On the other hand, a recent study, [99] pointed out that the seeding of pastures and legumes of traditional fallows has the potential to increase SOC and soil health regeneration in the long term in the Peruvian Central Andes. Such examples highlight the diversity of fallow practices across the Andes, where crop expansion into grasslands and upward shifts in cultivation [100] disrupt the balance between native ecosystems, rotations, and soil carbon storage. Thus, while fallow frequency emerges as a key predictor, its ecological implications require careful contextualization.
Climatic influences add further complexity. Precipitation has consistently been identified as a dominant feature for SOC modeling [22,89], and its interaction with CR is nonlinear. [17] found SOC increases under rotations were strongest in regions with intermediate mean annual temperature (8–15 °C) and precipitation (600–1000 mm), while SOC losses occurred below 120 mm annual rainfall. In the present study area, where annual precipitation ranges from 900 to 1300 mm, high SOC values were obtained (10.61 ± 2.95%, Table 1). This suggests that favorable climatic conditions may mask the effects of management, complicating the detection of SOC gains from sustainable practices [101]. Future modeling should therefore explicitly account for the interplay between climate and management to better resolve the contributions of CR features.
4.3. Limitations and implications
Although the Sentinel-2 time series effectively characterizes crop phenology [102], its usability is limited by cloud cover. In Chugay, 69.9% and 29.6% of annual scenes were cloud-covered during the main and secondary seasons, respectively, despite the five-day revisit time. Incorporating all-weather observations from non-optical platforms, such as Synthetic Aperture Radar, could address this limitation [103]. The straightforward accessibility of these data (e.g., via GEE; [104]) would facilitate easy integration. However, Synthetic Aperture Radar data require proper radiometric slope correction and the provision of invalid data masks over areas affected by layover and shadow, particularly in mountainous Andean regions [105]. Integrating Sentinel-1 and −2 time series shows promise [80], although robust multi-sensor workflows remain complex. Further studies are encouraged to advance a multi-sensor approach while maintaining a physiologically sound foundation. Phenological analysis could benefit from established tools like TIMESAT [106] or R packages such as phenofit [76] rather than custom scripts. Since agricultural management decisions are generally made at the parcel level, object-based, rather than per-pixel analysis, is recommended to reduce pixel noise and limit error propagation [80,107].
While EO information is currently easily accessible, greatly facilitating local and regional studies on crop mapping for sustainable development [108], the most critical gap for Andean crop mapping is reference data. Unlike other regions [29,109,110], the Andes lack comprehensive EO datasets linked with field cropping information. We encourage systematic data collection with quality assessment protocols. Since Andean agriculture is highly complex and parcel sizes are small, grouping similar crops into homogeneous categories can be necessary, as distinguishing between crops with similar spectral signatures or growth cycles is a common source of error. Field monitoring at multiple dates, along with a cropping calendar validated by farmers [111], and strategies for sparse data modeling are essential [57]. Regarding field data quality, methods such as hierarchical clustering or self-organizing maps [112] can help reduce spectral confusion and address data imbalance (e.g., grouping barley and oats).
In the study area, the mean values of SOC measured (10.6%) surpasses the values reported in highland Andean grasslands in Peru (4–6%; [13,15]). This highlights the importance of modeling SOC under this context to distinguish the crucial factors that affect C cycling in this ecosystem. Accurate prediction of SOC in the Andes requires not only technical advances in remote sensing and modeling but also strong coordination among stakeholders to ensure that field data collection and management practices are both scientifically robust and socially viable. The heterogeneity of mountain landscapes–shaped by steep environmental gradients, diverse management systems, and climate variability–makes SOC modeling particularly prone to uncertainty, underscoring the need for joint sampling protocols, standardized databases, and shared monitoring platforms developed by universities, non-governmental organizations, government agencies, and local communities [113–115]. Such collaborative approaches enable the integration of scientific insights by research institutions on practices like crop rotations and fallow periods, while non-governmental organizations and farmers contribute essential knowledge of local feasibility and cultural acceptance, thereby strengthening the credibility of SOC projections.
At the same time, robust monitoring, reporting, and verification protocols are indispensable to link SOC modeling with emerging carbon market mechanisms, where innovations in digital monitoring, reporting and verification, such as remote sensing, big data, and cloud infrastructures, are improving transparency, reducing transaction costs, and aligning local initiatives with international standards [116,117]. However, implementing these systems in Andean contexts faces constraints related to technical capacity, costs, and governance of data ownership, which highlights the need for pilot initiatives that couple multi-stakeholder coordination with open-access monitoring, reporting, and verification infrastructures [118]. Such integration not only reinforces the scientific basis of SOC monitoring but also positions Andean communities to participate in future carbon markets with greater equity and long-term sustainability.
Conclusion
To our knowledge, this study is the first to estimate multi-season, field-scale crop rotations and incorporate them as predictors in SOC modeling within the Andes. The results demonstrate that phenology-driven classification frameworks can reliably discriminate crops even in heterogeneous Andean agroecosystems characterized by small plots and diverse management practices. The integration of Sentinel-2 time series with phenological feature extraction and machine learning modeling effectively addresses many of these challenges. While studies considering fewer crop categories generally report higher accuracy, the successful classification of seven distinct classes in this study demonstrates both the scalability and robustness of the proposed approach. Beyond crop mapping, such reliable discrimination of crop rotations provides critical inputs for downstream SOC modeling.
Overall, the study reveals a hierarchy of predictors for SOC, with clay content exerting the strongest influence, followed by fallow frequency and precipitation. These results suggest that although soil properties remain indispensable, crop rotation and climate variables provide complementary information—particularly in contexts where soil data are limited or costly to obtain. Notably, fallow detection through remote sensing shows strong potential for improving SOC monitoring in Andean landscapes. Nevertheless, the diversity of fallow practices and the confounding influence of climate underscore the need for locally calibrated models. The integration of soil, rotation, and climate features therefore offers a more comprehensive understanding of SOC dynamics, with direct implications for designing sustainable land management strategies in fragile mountain agroecosystems.
Supporting information
S1 Table. Classification Features.
Classification features for multi-year cropland classification using Sentinel-2 imagery A and B time series (S2A, S2B) from the Harmonized Sentinel-2 Multispectral Instrument (MSI). Blue, green, red, near-infrared (NIR), shortwave infrared 1 (SWIR1), and shortwave infrared 2 (SWIR2) are surface reflectance values from Band 2 (blue, 496.6 nm (S2A)/492.1 nm (S2B)), Band 3 (green, 560 nm (S2A)/559 nm (S2B)), Band 4 (red, 664.5 nm (S2A)/665 nm (S2B)), Band 8 (NIR, 835.1 nm (S2A)/833 nm (S2B)), Band 11 (SWIR1, 1613.7 nm (S2A)/1610.4 nm (S2B)), and Band 12 (SWIR2, 2202.4 nm (S2A)/2185.7 nm (S2B)). f = 0.05.
https://doi.org/10.1371/journal.pone.0353966.s001
(DOCX)
S2 Table. Feature sets used in the soil organic carbon sensitivity analysis.
https://doi.org/10.1371/journal.pone.0353966.s002
(DOCX)
S3 Table. Kruskal–Wallis test and mean rank differences of Dunn’s post hoc test results.
Results with Holm’s correction method comparing Gaussian parameters across different cropland cover categories. NDVIMAX is the Normalized Difference Vegetation Index (NDVI) value at peak canopy cover time, is the time to reach peak growth, and is the full width at half maximum, signaling the length of the growing period. The significance codes are: *** for p < 0.001, ** for p < 0.01, * for p < 0.05,· for p < 0.1, and no symbol for .
https://doi.org/10.1371/journal.pone.0353966.s003
(DOCX)
S1 Fig. Boxplots of NDVI-derived phenological parameters.
(A) Maximum Normalized Difference Vegetation Index (NDVI) value at peak canopy cover time (NDVIMAX), (B) time to reach peak growth (), and (C) full width at half maximum (), estimated from Gaussian functions fitted to Sentinel-2 NDVI time series for each cropland type.
https://doi.org/10.1371/journal.pone.0353966.s004
(TIF)
S2 Fig. Classification of NDVI time series using Dynamic Time Warping (DTW).
The width of the violin diagrams represents the contribution of the corresponding cluster for the characterization of each crop.
https://doi.org/10.1371/journal.pone.0353966.s005
(TIF)
S3 Fig. Confusion matrix for the test set.
Test set corresponds to 30% of the full dataset (n = 788 instances).
https://doi.org/10.1371/journal.pone.0353966.s006
(TIF)
Acknowledgments
We sincerely thank all the farmers from the communities in the Chugay district who generously shared information for this study. We also thank the Google Earth Engine team for providing the free platform for Sentinel-2 data processing.
References
- 1. Amelung W, Bossio D, de Vries W, Kögel-Knabner I, Lehmann J, Amundson R, et al. Towards a global-scale soil climate mitigation strategy. Nat Commun. 2020;11(1):5427. pmid:33110065
- 2.
Lal R. World cropland soils as a source or sink for atmospheric carbon. Advances in Agronomy. Elsevier. 2001. p. 145–91. https://doi.org/10.1016/s0065-2113(01)71014-0
- 3. Newton P, Civita N, Frankel-Goldwater L, Bartel K, Johns C. What Is Regenerative Agriculture? A Review of Scholar and Practitioner Definitions Based on Processes and Outcomes. Front Sustain Food Syst. 2020;4.
- 4.
FAO. Food and Agriculture Organization of the United Nations Conservation Agriculture. FAO.org; 2021. https://www.fao.org/conservation-agriculture/en/
- 5. Khangura R, Ferris D, Wagg C, Bowyer J. Regenerative Agriculture—A Literature Review on the Practices and Mechanisms Used to Improve Soil Health. Sustainability. 2023;15(3):2338.
- 6. Qin J, Bian C, Duan S, Wang W, Li G, Jin L. Effects of different rotation cropping systems on potato yield, rhizosphere microbial community and soil biochemical properties. Front Plant Sci. 2022;13:999730. pmid:36247636
- 7. King AE, Blesh J. Crop rotations for increased soil carbon: perenniality as a guiding principle. Ecol Appl. 2018;28(1):249–61. pmid:29112790
- 8. Oswald A, De Haan S, Sanchez J, Ccanto R. The complexity of simple tillage systems. J Agric Sci. 2009;147(4):399–410.
- 9.
Orlove B, Godoy R. Sectoral Fallowing Systems in the Central Andes. In: 1986. https://papers.ssrn.com/abstract=2229453
- 10. Pestalozzi H. Sectoral Fallow Systems and the Management of Soil Fertility: The Rationality of Indigenous Knowledge in the High Andes of Bolivia. Mountain Research and Development. 2000;20(1):64–71.
- 11. de Haan S, Juárez H. Land use and potato genetic resources in Huancavelica, central Peru. Journal of Land Use Science. 2010;5(3):179–95.
- 12. Quiroz R, Ramírez DA, Kroschel J, Andrade-Piedra J, Barreda C, Condori B, et al. Impact of climate change on the potato crop and biodiversity in its center of origin. Open Agriculture. 2018;3(1):273–83.
- 13. Rolando JL, Dubeux JC Jr, Perez W, Ramirez DA, Turin C, Ruiz-Moreno M, et al. Soil organic carbon stocks and fractionation under different land uses in the Peruvian high-Andean Puna. Geoderma. 2017;307:65–72.
- 14. Rolando JL, Turin C, Ramírez DA, Mares V, Monerris J, Quiroz R. Key ecosystem services and ecological intensification of agriculture in the tropical high-Andean Puna as affected by land-use and climate changes. Agriculture, Ecosystems & Environment. 2017;236:221–33.
- 15. Carbajal M, Ramírez DA, Turin C, Schaeffer SM, Konkel J, Ninanya J, et al. From Rangelands to Cropland, Land-Use Change and Its Impact on Soil Organic Carbon Variables in a Peruvian Andean Highlands: A Machine Learning Modeling Approach. Ecosystems. 2024;27(7):899–917. pmid:39524473
- 16. Grados D, García S, Schrevens E. Assessing the potato yield gap in the Peruvian Central Andes. Agricultural Systems. 2020;181:102817.
- 17. Liu X, Tan S, Song X, Wu X, Zhao G, Li S, et al. Response of soil organic carbon content to crop rotation and its controls: A global synthesis. Agriculture, Ecosystems & Environment. 2022;335:108017.
- 18. Zou C, Pearce RC, Grove JH, Coyne MS. Conservation practices in tobacco production increase large aggregates and associated carbon and nitrogen. Soil Science Society of America Journal. 2015;79(6):1760–70.
- 19. Zou C, Li Y, Huang W, Zhao Gaokun, Pu G, Su J, et al. Rotation and manure amendment increase soil macro-aggregates and associated carbon and nitrogen stocks in flue-cured tobacco production. Geoderma. 2018;325:49–58.
- 20. Kellner JR, Clanton C, Demuth KM, Donovan M, Feng YK, Khim-Young M, et al. Digital soil mapping in support of voluntary carbon market programs in agricultural land. PLoS One. 2025;20(9):e0327895. pmid:40892935
- 21. McBratney AB, Mendonça Santos ML, Minasny B. On digital soil mapping. Geoderma. 2003;117(1–2):3–52.
- 22. Heuvelink GBM, Angelini ME, Poggio L, Bai Z, Batjes NH, van den Bosch R, et al. Machine learning in space and time for modelling soil organic carbon change. European J Soil Science. 2020;72(4):1607–23.
- 23. Sothe C, Gonsamo A, Arabian J, Snider J. Large scale mapping of soil organic carbon concentration with 3D machine learning and satellite observations. Geoderma. 2022;405:115402.
- 24. Hengl T, Mendes de Jesus J, Heuvelink GBM, Ruiperez Gonzalez M, Kilibarda M, Blagotić A, et al. SoilGrids250m: Global gridded soil information based on machine learning. PLoS One. 2017;12(2):e0169748. pmid:28207752
- 25. Lamichhane S, Kumar L, Wilson B. Digital soil mapping algorithms and covariates for soil organic carbon mapping and their implications: A review. Geoderma. 2019;352:395–413.
- 26. Wiesmeier M, Urbanski L, Hobley E, Lang B, von Lützow M, Marin-Spiotta E, et al. Soil organic carbon storage as a key function of soils - A review of drivers and indicators at various scales. Geoderma. 2019;333:149–62.
- 27. Garnot VSF, Landrieu L. Lightweight Temporal Self-Attention for Classifying Satellite Image Time Series. 2020.
- 28. Peña MA, Brenning A. Assessing fruit-tree crop classification from Landsat-8 time series for the Maipo Valley, Chile. Remote Sensing of Environment. 2015;171:234–44.
- 29.
Rustowicz R, Cheong R, Wang L, Ermon S, Burke M, Lobell D. Semantic segmentation of crop type in Africa: A novel misc and analysis of deep learning methods. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, 2019. 75–82.
- 30. Pradel W, Suarez V, Ramírez D, Rinza J, Briceño P, Otiniano R, et al. Scaling potential of mulching technologies in potato producers in Peruvian highlands. 2024.
- 31.
Castro A, Davila Arriaga C, Laura W, Cubas Saucedo F, Avalos G, López C. Climas del Perú: mapa de clasificación climática nacional. Servicio Nacional de Meteorología e Hidrología del Perú. 2021. http://repositorio.senamhi.gob.pe/handle/20.500.12542/1336
- 32. Aybar C, Fernández C, Huerta A, Lavado W, Vega F, Felipe-Obando O. Construction of a high-resolution gridded rainfall dataset for Peru from 1981 to the present day. Hydrological Sciences Journal. 2019;65(5):770–85.
- 33.
INRENA. Mapa de suelos del Perú. 1.a ed. Ministerio de agricultura. 1996.
- 34.
Bazán MAF, Luján Rojas JM, Llaque Fernández GI, Valderrama Puscan MW. Application of geospatial analysis for soil classification in La Libertad region, Peru. In: LACCEI 2022 Conference Proceedings, 2022. https://laccei.org/LACCEI2022-BocaRaton/meta/FP121.html
- 35. Briceño P, Ninanya J, Seminario JF, Otiniano R, Rinza J, Mestanza C, et al. Can regenerative agriculture enhance productivity, profitability, and reduce C emissions? A case study in Andean potato farming. European Journal of Agronomy. 2026;175:128020.
- 36.
Rinza J, Diaz E, Segura J, Ramirez D, Ninanya J, Otiniano R. Agronomic log of potato small farmers during the 2021-2022 campaign from Chugay district, La Libertad - Peru. International Potato Center. 2022. https://doi.org/10.21223/MGMJ1T
- 37.
Planet Labs PBC. Planet Application Program Interface: In Space for Life on Earth. https://api.planet.com 2022.
- 38. Hengl T. Finding the right pixel size. Computers & Geosciences. 2006;32(9):1283–98.
- 39. López I Losada R, Hedlund K, Haddaway NR, Sahlin U, Jackson LE, Kätterer T, et al. Synergistic effects of multiple “good agricultural practices” for promoting organic carbon in soils: A systematic review of long-term experiments. Ambio. 2025;54(11):1715–28. pmid:40423915
- 40. Wu Q. geemap: A Python package for interactive mapping with Google Earth Engine. JOSS. 2020;5(51):2305.
- 41.
Schwehr K. Sentinel-2: Cloud Probability in Earth Engine. Zenodo. 2020. https://doi.org/10.5281/zenodo.7411046
- 42. Hoyer S, Hamman J. xarray: N-D labeled Arrays and Datasets in Python. JORS. 2017;5(1):10.
- 43.
Ligges U, Short T, Kienzle P, Billinghurst D, Schnackenberg S, Billinghurst D. Package ‘signal’. 2024.
- 44.
Rouse JW, Haas RH, Well JA, Deering DW. Monitoring vegetation systems in the Great Plains with ERTS. NASA Technical Reports Server (NTRS). 1974. https://ntrs.nasa.gov/citations/19740022614
- 45. Leys C, Ley C, Klein O, Bernard P, Licata L. Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median. Journal of Experimental Social Psychology. 2013;49(4):764–6.
- 46. Reuss J, Macdonald J, Becker S, Richter L, Körner M. The EUROCROPSML time series benchmark dataset for few-shot crop type classification in Europe. Sci Data. 2025;12(1):664. pmid:40253409
- 47. Potgieter AB, Lawson K, Huete AR. Determining crop acreage estimates for specific winter crops using shape attributes from sequential MODIS imagery. International Journal of Applied Earth Observation and Geoinformation. 2013;23:254–63.
- 48.
Team RC. R: A language and environment for statistical computing. http://www.R-project.org/ 2022.
- 49. Levenberg K. A method for the solution of certain non-linear problems in least squares. Quart Appl Math. 1944;2(2):164–8.
- 50.
Signorell A, Aho K, Alfons A, Anderegg N, Aragon T, Arachchige C. DescTools: Tools for Descriptive Statistics. 2025.
- 51. Belgiu M, Zhou Y, Marshall M, Stein A. Dynamic Time Warping For Crops Mapping. Int Arch Photogramm Remote Sens Spatial Inf Sci. 2020;XLIII-B3-2020:947–51.
- 52. Sardá-Espinosa A. Time-series clustering in R using the dtwclust package. The R Journal. 2019.
- 53. Montero P, Vilar JA. TSclust: AnRPackage for Time Series Clustering. J Stat Soft. 2014;62(1).
- 54.
Bailly S, Giordano S, Landrieu L, Chehata N. Crop-rotation structured classification using multi-source sentinel images and LPIS for crop type mapping. In: IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium, 2018. 1950–3. https://doi.org/10.1109/IGARSS.2018.8518427
- 55. Quinton F, Landrieu L. Crop Rotation Modeling for Deep Learning-Based Parcel Classification from Satellite Time Series. Remote Sensing. 2021;13(22):4599.
- 56. Breiman L. Random forests. Machine learning. 2001;45(1):5–32.
- 57. Wang S, Azzari G, Lobell DB. Crop type mapping without field-level labels: Random forest transfer and unsupervised clustering techniques. Remote Sensing of Environment. 2019;222:303–17.
- 58. Mukherjee M, Khushi M. SMOTE-ENC: A Novel SMOTE-Based Method to Generate Synthetic Data for Nominal and Continuous Features. ASI. 2021;4(1):18.
- 59. Probst P, Wright MN, Boulesteix A. Hyperparameters and tuning strategies for random forest. WIREs Data Min & Knowl. 2019;9(3).
- 60.
Bischl B, Richter J, Bossek J, Horn D, Thomas J, Lang M. mlrMBO: A modular framework for model-based optimization of expensive black-box functions. In: 2017. https://doi.org/10.48550/arXiv.1703.03373
- 61. Wright MN, Ziegler A. ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R. J Stat Soft. 2017;77(1).
- 62. Roberts DR, Bahn V, Ciuti S, Boyce MS, Elith J, Guillera‐Arroita G, et al. Cross‐validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography. 2017;40(8):913–29.
- 63. Maxwell AE, Warner TA, Guillén LA. Accuracy Assessment in Convolutional Neural Network-Based Deep Learning Remote Sensing Studies—Part 1: Literature Review. Remote Sensing. 2021;13(13):2450.
- 64. Blair DC. Information Retrieval. Journal of the American Society for Information Science. 1979;30(6):374–5.
- 65. Laurence L, Heung B, Zhang J, Pennell T, Nyiraneza J, Strom H, et al. Integrating multi-year crop inventories as a proxy for soil management within a digital soil mapping framework for predicting nitrogen indices. Geoderma. 2024;448:116944.
- 66.
Soil Quality—Determination of Dry Bulk Density. ISO. 2017. https://www.iso.org/standard/68255.html
- 67.
Bazan TR. Manual de procedimientos de los análisis de suelos y agua con fines de riego. Instituto Nacional de Innovación Agraria. 2017. https://repositorio.inia.gob.pe/handle/20.500.12955/504
- 68.
Kaklamanos G, Aprea E, Theodoridis G. Mass spectrometry: principles and instrumentation. Chemical Analysis of Food. Elsevier. 2020. p. 525–52. https://doi.org/10.1016/b978-0-12-813266-1.00011-5
- 69. Liu Y, Chen S, Yu Q, Cai Z, Zhou Q, Bellingrath-Kimura SD, et al. Improving digital mapping of soil organic matter in cropland by incorporating crop rotation. Geoderma. 2023;438:116620.
- 70. Vergopolan N, Xiong S, Estes L, Wanders N, Chaney NW, Wood EF, et al. Field-scale soil moisture bridges the spatial-scale gap between drought monitoring and agricultural yields. Hydrol Earth Syst Sci. 2021;25(4):1827–47.
- 71. Guyon I, Eliseeff A. An introduction to variable and feature selection. J Mach Learn Res. 2003;3:1157–81.
- 72. Fisher RA. XV.—The Correlation between Relatives on the Supposition of Mendelian Inheritance. Trans R Soc Edinb. 1919;52(2):399–433.
- 73. Greenwell BM. Pdp: An R package for constructing partial dependence plots. The R Journal. 2017;9(1):421.
- 74. Jonsson P, Eklundh L. Seasonality extraction by function fitting to time-series of satellite sensor data. IEEE Trans Geosci Remote Sensing. 2002;40(8):1824–32.
- 75. Beck PSA, Atzberger C, Høgda KA, Johansen B, Skidmore AK. Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI. Remote Sensing of Environment. 2006;100(3):321–34.
- 76. Kong D, McVicar TR, Xiao M, Zhang Y, Peña‐Arancibia JL, Filippa G, et al. phenofit: An R package for extracting vegetation phenology from time series remote sensing. Methods Ecol Evol. 2022;13(7):1508–27.
- 77. Chandel NS, Tiwari PS, Singh KP, Jat D, Gaikwad BB, Tripathi H, et al. Yield Prediction in Wheat (Triticum aestivum L.) using Spectral Reflectance Indices. Current Science. 2019;116(2):272.
- 78. Chang K, Shen Y, Lo J. Predicting Rice Yield Using Canopy Reflectance Measured at Booting Stage. Agronomy Journal. 2005;97(3):872–8.
- 79. Zhao D, Reddy KR, Kakani VG, Read JJ, Koti S. Canopy reflectance in cotton for growth assessment and lint yield prediction. European Journal of Agronomy. 2007;26(3):335–44.
- 80. De Vroey M, Radoux J, Defourny P. Classifying Sub-Parcel Grassland Management Practices by Optical and Microwave Remote Sensing. Remote Sensing. 2022;15(1):181.
- 81. Morais TG, Domingos T, Falcão J, Camacho M, Marques A, Neves I, et al. Permanent pastures identification in Portugal using remote sensing and multi-level machine learning. Front Remote Sens. 2024;5.
- 82. Vanpoucke K, Heremans S, Buls E, Somers B. Field-Level Classification of Winter Catch Crops Using Sentinel-2 Time Series: Model Comparison and Transferability. Remote Sensing. 2024;16(24):4620.
- 83. Wang R, Feng Q, Jin Z, Ma K, Zhang Z, Liang T. Identification and Area Information Extraction of Oat Pasture Based on GEE—A Case Study in the Shandan Racecourse (China). Remote Sensing. 2022;14(17):4358.
- 84. Ashourloo D, Nematollahi H, Huete A, Aghighi H, Azadbakht M, Shahrabi HS, et al. A new phenology-based method for mapping wheat and barley using time-series of Sentinel-2 images. Remote Sensing of Environment. 2022;280:113206.
- 85. Htitiou A, Boudhar A, Lebrini Y, Hadria R, Lionboui H, Elmansouri L, et al. The Performance of Random Forest Classification Based on Phenological Metrics Derived from Sentinel-2 and Landsat 8 to Map Crop Cover in an Irrigated Semi-arid Region. Remote Sens Earth Syst Sci. 2019;2(4):208–24.
- 86. Delaune T, Lago PD, Hijbeek R, Descheemaeker K, Masolele R, Andersson J. Can Sentinel-2 imagery support the identification of legume fields in maize-dominated smallholder farming systems?. European Journal of Agronomy. 2025;170:127764.
- 87. Ibrahim ES, Rufin P, Nill L, Kamali B, Nendel C, Hostert P. Mapping Crop Types and Cropping Systems in Nigeria with Sentinel-2 Imagery. Remote Sensing. 2021;13(17):3523.
- 88. Poggio L, de Sousa LM, Batjes NH, Heuvelink GBM, Kempen B, Ribeiro E. SoilGrids 2.0: Producing soil information for the globe with quantified spatial uncertainty. SOIL. 2021;7(1):217–40.
- 89. Chen S, Arrouays D, Leatitia Mulder V, Poggio L, Minasny B, Roudier P, et al. Digital mapping of GlobalSoilMap soil properties at a broad scale: A review. Geoderma. 2022;409:115567.
- 90. Prout JM, Shepherd KD, McGrath SP, Kirk GJD, Hassall KL, Haefele SM. Changes in organic carbon to clay ratios in different soils and land uses in England and Wales over time. Sci Rep. 2022;12(1):5162. pmid:35338205
- 91. Six J, Conant RT, Paul EA, Paustian K. Stabilization mechanisms of soil organic matter: Implications for C-saturation of soils. Plant and Soil. 2002;241(2):155–76.
- 92. O’Rourke SM, Angers DA, Holden NM, McBratney AB. Soil organic carbon across scales. Glob Chang Biol. 2015;21(10):3561–74. pmid:25918852
- 93. Varvel GE. Soil Organic Carbon Changes in Diversified Rotations of the Western Corn Belt. Soil Science Soc of Amer J. 2008;72(2):567–8.
- 94. Yang C, Shen F, Li X, Cui W, Zhang L, Yang L, et al. Spatio-temporal mapping reveals changes in soil organic carbon stocks across the contiguous United States since 1955. Commun Earth Environ. 2025;6(1).
- 95. Paustian K, Lehmann J, Ogle S, Reay D, Robertson GP, Smith P. Climate-smart soils. Nature. 2016;532(7597):49–57. pmid:27078564
- 96. Iliquín Trigoso D, Salas López R, Rojas Briceño NB, Silva López JO, Gómez Fernández D, Oliva M, et al. Land Suitability Analysis for Potato Crop in the Jucusbamba and Tincas Microwatersheds (Amazonas, NW Peru): AHP and RS–GIS Approach. Agronomy. 2020;10(12):1898.
- 97. Arce A, de Haan S, Juarez H, Burra DD, Plasencia F, Ccanto R, et al. The Spatial-Temporal Dynamics of Potato Agrobiodiversity in the Highlands of Central Peru: A Case Study of Smallholder Management across Farming Landscapes. Land. 2019;8(11):169.
- 98. Alavi-Murillo G, Diels J, Gilles J, Willems P. Soil organic carbon in Andean high-mountain ecosystems: importance, challenges, and opportunities for carbon sequestration. Reg Environ Change. 2022;22(4).
- 99. Meza K, Vanek S, Olivera E, Huaraca J, Scurrah M, Ccanto R, et al. Perennial forage-based fallows support long-term soil carbon and nutrient trajectories in smallholder rotations of the Peruvian Andes. Agriculture, Ecosystems & Environment. 2025;394:109884.
- 100. Tovar C, Seijmonsbergen AC, Duivenvoorden JF. Monitoring land use and land cover change in mountain regions: An example in the Jalca grasslands of the Peruvian Andes. Landscape and Urban Planning. 2013;112:40–9.
- 101. Olson KR. Soil organic carbon sequestration, storage, retention and loss in U.S. croplands: Issues paper for protocol development. Geoderma. 2013;195–196:201–6.
- 102. Gascon F. In: 2018. 8166–8.
- 103. Zhou X, Wang J, He Y, Shan B. Crop Classification and Representative Crop Rotation Identifying Using Statistical Features of Time-Series Sentinel-1 GRD Data. Remote Sensing. 2022;14(20):5116.
- 104. Gorelick N, Hancher M, Dixon M, Ilyushchenko S, Thau D, Moore R. Google Earth Engine: planetary-scale geospatial analysis for everyone. Remote Sensing of Environment. 2017;202:18–27.
- 105. Vollrath A, Mullissa A, Reiche J. Angular-Based Radiometric Slope Correction for Sentinel-1 on Google Earth Engine. Remote Sensing. 2020;12(11):1867.
- 106.
Eklundh L, Jönsson P. TIMESAT 3.3 with seasonal trend decomposition and parallel processing software manual. 2017.
- 107. Belgiu M, Csillik O. Sentinel-2 cropland mapping using pixel-based and object-based time-weighted dynamic time warping analysis. Remote Sensing of Environment. 2018;204:509–23.
- 108. Burke M, Driscoll A, Lobell DB, Ermon S. Using satellite imagery to understand and promote sustainable development. Science. 2021;371(6535):eabe8628. pmid:33737462
- 109.
Sykas D, Papoutsis I, Zografakis D. Sen4AgriNet: A harmonized multi-country, multi-temporal benchmark misc for agricultural earth observation machine learning applications. In: 2021. 5830–3. https://doi.org/10.1109/IGARSS47720.2021.9553603
- 110. Weikmann G, Paris C, Bruzzone L. TimeSen2Crop: A Million Labeled Samples Dataset of Sentinel 2 Image Time Series for Crop-Type Classification. IEEE J Sel Top Appl Earth Observations Remote Sensing. 2021;14:4699–708.
- 111. Gurgiser W, Juen I, Singer K, Neuburger M, Schauwecker S, Hofer M, et al. Comparing peasants’ perceptions of precipitation change with precipitation records in the tropical Callejón de Huaylas, Peru. Earth Syst Dynam. 2016;7(2):499–515.
- 112. Simoes R, Camara G, Queiroz G, Souza F, Andrade PR, Santos L, et al. Satellite Image Time Series Analysis for Big Earth Observation Data. Remote Sensing. 2021;13(13):2428.
- 113. Bossio DA, Cook-Patton SC, Ellis PW, Fargione J, Sanderman J, Smith P, et al. The role of soil carbon in natural climate solutions. Nat Sustain. 2020;3(5):391–8.
- 114. Reed MS, Graves A, Dandy N, Posthumus H, Hubacek K, Morris J, et al. Who’s in and why? A typology of stakeholder analysis methods for natural resource management. J Environ Manage. 2009;90(5):1933–49. pmid:19231064
- 115. Smith P, Soussana J-F, Angers D, Schipper L, Chenu C, Rasse DP, et al. How to measure, report and verify soil carbon change to realize the potential of soil carbon sequestration for atmospheric greenhouse gas removal. Glob Chang Biol. 2020;26(1):219–41. pmid:31469216
- 116. Batjes NH, Ceschia E, Heuvelink GBM, Demenois J, le Maire G, Cardinael R, et al. Towards a modular, multi-ecosystem monitoring, reporting and verification (MRV) framework for soil organic carbon stock change assessment. Carbon Management. 2024;15(1).
- 117. Körner MF, Leinauer C, Ströher T, Strüker J. Digital Measuring, Reporting, and Verification (dMRV) for Decarbonization. Business & Information Systems Engineering. 2025.
- 118.
World Bank. Technical guidance note on standardizing digital MRV in carbon markets: system evaluation criteria and hotspots assessment. World Bank. 2025. https://openknowledge.worldbank.org/entities/publication/397c4e52-445a-4cf4-89df-f2e61373a524