Remote Sensing-Based Assessment of Hayfield Productivity and Forage Quality in Mongolia’s Steppe Ecosystems

In Bornuur sub-province, agricultural intensification and livestock pressure have driven land degradation and reduced hay productivity to only 565.7 tons from 7,443.2 ha in 2020 (Callaghan et al., 2024; Danzhalova et al., 2023; Jamsranjav et al., 2018; NSOM, 2025).

Fieldwork was conducted on 15–16 August 2020, in Bornuur sub-province, Tuv Province, Mongolia. It focused on the major forage production areas. Twenty sampling sites were selected using a stratified random design to capture variation in land use intensity and topography.

At each site, aboveground biomass was collected from 1 m² quadrats, clipped at ground level, oven-dried at 65°C, and weighed to determine the dry matter yield (g/m²), in line with national monitoring protocols (NSOM, 2025).

Plant species composition was recorded using the Hayfield Plant Identifier manual (Zhang et al., 2021) to identify dominant and subdominant species. Subsamples of dominant species were analysed at the Mongolian University of Life Sciences using standard methods to determine crude protein, NDF, ADF, and in vitro dry matter digestibility (Ren et al., 2016).

This approach enabled integrated assessment of both forage yield and nutritive quality under existing land use conditions.

Topographic analysis was based on the ALOS PALSAR DEM (12.5 m resolution) from JAXA, pre-processed for GIS use (Rosenqvist et al., 2007). Slope, aspect, and elevation were derived as key terrain variables. Slope influences erosion, and vegetation patterns (Hou et al., 2014), aspect controls microclimate conditions such as solar radiation and soil moisture (Liu et al., 2012), and elevation affects temperature and precipitation gradients. These variables were used to assess topographic effects on forage yield and quality.

Sentinel-2 (10 m) imagery acquired on 21 July 2020 (ID: S2B_MSIL2A_20200721T035539_N0500_R004_T48UWU) – 25 days prior to the field campaign – was obtained from Google Earth Engine as the only cloud free scene within a 30day window. To minimise phenological discrepancies, we assumed relatively stable vegetation conditions between late July and midAugust in this steppe ecosystem (Drusch et al., 2012). The July 21 image was selected as the closest cloud-free acquisition to the August 15–16 field campaign date. Relevant spectral bands, including red (B4), green (B3), blue (B2), near-infrared (B8), and shortwave infrared (B11), were used for analysis and vegetation index calculation. Sentinel-2 Level-2A images are atmospherically corrected using the Sen2Cor algorithm. Pixels with cloud cover >10% were excluded using the QA60 band.

Among the indices, the NDVI was calculated using the red and near-infrared bands to assess vegetation condition. It ranges from 1 to 1, with higher values indicating greater vegetation density and vigour (Rouse et al., 1974). (1) NDVI=NIR−REDNIR+RED {\rm{NDVI}} = {{\left( {{\rm{NIR}} - {\rm{RED}}} \right)} \over {\left( {{\rm{NIR}} + {\rm{RED}}} \right)}} where:

• RED – red band, NIR – near-infrared band

The Enhanced Vegetation Index (EVI) improves upon NDVI by reducing atmospheric effects, canopy background noise, and soil influences. It is especially effective in areas with dense vegetation and uses additional correction factors to enhance sensitivity to canopy structure and background conditions. (2) EVI=G×NIR−REDNIR+C1×RED−C2×BLUE+L {\rm{EVI}} = G \times {{\left( {{\rm{NIR}} - {\rm{RED}}} \right)} \over {\left( {{\rm{NIR}} + {C_1} \times {\rm{RED}} - {C_2} \times {\rm{BLUE}} + L} \right)}} where:

• G – gain factor, C1 and C2 – coefficients for the red and blue bands, L – canopy background adjustment factor (Huete, 1988).

The Soil-Adjusted Vegetation Index (SAVI) is designed for areas with sparse vegetation and reduces the influence of soil brightness on spectral reflectance. It modifies NDVI by incorporating a soil adjustment factor (L), improving accuracy in low vegetation cover conditions. (3) SAVI=NIR−REDNIR+RED+L×1+L {\rm{SAVI}} = {{\left( {{\rm{NIR}} - {\rm{RED}}} \right)} \over {\left( {{\rm{NIR}} + {\rm{RED}} + L} \right)}} \times \left( {1 + L} \right)

• L – is typically set to 0.5 in regions with sparse vegetation (Huete, 1988).

The Normalised Difference Water Index (NDWI) is used to assess vegetation water content and detect surface water bodies. It is calculated using near-infrared and shortwave infrared bands, where higher values indicate greater moisture content or the presence of water (Gao, 1996). (4) NDWI=(NIR−SWIR)(NIR+SWIR) {\rm{NDWI}} = {{({\rm{NIR}} - {\rm{SWIR}})} \over {({\rm{NIR}} + {\rm{SWIR}})}} where:

• NIR – near-infrared band, SWIR – Shortwave-infrared band

Pearson’s correlation analysis was applied to quantify the relationships between vegetation indices, topographic variables, and field-measured forage yield and nutritive quality (Huete, 1988; Rouse et al., 1974). Linear regression models were developed to assess the predictive performance of vegetation indices for biomass yield (Booth and Tueller, 2003). Descriptive statistical measures, including mean, standard deviation, skewness, and kurtosis were computed to characterise spatial variability in the dataset (Tabachnick, 2013).

