The current work was conducted in the Béni Amir irrigated sub-perimeter of the Tadla Plain, covering 674 km² in central Morocco (Figure 1). The Tadla Plain is a major hydro-agricultural region of Central Morocco. Formerly barren and exploited for pastoral purposes, this region has become a fertile agricultural area following the installation of an irrigation network, and nowadays contributing to a large proportion of national agricultural production, up to 30% for sugar beets, 12% for fodder, 11% for citrus fruits and olives, 10% for market gardening, 6% for cereals and 10% for milk (11). Since the development and commissioning of the irrigated perimeter, the salt-affected area of agricultural lands has been steadily increasing (10). Irrigation water sources include both shallow groundwater and surface water from the Oum Er Rbia River. Reported salinities are on the order of about 3.2 g/L for groundwater and about 1.3 g/L for the Oum Er Rbia surface water (8). In practice, some farmers also blend surface water with groundwater at the field or parcel scale (24), which can further vary the salinity of applied irrigation water across the perimeter. The irrigation network in the region, managed by the Regional Office for Agricultural Development of Tadla (ORMVAT) and supplied by the Oum Er Rbia River, primarily relies on surface (gravity) irrigation, which leads to considerable water losses due to inefficient infrastructure and evaporation. Although some farmers use sprinkler systems, a growing number are shifting toward drip irrigation to improve water-use efficiency (25).

The geology of the region is characterized by a vast syncline filled with sedimentary deposits accumulated during the Cretaceous and Tertiary eras (26).

The pedology reveals considerable soil heterogeneity in the Tadla Plain. Kastanozems, which cover 43% of the land, are rich in organic matter and support soil fertility and agricultural productivity. On the other hand, Leptosols, accounting for 32% of the study area, contain high levels of calcium and magnesium, which influence soil chemistry and plant growth. Nitisols-Alisols, cover 18% of the area. Other various soil types characterize the remaining 7% of the area (27).

The climate of the study area is arid to semi-arid, with annual precipitation ranging from 150 to 450 mm (28). The dry season, which extends from April to October, is characterized by minimal rainfall, typically between 0 and 50 mm per month. In contrast, the rainy season, which occurs from November to March, accounts for approximately 70% of the total (29). Temperatures exhibit significant seasonal variations, with a maximum of 46°C in August and a minimum of -6°C in January, resulting in an annual average of 20°C. The yearly average evaporation is approximately 1800 mm, nearly six times the annual cumulative rainfall (30). The average altitude ranges from 350 m to 500 m, and the overall slope is less than 6°, with the lowest point located at Sidi-Driss (31).

The methodology used in this study is summarized in Figure 2. The ground truth data were collected by sampling soil from 43 georeferenced locations between October 28 and 31, 2021, using a stratified random sampling approach to ensure a representative assessment of soil variability across the Béni Amir sub-perimeter. The sampling design was informed by previous studies in the Tadla plain, which identified significant gradients. Each sampling point corresponded to a homogeneous 30 m × 30 m area, matching the spatial resolution of the Landsat-8 OLI reflective bands used in this study. Field conditions during sampling were largely post-harvest, with minimal vegetation cover, ensuring that satellite reflectance captured soil rather than canopy characteristics. The electrical conductivity (ECe) of each soil sample (0–10 cm) was analyzed in the laboratory using the saturated paste extract method, as described by Rhoades (32). The sampled soils were classified into five salinity classes (Supplementary Table S1 in the Supplementary Material) in accordance with Ivushkin et al. (33). To support model training and testing, given the limited dataset, the samples were resampled to 144 sample instances using a controlled data augmentation strategy detailed in Section 2.3.

