The study area (Shimbay district) is situated in the northwestern part of Uzbekistan (Figure 1) along the Kegeyli Channel, on the right lower riches of the Amu Darya River in the front part of the Aral Sea delta, covering a 1405 km² area. It is located between 43°20′20.338″N and 59°58′16.041″E (Jumaniyazov et al., 2026). Shiba district is characterized as a flat plain with an absolute elevation of 52–67 m. The climate of the region is sharply continental with a long summer period (+40–45˚C) and a short winter period (−10–12˚C), where the groundwater level is close to the Earth’s surface. The annual precipitation is 100 mm, while evaporation exceeds 2400 mm, demonstrating a typical arid zone. Groundwater is typically shallow and close to the land surface, promoting the capillary rise of saline groundwater and accelerating salt deposition in the root zone.

The district’s agricultural territory is dominated by irrigated croplands fed by the Kegeyli Canal, which are largely used for growing for wheat, rice, maize, and pasture. Alluvial meadows and meadow-alluvial soils have low organic matter content and are highly susceptible to secondary salinization under long-term irrigation (Jumaniyazov et al., 2024). Insufficient natural drainage and an inefficient drainage network contribute to waterlogging and soluble salt accumulation in agricultural areas. Soil salinity is one of the most significant environmental restrictions to agricultural productivity in the region, with moderately and extremely saline soils found across vast expanses of irrigated agriculture. The natural vegetation consists primarily of desert shrubs and sparse halophytic species adapted to saline and arid conditions.

Satellite image processing was conducted using the cloud-based Google Earth Engine (GEE) platform. Sentinel-2 Level-2A surface reflectance imagery from the COPERNICUS/S2_SR_HARMONIZED collection was used to construct a multi-temporal dataset from 2017 to 2025. Images were selected for the main growing season (1 April–31 October) to represent the period of maximum vegetation activity and agricultural development in the Shimbay district.

The image collection was filtered using the CLOUDY_PIXEL_PERCENTAGE metadata filter according to three criteria: (i) the boundary of the study area, (ii) the specified growing season for each year, and (iii) cloud contamination. Images exceeding the predefined cloud threshold were excluded to minimize atmospheric effects. Annual seasonal median composites were then generated from all cloud-filtered Sentinel-2 images acquired during the growing season to reduce residual cloud contamination and provide a consistent spectral representation of the study area.

Following image preprocessing, the Sentinel-2 spectral bands were used to derive the predictor variables employed in the modeling process. These included vegetation, salinity, and soil moisture indices, namely, the NDVI, NDSI, and soil moisture index (SMI). All predictor layers were generated at a spatial resolution of 10 m and exported in the GeoTIFF format using the WGS84 geographic coordinate system (EPSG:4326).

Field observations consisted of 20 sampling locations where soil EC was measured to represent soil salinity. The predictor values extracted from the satellite-derived layers were linked with the corresponding EC measurements to establish the modeling dataset used for statistical analysis.

Three modeling approaches were evaluated for soil salinity prediction: MLR, RF, and PLSR. The model performance was assessed using the coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE). All image preprocessing, predictor extraction, statistical analyses, and model implementation were performed using the cloud-based GEE platform (Google LLC, Mountain View, CA, USA), while QGIS (version 3.34, QGIS Development Team, Open Source Geospatial Foundation, Beaverton, OR, USA) was used for spatial visualization and map preparation. Subsequently, the processed satellite dataset and derived predictor variables were used to analyze the spatial and temporal dynamics of soil salinity and evaluate the predictive performance of the selected modeling approaches across the agricultural lands of the Shimbay district.

To identify the land cover classes of Shimbay district for soil salinity assessment of croplands, the European Space Agency (ESA) WorldCover dataset with 10-meter spatial resolution was obtained. According to (Zanaga et al., 2021), the land cover map from ESA WorldCover has 11 classes. Based on these data, a land use and land cover map of the study area was developed. Due to the natural environmental conditions and the time frame for space images, the study area contained only 8 classes of land cover maps, as shown in Figure 2.

As recognized by (Avdan et al., 2022; Khongnawang et al., 2022) to regulate soil salinity, foundational data are essential, especially regarding the geographic dispersion of the electrical conductivity of a saturated soil paste (ECe − ds/m). According to investigations (Farahmand & Sadeghi, 2020), the EC values of saline soils typically exceed 4 dS/m. Soil samples were obtained from 20 locations (0–30 cm) in November 2025, which were uniformly distributed among agricultural fields covering a 758 km² area (Figure 3). The collected samples were derived from agricultural lands cultivated with corn, rice, and wheat, and EC was measured at the laboratory of the Institute of Agriculture and Agrotechnologies of Karakalpakstan.

Based on the literature review, a previous study on soil salinity modeling assessment (Abdikairov et al., 2024) identified the most appropriate models for soil salinity assessment including Multiple Linear Regression (MLR), RF and Partial Least Square Regression (PLSR). These models were chosen because they reflect typical statistical regression, multivariate regression, and ML methodologies used in soil salinity research based on remote sensing.

MLR was used as a standard statistical method to assess the linear correlations between soil EC and specific spectral characteristics acquired from Sentinel-2 data. The MLR model, which predicts soil salinization as a linear combination of predictor factors, is often employed in ecological and soil research because of its clarity and convenience. MLR is an easy-to-understand method that can be implemented rapidly and efficiently without requiring significant ML techniques (Ngabire et al., 2022). A linear equation is used to predict the relationship between many explanatory variables and a response variable (Shahabi et al., 2017). The relationship between predictors and EC in the models is described by Equation 1, which was developed by (Taghadosi et al., 2019a). where yi is a model-estimated parameter correlated with EC outcomes. The K predictor elements (X1,i,X2,i…,Xk,i) indicates satellite feature indices. The unknown coefficients (a1, a2…,ak) are found during the analysis, while ei represents the error component for every data point.

