The study area comprises the 2,324 ha Sardeta ravine basin, located in the municipalities of Montesusín and Robres (Huesca province, Aragón, Spain), within the Los Monegros region, approximately 28 km south of Huesca, with an elevation ranging between 300 and 350 m above sea level. This location within Spain is shown in Figure 2. The area is characterized by an arid climate typical of the Ebro Valley and Lower Aragón, with hot summers, cold winters, and scarce rainfall throughout most of the year. In Montesusín, land consolidation and irrigation modernization have been carried out (mostly from surface irrigation to sprinkler irrigation), while in Robres this has not yet been done, resulting in the predominance of small, terraced plots with surface irrigation.

Crop distribution was characterized using official agricultural statistics from the SINGEAR database (Gobierno de Aragón, 2026), with mean values for 2019–2023. In both municipalities, winter cereals are dominated by barley, which represents 33.23% and 47.17% of the cultivated area in Robres and Montesusín, respectively, compared with 18.01% and 11.47% for wheat. In irrigated areas of Montesusín, alfalfa (25.61% of the total area) and maize (17.58%) are also widely cultivated. Maize is frequently grown as a second crop following barley, representing 48.34% of the barley area, whereas double cropping after wheat is considerably less common (16.02% of the wheat area) and in some years accounts for less than 3% of the cultivated surface. This difference is mainly associated with the longer growth cycle of wheat, which delays the establishment of the second crop and reduces its agronomic viability. These patterns indicate that barley and maize are the predominant cropping systems in the study area and were therefore selected as representative crops for the analysis.

The soils of the basin developed on limestone sediments and are classified as Xerorthent (poorly developed) and Calcixerept (moderately developed). A soil map of the basin is shown in Figure 2, identifying seven distinct homogeneous edaphic units (Usón Murillo et al., 2023). Table 1 describes these units. Due to the complexity and low representativeness of units 2 and 3 (Shallow Soils and Saline Association), these units were excluded from the variability analysis. Additionally, units 4 and 6 (Valley Bottom and Floodplain) were grouped into a single unit belonging to the same soil subgroup.

Agronomic practices, including sowing dates, irrigation management, and fertilization, were uniform within each plot; therefore, the observed spatial variability in crop development can be primarily attributed to differences in soil properties.

The applied methodological framework was designed to characterize soil-driven crop variability under irrigated semi-arid conditions through a progressive analysis from plot-scale assessment to controlled within-plot evaluation.

Among available remote sensing approaches, the Normalized Difference Vegetation Index (NDVI) was selected due to its simplicity, robustness, and wide availability from optical satellite imagery. NDVI provides a consistent and operational indicator of crop development and spatial variability while requiring minimal data processing, making it suitable for cost-effective and transferable agricultural applications. It has been widely used as an indirect indicator of soil-driven crop variability, showing strong relationships with crop productivity, soil physicochemical properties, and spatial yield patterns (Amezketa Lizarraga et al., 2011; Gopp et al., 2017; Gopp and Savenkov, 2019; Lobell et al., 2003). Although NDVI may exhibit saturation under dense canopy conditions, particularly during peak crop development, this limitation does not compromise its suitability for the objectives of this study. The analysis focuses on spatial variability and temporal dynamics rather than absolute biomass estimation and includes phenological stages where saturation effects are less critical. Moreover, the use of dispersion metrics (standard deviation) enhances the detection of within-plot heterogeneity even under high canopy cover conditions.

Given the objective of developing a practical and transferable methodology for within-plot variability assessment, NDVI derived from optical satellite imagery was selected as the most efficient and cost-effective indicator. Its wide availability, minimal processing requirements, and low implementation cost make it particularly suitable for research and operational applications under limited resource conditions.

NDVI time series were derived from Sentinel-2 satellite Level-2A surface reflectance imagery using Sentinel Hub EO Browser (Sentinel Hub, 2026). Sentinel-2 provides multispectral observations with a spatial resolution of 10 × 10 m per pixel, enabling the characterization of crop development at the plot scale. For each study plot, polygon boundaries previously defined in QGIS (version 3.18) were uploaded to EO Browser to extract NDVI values for all available satellite acquisitions during the 2019–2023 period. NDVI was computed from the red and near-infrared spectral bands and exported as pixel-level data in comma-separated values (CSV) format.