Concurrently, Landsat imagery acquired on November 12, 2021(cloud cover = 0.11%), downloaded from (https://earthexplorer.usgs.gov/, accessed on December, 1st 2021), and preprocessed in QGIS 2.18.0. The image has been radiometrically calibrated and atmospherically corrected using the Dark Object Subtraction (DOS) (23, 34). The data was geometrically corrected to align with the collected ground truth data, facilitating seamless integration for comparative analysis.

Furthermore, spectral bands spanning from the visible, including Blue (B2), Green (B3), and Red (B4), to the short-wave infrared wavelengths, SWIR1 (B6) and SWIR2 (B7), are used as recommended by previous studies (35, 36). Their efficiency in identifying and mapping salt-affected soils underscores their essential role in assessing soil degradation in both agricultural and natural landscapes. Spectral indices were calculated, including soil salinity indices (37, 38), intensity indices (39), brightness indices (38), and the land degradation index (LDI) (7, 40) (Supplementary Table S2, Supplementary Figure S1 in the Supplementary Material).

Data preprocessing was performed in a Visual Studio Code environment using the Python programming language, involving the application of resampling techniques to overcome class imbalances, scaling methods to normalize spectral and laboratory data, and principal component analysis (PCA) to address multicollinearity among predictor variables, and to enhance the stability and performance of the models. The first five principal components (PC1-PC5) were subsequently used as input variables for all ML models. This transformation ensured that all predictors were orthogonal, and representative of the main spectral variance associated with soil salinity.

Additionally, a systematic data-splitting approach has been applied, dividing the data into training (70%) and testing (30%) subsets, with 20 (folds) runs using different random seeds to assess model stability. It is worth noting that the resampling method was applied only to the training subset. The testing subset was left untouched to validate the models’ performance. The ML model’s performance was evaluated using three metrics: the coefficient of determination (R2), Mean Absolute Error (MAE), and Root Mean Square Error (RMSE). In each iteration, the metrics were calculated to determine the models’ performance and stability. Additionally, 95% Confidence Intervals (CI) for R² were also estimated to quantify the uncertainty of the models (41, 42). Additionally, statistical differences among models were evaluated using the Friedman test (43, 44), followed by the Nemenyi post-hoc test for pairwise comparisons (45). Once validated, the four ML models with median performance metrics were deployed to generate soil salinity maps, enabling a comparative spatial assessment of their predictive performance and providing valuable insights for sustainable land management.

Given the limited sample size characteristic of field campaigns in data-scarce regions, machine learning models are highly susceptible to overfitting and learning the noise in the data rather than the actual insights regarding the target variable. To mitigate this risk and enhance model robustness, a data augmentation strategy employing bootstrapping with noise is used (46). This technique involves resampling the original training data samples with replacement (bootstrapping), where the target variable (ECe) belongs to different intervals. Next, a small amount of random noise drawn from a Gaussian distribution was injected into the data. According to Aksoy et al. (15), the use of random oversampling techniques enables the model to learn a more generalizable function of the target variable, rather than memorizing individual data points.

The dataset was normalized using a standard-scaling method (Equation 1). The standardized values were subsequently rescaled to the interval of (–1, 1) to enhance numerical stability and facilitate the convergence of ML models, as recommended in previous studies (47, 48).

where xscaled  is the scaled variable, xi is the ith sample variable, xmean is the average value, and σ is the standard deviation. The mean and standard deviation of the training set were computed before the resampling phase, and were used on the resampled training subset, as well as the testing subset.

Among the variables used, most spectral bands and indices are inherently correlated, and such multicollinearity can distort regression coefficients, inflate variance, and reduce the stability of predictive models (49). PCA addresses this issue by transforming the original set of correlated variables into a smaller number of orthogonal (uncorrelated) principal components, each representing a linear combination of the original features (50). In our case, we retained enough components to preserve up to 99% of the variance in the data, resulting in five non-collinear components that capture nearly all the information in the original data. This transformation not only enhances model stability and predictive performance but also mitigates the risk of overfitting that may arise when highly collinear predictors are used in ML algorithms (49).

The algorithms used in this research fall into the category of supervised ML algorithms. The ground truth corresponding to the samples was provided as input to the model during the calibration phase. Four ML models were implemented: KNN, SVR, RF, and MLP. In addition to being widely used in literature, these models use a learning approach to identify hidden patterns in the data and correlate the input and output variables (22, 51).

To validate the performance of the models, three evaluation metrics were used, namely: R²,MAE, and RMSE, as described by Equations 2–4, respectively.

Where yi is the observed value, y^i is the predicted value, y ¯is the mean value, and n s the number of samples in the testing dataset.

The statistical analysis parameters for the target variable (ECe) and independent variables are illustrated in Figure 3. The ECe values in the study area ranged from 0.84 to 10.28 dS/m, with a standard deviation (Std) of 2.29 dS/m, reflecting considerable variation in soil salinity levels. Spatially, high EC values were concentrated downstream (southwest), while lower values were observed towards the northeast region.

Derived spectral bands and indices varied between 0 and 1.5. Among them, LDI and OLI showed the greatest variability, with a Std of 0.28 and 0.09, respectively. Meanwhile, other variables exhibited lower variability, with a Std ranging between 0.01 and 0.06.

The correlation matrix (Figure 4) revealed a strong positive relationship between the target variable EC and OLI and LDI, as well as bands 7 and 6, indicating their potential relevance to salinity variability. Conversely, BI and Int2 showed weak correlation with ECe, suggesting limited direct predictive value. However, several predictors exhibited strong intercorrelations, particularly Int1 with SI, SI1, B4, B3, B2, B7, and B6, indicating substantial multicollinearity among the spectral variables.

Given these interrelationships, PCA was applied to standardized predictors, transforming the original 13 correlated variables into orthogonal principal components (PCs). The loading pattern (Supplementary Figure S2a) indicates that PC1 carries broadly uniform positive contributions across predictors, PC2 is dominated by BI (≈0.75), PC3 by OLI-SI/LDI/OLI (≈0.58–0.59), PC4 by B6 (≈0.53) with a negative OLI contribution (≈−0.41), and PC5 by LDI (≈0.70) with a negative OLI contribution (≈−0.60). The scree and cumulative variance plots (Supplementary Figure S2–c) indicate that PC1–PC5 account for approximately 99% of the total variance.

Model training was performed using 70% (80 samples) of the dataset, with the remaining 30% (11 samples) reserved for testing. A comparison of the predicted and measured ECe values revealed that the models exhibit comparable performances during the training phase (Figure 5, Table 2). As previously mentioned, the models that yielded median scores were retained for performance representation. The KNN model yielded the lowest accuracy (R² = 0.91), while the RF, SVR, and MLP models achieved the highest accuracy (R² = 0.99). Overall, the predicted values generally aligned well with the measured ECe values, indicating satisfactory model training (Figure 6).

However, during the testing phase, performance differences become more evident (Table 3). The RF and MLP models resulted in relatively low accuracy (R² = 0.53; RMSE = 0.75 dS/m; and R² = 0.45; RMSE = 0.90 dS/m, respectively), failing to generalize to the independent test dataset, indicating overfitting, especially in the case of MLP. In contrast, the SVR model exhibited moderate predictive accuracy, with an R² of 0.59 and an RMSE of 0.62 dS/m, suggesting slight overfitting (Figure 7). On the other hand, KNN achieved the best results, with an R² score of 0.76 and an RMSE of 0.40 dS/m.

Although the overall results are relatively low (except for KNN), the models still follow the general pattern of EC variation and predict the peaks and troughs observed in the measured data (Figure 8). Additionally, sudden variations in EC lead to predictions that overshoot or undershoot the target value; these examples occur around the third (sudden decrease), fourth (sharp increase), and fifth peaks. These peaks were not reproduced by any of the models, with a greater overestimation and underestimation in MLP and RF, respectively. On the much smoother datapoints, the models tend to closely follow the EC data, with KNN producing the best visual results (8th to 12th datapoints). At the same time, MLP and SVR exhibit erratic behavior (9th-12th datapoints). Generally, the chart shows that while all models effectively capture the temporal dynamics of EC, MLP and SVR exhibit significant instability. In contrast, KNN provides the most visually consistent predictions with the observed measurements.

To strengthen the comparative evaluation of the tested models, additional statistical analyses were performed. First, the 95% CI for R² was calculated to assess the variability and robustness of model performance (Table 4). Based on results from k-fold experiments with a 20-fold split (Figure 9), the KNN model achieved the highest predictive accuracy, with a mean R² of 0.75 (95% CI: 0.73–0.77). Conversely, the MLP recorded the lowest performance (mean R² = 0.45, CI: 0.42–0.52). To assess whether differences among models were statistically significant, a Friedman test was conducted, yielding a χ²(3) statistic of 22.90 and a p-value of 0.0 (<0.05), confirming substantial differences in performance rankings. A subsequent Nemenyi post-hoc test revealed that the performances of RF, SVR, and MLP were statistically comparable (p > 0.05), with SVR and MLP being marginally different (p=0.04). In contrast, KNN differed significantly from the other models (p < 0.05) (Table 4). These findings reinforce the conclusion that KNN indeed performs best among the models used for soil salinity in the study area.

The maps generated (Figure 10) show a progressive increase in soil salinity from upstream to downstream (southwest). The spatial analysis of soil salinity classes shows distinct predictive behaviors under different ML models (Figure 11). The KNN classified the most significant portion of areas as moderately saline (49%) and showed a substantial share in non-saline soils (20%), but also produced the highest proportion of strongly saline (8.4%) and extremely saline (0.6%) areas. SVR and RF yielded comparable distributions, with moderately saline classes covering 41% and 32% of the regions, respectively, while strongly saline soils represented 6% (SVR) and 10% (RF), and extremely saline areas remained limited (1% and 0.8%). In contrast, MLP predicted the highest proportion of slightly saline soils (48.5%) and a similar share of moderately saline areas (30%), with only 1.4% being strongly saline. Overall, while KNN provided the most balanced and accurate classification, MLP and RF emphasized slightly to moderately saline conditions, and SVR maintained intermediate estimates across classes.