Multivariate regression analysis was conducted to evaluate the combined predictive power of all variables. LASSO (Least Absolute Shrinkage and Selection Operator) regularisation was applied for variable selection using 5-fold cross-validation to prevent overfitting, given the limited sample size (n = 20) (Hastie, 2015; Tibshirani, 1996). Model performance was assessed using R², Root Mean Square Error (RMSE), and Mean Absolute Error (MAE) (James et al., 2013). SHAP (SHapley Additive exPlanations) analysis was performed to quantify the marginal contribution of each variable to yield predictions (Lundberg and Lee, 2017). All statistical analyses were performed using RStudio (version 2024.04.0+735) with the glmnet, caret, and shapviz packages (Friedman et al., 2010; R Core Team, 2024).

During the final ten days of August 2020, vegetation surveys were conducted at 20 designated sampling points within the hayfields of Bornuur sub-province (Fig. 3).

Fig. 3.

Field work points location

Source: own elaboration.

The lower ecotones of mixed larch – pine – birch – poplar forests in Bornuur sub-province support diverse hayfield vegetation types. These include various combinations of forb, grass, wormwood, and legume communities (M1–M11, M13–M20), reflecting high compositional heterogeneity across sites. In contrast, hayfields in the Boroo River valley are mainly dominated by grass species (M12), indicating a simpler vegetation structure (Table 1).

Following Gao et al. (2015), hayfields were classified into high- (>12 c/ha), medium- (9–11 c/ha), and low-yield (<8 c/ha) systems, with community composition varying accordingly (Table 1).

Topographic parameters (elevation, slope, and aspect) for the Bornuur sub-province were derived from the ALOS PALSAR DEM (12.5 m resolution) (Fig. 4), and descriptive statistics were calculated from 20 field sampling points (Table 2).

Fig. 4.

Land surface factors map of Bornuur sub-province. a – elevation, b – slope, c – aspect

Source: own elaboration.

Topographic parameters derived from ALOS PALSAR DEM are summarised in Table 2.

Aspect showed a moderate negative correlation with yield (r = −0.44, p = 0.051), approaching but not reaching statistical significance at the α = 0.05 level. This suggests a potential microclimatic influence, where southfacing slopes may experience greater water stress due to higher solar radiation (Liu et al., 2012).

Vegetation indices were derived from cloud-free Sentinel-2 imagery (August 2020) for Bornuur sub-province, and descriptive statistics were calculated for 20 field sites (Fig. 5).

Fig. 5.

Vegetation indices map of Bornuur sub-province. a – EVI, b – NDWI, c – NDVI, d – SAVI

Source: own elaboration.

Notably, the NDWI showed very low spatial variability across the 20 sampling sites (SD = 0.027, coefficient of variation = 5.0%), indicating uniformly dry surface conditions during the satellite overpass. This low variance limits the discriminatory power of the NDWI for explaining yield differences among sites.

Vegetation indices derived from cloud-free Sentinel-2 imagery (August 2020) were analysed across 20 field sites in Bornuur sub-province (Fig. 6).

Fig. 6.

Correlation matrix of variables *p < 0.05, **p < 0.01, ***p < 0.001. Values without asterisks are not statistically significant

Source: own elaboration.

Linear regression analysis assessing the relationship between forage yield and vegetation indices, along with surface environmental variables, demonstrated strong predictive performance for the NDVI and EVI (Fig. 7–8). The coefficient of determination (R²) was 0.78 for NDVI and 0.58 for EVI, indicating substantial explanatory power in modelling yield variability. These results highlight the effectiveness of satellite-derived NDVI and EVI as reliable proxies for estimating forage crop productivity, supporting their application in precision agriculture and data-driven land management strategies.

Fig. 7.

Linear regression plots of forage yield versus vegetation indices (NDVI, EVI, SAVI, NDWI)

Source: own elaboration.

Fig. 8.

Linear regression analysis: Yield and land surface factors

Source: own elaboration.

To address the combined effect of vegetation indices and terrain variables on predicting forage yield (Research Question 3), we compared two models: a FULL model incorporating all predictors (NDVI, EVI, SAVI, NDWI, elevation, slope, aspect, total protein, general nutrition, total ash, harvest residue) and a LASSO-regularised model for variable selection. Model performance was evaluated using 5-fold cross-validation.

The FULL model achieved an R² of 0.879, RMSE of 0.853, and MAE of 0.653 (Table 4). The LASSO model showed comparable performance (R² = 0.869, RMSE = 0.889, MAE = 0.677), indicating that the full set of predictors does not substantially overfit despite the limited sample size (n = 20).

Variable importance analysis (Table 5) revealed that the NDVI had the strongest positive effect on yield prediction (estimate = 18.83, p = 0.022), confirming its dominant role. The EVI and NDWI showed moderate positive effects but were not statistically significant at α = 0.05. Slope exhibited a negative but non-significant effect (estimate = −0.289, p = 0.266), while terrain variables (elevation, aspect) and forage quality metrics (protein, general nutrition) showed negligible contributions.

SHAP (SHapley Additive exPlanations) analysis further confirmed that the NDVI made the greatest absolute contribution to yield predictions (mean SHAP = 1.567), followed by EVI (0.698) and slope (0.680) (Table 6). All other variables contributed minimally (< 0.25).