The analysis was conducted using agricultural plots for which detailed soil characterization was available, allowing the direct linkage between crop development and soil properties. A total of 12 plots were initially selected based on soil profile information; in four plots complete soil profile pits were excavated and sampled by horizons, while the remaining plots were characterized through auger sampling. For the satellite-based analysis, only plots located within the most representative soil units of the study basin were considered. Two edaphic units characterized by low spatial representativeness were excluded, resulting in a final dataset of 11 plots distributed across Platforms (three plots), Glacis–Slope (four plots), Redoxmorphic Floodplain (one plot), and Valley Bottom–Floodplain (three plots).

To ensure the reliability of the NDVI time series, cloud-contaminated observations were removed using the cloud detection algorithm implemented in EO Browser, which is based on the Braaten–Cohen–Yang method (Sentinel Hub, 2026). Only satellite acquisitions with 0% cloud coverage within each plot were retained to minimize potential contamination of NDVI values. Although this strict selection may reduce the temporal density of the image series and potentially lead to some degree of under-sampling in the NDVI time series, the remaining observations were sufficient to identify the main phenological stages analyzed in this study.

Pixel-level NDVI data were imported into QGIS (version 3.18), where zonal statistics were computed for each plot and acquisition date. Mean, median, standard deviation, and the 10th and 90th percentiles were calculated from all pixels within each plot, transforming spatially distributed information into parcel-level metrics describing crop development and within-plot variability. These statistics were associated with the acquisition date of each satellite image, enabling the construction of parcel-level NDVI time series for analyzing crop temporal dynamics. NDVI dispersion metrics, primarily represented by the standard deviation, were used to characterize within-plot spatial heterogeneity. These indicators enhance the detection of soil-driven variability and remain sensitive to spatial differences in crop performance even when NDVI values approach saturation.

The resulting NDVI time series were analyzed to characterize crop dynamics and identify crop cycles throughout the study period. The temporal evolution of NDVI values enabled the differentiation between perennial and annual crops based on their distinct growth patterns. Perennial crops such as alfalfa exhibited multiple growth cycles within a year, resulting in recurrent NDVI peaks that hindered the consistent identification of phenological stages (Figure 3A). In contrast, annual cereal crops showed a single, well-defined growth cycle, typically following the sequence barley–second maize crop (Figure 3B). Consequently, the analysis focused on winter cereals (mainly barley) and summer crops (mainly maize), while perennial crops were excluded from further analysis.

Crop phenological stages were determined from NDVI time series using a threshold-based approach combined with a temporal consistency criterion to identify the start of the season (SOS), maximum development (MAX), and end of the season (EOS). Although previous studies define phenological transitions based on relative thresholds of seasonal NDVI amplitude (Ren et al., 2017; Yin et al., 2025), these approaches do not fully account for signal variability associated with intensive irrigation management, heterogeneous crop emergence, or transient disturbances commonly observed in semi-arid irrigated systems. Therefore, threshold values and detection criteria were adapted to local agronomic conditions to ensure robust identification of crop development stages.

In the study area, crops are intensively managed under irrigation, typically reaching high NDVI values and exhibiting rapid senescence following the cessation of irrigation prior to harvest. Moreover, NDVI temporal series may present short-term fluctuations caused by soil background effects, weed growth during early stages, spatial heterogeneity within plots, or temporary irrigation interruptions. These factors may produce transient NDVI increases or decreases that do not represent actual phenological transitions.

Early vegetation signals detected before crop establishment may represent weed growth or residual vegetation rather than the primary crop. This can confound the identification of phenological transitions when using vegetation index time series (Sarvia et al., 2026). However, crop establishment in the study area generally occurs rapidly due to warm conditions and intensive agronomic management, particularly under double-cropping systems, and the widespread use of pre-emergence herbicides typically limits early weed infestation.

To avoid misclassification of phenological transitions, a temporal consistency criterion based on consecutive observations was defined:

The SOS stage was determined as the fourth consecutive observation following three successive NDVI increases with values exceeding 0.4 (approximately 40% of seasonal amplitude). This threshold was selected to ensure sustained crop establishment while minimizing the influence of early signal variability. The criterion was particularly relevant for maize, where typical row spacing (∼0.7 m) results in substantial soil exposure and irregular NDVI values during initial growth stages.

The MAX stage corresponded to the maximum NDVI value recorded during the crop cycle, representing peak canopy development. Although NDVI may saturate under dense canopy conditions, this stage provides a consistent reference for crop maturity under irrigated management.