PLSR was employed to address potential convergence and multidimensionality difficulties using spectral predictors. PLSR reduces the number of initial variables to fewer implicit components while increasing the correlation between prediction factors and responses. This approach is most successful when the spectral bands and indices have a high intercorrelation coefficient and the total number of forecasters is greater than the number of samples. The PLSR model technique applies a simple multidimensional structure to link two separate data matrix variables: the variable of prediction X and the variable of responder Y (Nawar et al., 2014). According to (Sahbeni, 2021), PLSR facilitates the integration of spectral bands and spectral enhancers derived from the same image into a unified model, thus mitigating the collinearity problem prevalent in conventional regression analysis.

RF, a non-parametric group ML approach, was used to record the intricate non-linear interactions between spectral features and soil salinity. RF creates numerous choice trees from bootstrapping examples and randomly selected subsets of predictor parameters. The final projection is formed by combining the results of each tree. The RF model is well-known for its resistance to overfitting, capacity to manage non-linear responses, and superior prediction accuracy in ecological modeling scenarios. RF is a type of supervised ML that is effective for addressing regression and classification issues while identifying non-linear correlations between the desired outcome and input parameters (Zarei et al., 2021). In many cases, it has been applied to evaluate soil salinity, including in the Tarim Basin of China, based on its well-established proficiency in managing excessive dimension and interference in input information (L. Ma et al., 2017).

We used the following variables to develop the models: EC (Figure 3), NDVI, NDSI, SMI (NDMI), and land surface temperature (LST) (Figure 4). The calculation of these variables is described in Table 1.

To evaluate the accuracy of the models, the soil salinity collection was randomly divided into training (70%: 14 soil datasets) and testing (30%: 6 soil datasets). The accuracy of the RF, PLSR, and MLR models was measured using three reliable indicators: coefficient of determination (R²), RMSE, and MAE. The values of the abovementioned metrics range from 0 to +1. R² coefficient of determination represents the statistical connection between the measured and projected values. The RMSE is the square root of the average of the proportional discrepancies among the measured and projected variables. The MAE of the model prediction is calculated from the observed data, regardless of the error direction. These statistics are shown in Equation (2) by (Y. Ma & Tashpolat, 2023), Equation (3) by Fan et al. (2015)and Equation (4) by Hoa(Hoa et al., 2019). where yi signifies the measurements using an average of y̅, y^ refers the projected outcomes, and N represents the total number of observations.

As previously mentioned, 20 soil samples (0–30 cm) were collected from 20 representative locations across the irrigated agricultural lands of the study area, which cover an area of 758 km² area. A stratified spatial sampling technique was used to ensure that the sampling points accurately represented the primary agricultural areas and the existing heterogeneity in land cover and salinity conditions throughout the study area. The sampling sites were spread over the area to maximize spatial coverage while considering field accessibility and logistical restrictions. Table 2 presents soil EC measurements, together with the coordinates of the sampling locations. To facilitate the spatial analysis and visualization of salinity patterns, the measured EC values (0–7.32 dS m⁻¹) were grouped into four salinity classes based on the observed range of field measurements: non-saline (0–2.11 dS m⁻¹), slightly saline (2.11–4.27 dS m⁻¹), moderately saline (4.28–6.11 dS m⁻¹), and highly saline (6.12–7.32 dS m⁻¹). The dominant salts in the study area were sodium and chloride compounds.

This study demonstrated the efficacy of using remote sensing data and modeling tools to estimate soil salinity dynamics in the dry Shimbay district of Karakalpakstan. The multi-temporal examination of Sentinel-2 data indicated significant interannual fluctuations in soil salinity, with NDSI values ranging from 0.011 to 0.058 between 2018 and 2025. The findings revealed a general upward trend in salinization in recent years, particularly in 2024-2025, indicating increased salt accumulation.

Environmental parameters revealed that moisture deficiency and temperature conditions had a significant influence on soil salinity. Consistently negative SMI measurements revealed chronic soil dryness, whereas elevated LST values (up to 32.9°C) increased evaporation and salt concentrations. Despite increased salinity stress, vegetation conditions, as measured by NDVI, remained relatively steady, indicating moderate crop cover. Among the evaluated models, PLSR achieved the highest predictive performance (R² = 0.934, RMSE = 0.398, MAE = 0.329), followed by RF (R² = 0.810, RMSE = 0.966, MAE = 0.728) and MLR (R² = 0.677, RMSE = 0.902, MAE = 0.687). These findings suggest that PLSR is well suited for modeling soil salinity under the study area’s environmental conditions.

Although the proposed approach showed strong predictive performance, the results should be interpreted in light of the study’s limitations, including the limited number of field samples, the reliance on Sentinel-2 optical imagery, the use of a single 70/30 train-test split for model validation, the absence of uncertainty analysis, and the focus on a single study area. Consequently, further validation using larger datasets, more rigorous validation strategies (e.g., k-fold cross-validation), uncertainty assessment, multi-season observations, and different agro-ecological regions is required before broader application.

The findings highlight the potential of combining satellite-derived indices with advanced modeling techniques to support the assessment of regional soil salinity and sustainable land management in arid environments. Future research should incorporate additional environmental variables, expand spatial and temporal coverage, and evaluate model performance under diverse climatic and agricultural conditions.