The EOS stage was defined as the fourth consecutive observation following three successive NDVI decreases after the maximum value, when NDVI values fall below 0.9 (approximately 90% of seasonal amplitude), representing the onset of crop senescence. This criterion captures the initial decline in vegetation activity rather than advanced senescence, allowing the detection of spatial differences in crop drying related to soil water availability and water-holding capacity. The use of consecutive observations further ensured that detected transitions reflected sustained phenological changes rather than temporary NDVI fluctuations.

Although this approach enables comparison among soil units at the plot scale, differences in agronomic management between farmers may partially obscure soil-driven crop responses. This limitation is addressed through a controlled within-plot analysis presented in the following section.

Statistical analyses were performed using IBM SPSS Statistics (version 26). Analysis of variance was used to evaluate differences in NDVI-derived metrics among soil units, and mean comparisons were conducted using the Tukey B test when significant differences were detected (p < 0.05).

To minimize the influence of agronomic management and isolate soil-driven crop responses, we conducted a controlled within-plot analysis. This approach allows the isolation of soil effects under homogeneous management conditions.

A large agricultural plot (21 ha) located in the municipality of Montesusín (41.887777°N, 0.386029°W) was selected because it includes two representative soil units (Platform and Glacis–Slope) with contrasting edaphic characteristics. The field is managed uniformly with respect to sowing date, fertilization, irrigation, phytosanitary treatments, and harvesting operations. Under these homogeneous management conditions, differences in crop development can be primarily attributed to soil variability.

To characterize crop performance within each soil unit while avoiding non-representative areas such as field edges, headlands, or traffic zones, six polygons of approximately 1,000 m² were randomly delineated within each unit, resulting in twelve sampling areas. This random selection ensured that the sampled zones represented dominant soil conditions while minimizing localized disturbances. Figure 4 illustrates the spatial distribution of the sampling polygons within the study plot, with striped areas corresponding to the Glacis–Slope unit and dotted areas to the Platform unit.

Crop cycles within the selected plot were first identified using Sentinel-2 NDVI time series derived from EO Browser, following the procedure described in Section NDVI-Based Assessment of Plot-Scale Variability Among Soil Units. Phenological stages (SOS, MAX, and EOS) were determined from these time series. Subsequently, high-resolution PlanetScope imagery (3 m spatial resolution) corresponding to the identified dates was acquired, allowing these stages to be analyzed at a finer spatial resolution. This two-step approach enabled precise temporal identification of crop development while improving the spatial characterization of within-plot variability.

NDVI values were extracted for all crop cycles recorded between 2019 and 2023, including winter crops (barley) and summer crops (maize). The use of a consistent phenological framework enabled comparison of crop responses between soil units under identical management conditions across multiple growing seasons.

Statistical analyses were performed using IBM SPSS Statistics (version 26). Data distribution was evaluated prior to analysis, and a two-factor univariate analysis of variance was conducted to assess the effects of soil unit and crop type on NDVI median values and NDVI standard deviation at the three phenological stages.

To further interpret the soil-driven variability observed under homogeneous management conditions, an analysis was conducted to identify the soil properties most strongly associated with crop performance, using NDVI-derived parameters as indicators of crop response. This approach supports the interpretation of spatial variability in terms of soil–crop interactions.

The analysis incorporated variables grouped into three categories: soil properties (soil unit, clay and sand content, water table depth, depth of redoximorphic features, surface horizon organic matter, and salinity indicators); plot characteristics (plot area, mean elevation, slope, and irrigation system); and crop-related variables (crop type, growing season, and NDVI-derived indicators).

Spearman’s rank correlation analysis was used to evaluate the relationships among the variables. Although most variables followed a normal distribution, Spearman’s coefficient was selected to assess monotonic relationships, which are expected in soil–crop systems where responses may not follow strictly linear patterns. This approach provides a robust assessment of directional associations while accommodating potential non-linear trends. Correlation coefficients (ρ) were considered significant at the 95% confidence level.

In addition, principal component analysis (PCA) was performed to explore multivariate relationships and identify the main factors explaining variability within the soil–crop system. This approach facilitates the interpretation of complex interactions by reducing data dimensionality and highlighting the variables that contribute most strongly to soil–crop system variability.

To evaluate crop response across soil units at the plot scale, NDVI-derived metrics were analyzed at key phenological stages (SOS, MAX, and EOS). In addition, the seasonal MEAN NDVI was considered to provide an integrated indicator of crop performance throughout the growing cycle. NDVI median values were used as an indicator of overall crop development, whereas the NDVI standard deviation was used to characterize within-plot spatial heterogeneity. A two-factor analysis of variance was conducted to assess the effects of soil unit and crop type on NDVI metrics across growth stages.

NDVI values did not differ significantly among soil units or between crop types across the evaluated growth stages. These results indicate that none of the soil units imposed significant limitations on crop development and that growth patterns were comparable between summer and winter crops. Significant differences were observed only among phenological stages, with the lowest values at emergence, a peak at maximum development, and a decline during senescence (Figure 5).

The absence of significant differences in NDVI median values (0.6577 for summer crops and 0.6502 for winter crops) is consistent with NDVI saturation under dense canopy conditions (Tucker, 1979). This occurs because red reflectance saturates earlier than near-infrared, reducing sensitivity at high biomass levels (Salvador-Castillo et al., 2021), a behavior linked to the strong absorption of chlorophyll in the red spectral region (Gitelson et al., 2003). In this context, the limited discrimination at maximum canopy development reinforces the importance of focusing on SOS and EOS stages, as well as on dispersion metrics, to better capture soil-driven variability.

In contrast, NDVI standard deviation differed significantly among soil units (Figure 6). The Glacis Slope unit exhibited the highest variability, whereas the Platform unit showed the lowest values, with Floodplain and Redoximorphic Floodplain displaying intermediate behavior (letter groupings in Figure 6). A tendency toward higher variability in winter crops compared with summer crops was also observed (p = 0.071), suggesting greater spatial heterogeneity under less uniform water supply conditions. This pattern may reflect more homogeneous irrigation management in maize, whereas barley shows greater dependence on rainfall, leading to increased spatial variability in crop response; alternatively, it could also be influenced by earlier or stronger NDVI saturation under the denser maize canopy, potentially reducing the sensitivity of the index to spatial differences.

The crop developmental stage also influenced spatial variability. NDVI standard deviation increased from SOS (0.043) to EOS (0.061), with the highest variability observed during senescence. This suggests that crop senescence is the most suitable stage for detecting soil-driven differences in crop response. In contrast, variability at SOS was more strongly affected by differences in sowing dates and early establishment conditions.

Variability in crop emergence timing among plots was assessed by comparing the SOS dates across all evaluated crops. Substantial differences were observed, with emergence varying by up to 3 months for winter crops and by approximately 2 months for summer crops. These differences reflect variations in sowing dates and management practices among plots and may introduce uncertainty when interpreting soil-driven crop responses. To minimize these effects, subsequent analyses were conducted under homogeneous management conditions.

To isolate the influence of soil properties from management effects, a plot with homogeneous agronomic management was analyzed in detail. Uniform irrigation, fertilization, and sowing practices across the field enabled evaluation of soil-driven variability under consistent conditions.

The cropping sequence within the selected plot over the last 5 years was analyzed to provide temporal context for the study. Figure 7 shows the growth cycles of representative winter and summer crops, including barley (2021 and 2023) and maize (2019–2023), together with the timing of the SOS, MAX, and EOS stages. Main-crop maize and double-crop maize following barley are distinguished, with the latter exhibiting a delayed growth cycle due to later planting. This temporal framework confirms the consistency of phenological development across seasons and supports the use of these stages to evaluate soil-driven variability.

Crop response is discussed across the key phenological stages.

To better understand the mechanisms underlying the observed variability, relationships among soil properties and crop response were examined using correlation analysis and principal component analysis (PCA).

Correlation analysis using Spearman’s rank coefficients was conducted to explore relationships among soil properties, plot characteristics, and crop response indicators derived from NDVI metrics. Only statistically significant relationships are presented to facilitate interpretation of the main drivers of crop variability, as summarized in Figure 11. Relationships are described below according to three variable groups.

Soil-related variables showed strong interconnections reflecting their influence on soil hydrological behavior and crop development. Soil unit was associated with several edaphic and groundwater-related attributes, confirming its role as an integrative indicator of soil conditions across the watershed. Texture-related properties exhibited contrasting patterns consistent with their control over soil water dynamics: sand content was associated with conditions favoring permeability and drainage, whereas clay content reflected greater water retention and a tendency toward salt accumulation. Organic matter content emerged as a stabilizing factor, being associated with improved soil structure and more homogeneous crop development. Relationships involving salinity and redoximorphic features further indicate the influence of groundwater conditions and soil moisture regimes.

Topographic variables related to plot morphology were closely linked to soil physical properties and hydrological behavior. Slope was associated with increased crop variability, suggesting uneven soil moisture redistribution along the terrain. Elevation showed relationships consistent with differences in drainage conditions and soil development across the landscape. The irrigation system was associated with variability patterns, reflecting differences in water application uniformity and management intensity.

Crop-related indicators derived from NDVI confirmed the sensitivity of remote sensing metrics to soil-driven heterogeneity. NDVI standard deviation, used as an indicator of within-plot crop variability, showed strong monotonic associations with terrain and soil moisture–related variables. In contrast, NDVI mean values were inversely related to variability, indicating that greater crop vigor is associated with more homogeneous canopy development.

Collectively, the correlation patterns indicate that spatial variability in crop development is primarily controlled by soil-water dynamics. The interaction of texture, topography, groundwater influence, salinity, and organic matter governs soil moisture distribution, which in turn determines crop response. These relationships support the use of NDVI standard deviation as a robust indicator of soil-driven heterogeneity in agricultural landscapes.

Principal Component Analysis (PCA) was performed to visualize how variables group together and to provide a clearer interpretation of the relationships identified in the correlation analysis. This multivariate approach complements the correlation results by revealing underlying gradients that control crop variability.

The first two principal components explained 37.33% of the total variance, reflecting the complexity of the system and the diversity of variables considered, including soil properties, terrain attributes, and crop indicators collected across different crops and years.

Component 1 was primarily associated with textural fractions and elevation, with clay and sand showing opposite contributions. This axis represents a gradient of soil texture and permeability that influences water movement and drainage conditions.

Component 2 was strongly related to soil unit, slope, and organic matter content, representing a gradient of soil heterogeneity and moisture retention capacity across the landscape.

The loading plot (Figure 12) reveals meaningful groupings of variables that support the interpretation of soil hydrological processes. Variables located close to each other are positively correlated, whereas those positioned on opposite sides of the origin are negatively correlated. Sand content clustered with water table depth, reflecting high permeability and enhanced vertical water movement. The irrigation system was grouped with salinity, indicating that areas affected by soil salinity tend to maintain traditional surface irrigation practices. Elevation was clustered with redoximorphic feature depth and organic matter content, consistent with the coarser texture and improved drainage conditions of the Platform unit. The soil unit was grouped with NDVI standard deviation, reinforcing the role of this indicator as a proxy for soil-driven crop variability. To facilitate interpretation, the loadings of the variables on the first two principal components are presented in Table 2.

Overall, the PCA highlights the dominant influence of soil texture, terrain, and groundwater on soil moisture dynamics, revealing patterns consistent with the associations identified in the correlation analysis. These processes govern crop spatial variability, particularly in finer-textured areas where salinity accumulation and redoximorphic conditions are more likely to occur.

Taken together, these results lead to the following conclusions.

NDVI-derived metrics proved effective for assessing crop development and spatial variability at the plot scale. NDVI mean values reflected overall crop performance, whereas NDVI standard deviation captured within-plot variability associated with soil heterogeneity.

The use of key phenological stages (SOS, MAX, and EOS) simplified the analysis and facilitated interpretation of crop response. Variability in agronomic management partially masked soil-related differences at the watershed scale, but these became clearer under homogeneous management conditions.

NDVI standard deviation revealed differences in crop development between soil units, with more homogeneous growth in the Platform unit and greater heterogeneity in the Glacis Slope unit. Crop senescence (EOS) showed the greatest spatial variability, making it the most suitable stage for detecting soil-driven differences, while emergence timing (SOS) also showed potential for supporting precision management of the current crop.

Taken together, these findings demonstrate that NDVI variability metrics provide a robust approach for identifying soil-driven crop heterogeneity and supporting site-specific management.

The main limitations relate to differences in agronomic management and to the inherent constraints of NDVI under dense canopy conditions. Future research should evaluate additional vegetation indices and sensing approaches, including the integration of synthetic aperture radar and optical observations, to provide complementary information on soil moisture and structural properties linked to spatial variability. This methodology should also be tested across diverse agroecosystems to further support precision agriculture and sustainable land management.