An Interpretable Machine Learning Framework for Unravelling the Dynamics of Surface Soil Moisture Drivers

1. Introduction

Surface soil moisture (SSM) influences climate processes by governing the distribution of precipitation into runoff, evapotranspiration, and infiltration, as well as by affecting the partitioning of incoming energy into latent and sensible heat fluxes (Seneviratne et al., 2010; Zhang et al., 2021; Han et al., 2023). SSM also plays a crucial role in the global hydrological cycle and understanding water resource management, flood generation, and climate changes at local and global scales (Sure and Dikshit et al., 2019; Yang et al., 2024). The significance of SSM becomes more pronounced in dry and semi-arid regions, where water scarcity is a prominent concern (Schwinning and Sala, 2004; Cosh et al., 2008). In these regions, changes in SSM can lead to a cascading effect, affecting the rates of groundwater recharge (Chen and Hu, 2004; Parizi et al., 2020), vegetation dynamics (D'Odorico et al., 2007), and regional climate patterns (Seneviratne et al., 2010). Hence, precise monitoring and prediction of SSM and the investigation of factors affecting it are indispensable for sustainable water resource management, risk assessment, and mitigation associated with drought and other hydrometeorological hazards in these regions (Cai et al., 2019; Lagos et al., 2020; Nikraftar et al., 2021a).

However, one significant bottleneck in SSM monitoring, particularly in arid and semi-arid regions, is the scarcity of high-quality and reliable in-situ data (Gruber et al., 2013; Everson et al., 2017; Rasheed et al., 2022). While in-situ sensors provide localized, high-quality measurements, their geographic coverage is often limited due to financial, logistical, and accessibility constraints (Xu et al., 2021a). In many countries, particularly those with vast arid and semi-arid regions, installing and maintaining in-situ networks is challenging (Gruber et al., 2013; Everson et al., 2017; Dorigo et al., 2021). Thus, validating global SSM datasets with in-situ measurements is becoming increasingly important (Jamei et al., 2020). Country-specific validation of these global datasets allows for the calibration of models to regional conditions, enhancing their reliability and applicability for resource management strategies (Dorigo et al., 2021; Guevara et al., 2021).

Understanding the factors that influence SSM is key to effective monitoring and management (Xu et al., 2021b; Rasheed et al., 2022). SSM is influenced by many factors such as precipitation (Cho and Choi, 2014; Fu et al., 2022), texture and organic matter content of soil (Cosby et al., 1984; Jawson and Niemann, 2007; Wang and Franz, 2015; Han et al., 2023), topography (Perry and Niemann, 2007), vegetation (Jiang et al., 2008), and groundwater (Meng et al., 2022). Seasonal trends, land-use changes, and evapotranspiration rates are also critical elements affecting SSM dynamics (Jung et al., 2010; de Queiroz et al., 2020). Understanding these variables and their interactions is vital for developing strategies for SSM management, particularly in arid and semi-arid regions where each drop of water counts (Wang et al., 2019). Insight into these factors can improve the precision of hydrological models and contribute to the development of adaptive management strategies for water resources (Grayson et al., 1997).

The advent of machine learning techniques has opened new avenues for hydrological modelling, particularly in the exploration of complex, non-linear relationships among factors influencing SSM (Grayson et al., 1997; Ali et al., 2015; Boueshagh and Hasanlou, 2019; Adab et al., 2020). Machine learning models can handle large datasets and account for complex interactions, thereby offering a more nuanced understanding of hydrological processes (Ahmad et al., 2010). Machine learning techniques such as Random Forest, Support Vector Machines, and Neural Networks have shown promise in hydrological applications, including SSM prediction and classification (Carranza et al., 2021). Nevertheless, numerous water scientists hesitate to adopt machine learning techniques due to their perceived "black box" nature, as they make it challenging to grasp how the model leverages input variables for prediction (Wang et al., 2022). Recently, SHapley Additive exPlanations (SHAP), a model-agnostic game-theoretic approach, has found successful application in the interpretation of machine learning models (Lundberg and Lee, 2017; Liu et al., 2022). Unlike conventional methods that only quantify the influence of input variables on the model output, SHAP can reveal whether each variable exerts a positive or negative impact on the model (Wang et al., 2022).

While previous research has utilized factors that influence SSM and employed machine learning algorithms for estimation (e.g., Carranza et al., 2021; de Oliveira et al., 2021), most studies focus on specific hydrologic conditions or operate within a limited catchment scale, resulting in models with confined applicability. The lack of diverse, large-scale data and the limited scope of prior research mean that we are yet to develop models that are broadly applicable across different hydro-climatic conditions. On the other hand, due to the spatiotemporal variability in SSM, it's crucial to understand how spatial non-stationarity of pivotal environmental factors affects SSM across different seasons. To fill this gap, we introduce a framework that combines the model-agnostic SHAP technique with a two-step clustering analysis to provide spatial interpretations for machine learning models in Iran’s 609 catchments with diverse hydro-climatic conditions. These spatial interpretations can provide guidance and policy recommendations for the effective management of SSM within 609 various catchments across Iran.

2. Materials and Methods

2.1. Study Area

Iran with an area of approximately 1,648,195 km² is located in the Middle East and faces challenges related to water resource scarcity (Nikraftar et al., 2023). Iran is typically characterised by an arid and semi-arid climate, with an average annual precipitation of approximately 250 mm (Rahmani et al., 2016). As indicated by climate classification (Figure 1), a significant portion of the country falls under the category of a warm-dry climate, naturally resulting in SSM deficits (IMO, 2019). In recent years, SSM deficits in Iran have been exacerbated by factors that include climate change, excessive groundwater extraction, mismanagement of surface water, and inefficient irrigation practices (Ashraf et al., 2021). Such SSM shortages can lead to detrimental environmental consequences, including desertification, soil degradation, dust storms, wind erosion, and the degradation of air and water quality (Sivakumar and VStefanski, 2007). Therefore, it is imperative to monitor SSM and to study the factors influencing it in Iran to ensure effective environmental management (Gheybi et al., 2019; Jamei et al., 2020). Nonetheless, SSM remains inadequately monitored in numerous regions of Iran, and the existing measurements (as depicted in Figure 1) lack sufficient temporal and spatial resolution (Rahmani et al., 2016).

Under these circumstances, remote sensing data can offer a viable solution for monitoring SSM in Iran (Ranjbar Saadatabadi et al., 2021). Consequently, it is essential to validate and assess global SSM products using in-situ stations to gauge their real-world effectiveness before using them across all climatic regions in Iran. In this study, after validation of the SSM products, the environmental factors affecting SSM were investigated across 609 catchments in Iran. The total drainage area of these catchments is 1,648,195 km² (IWRMC, 2023). These catchments are distributed across diverse climates, spanning from cold-dry to warm-humid, as illustrated in Figure 1. They exhibit different topographic characteristics, geological compositions, soil varieties, and vegetation, leading to diverse hydrological conditions throughout the country. Additionally, the land cover map obtained from ESA (ESA-World Cover, 2020) with a 10 m resolution reveals that bare/sparse vegetation, grasslands, croplands, and forests constitute the predominant land cover in the studied catchments, accounting for 64.2%, 18.9%, 11.8%, and 1.70%, as given in Table S1. Therefore, validating SSM products through in-situ stations and utilising an appropriate dataset to investigate the environmental factors affecting SSM in this diverse environment facilitates the advancement of sustainable SSM management practices.

2.2. Datasets

2.2.1. SSM Data

The in-situ SSM data from agrometeorological stations were gathered from the Iran Meteorological Organization, IMO (IMO, 2023) for the period spanning February 1, 2006, to March 31, 2023. Unfortunately, SSM data is not widely available in numerous regions of Iran and have several gaps (Rahmani et al., 2016; Gheybi et al., 2019). Hence, a total of 42 stations were chosen throughout Iran, distributed as follows: 13 stations in cold-dry climates, 2 in cold-humid, 11 in warm-dry, 3 in warm-humid, 8 in moderate-dry, and 5 in moderate-humid (Figure 1). Considering previous studies (e.g., Colliander et al., 2017; Xu, 2021a; Tian & Zhang; 2023), here we utilized three primary SSM datasets, namely Soil Moisture Active Passive (SMAP), Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2), and Climate Forecast System, Version 2 (CFSv2), to validate SSM at a depth of 0-5 cm using in-situ data (Table 1).

The SMAP satellite mission by NASA was launched on January 31, 2015, to globally map soil moisture and the freeze/thaw state of landscapes (Wu et al., 2020). While the satellite initially utilized both an L-band radar and an L-band radiometer, the radar instrument encountered a failure after approximately 11 weeks of operation, and the production of soil moisture data continues using only the radiometer data (Colliander et al., 2017). SMAP measurements provide direct sensing of SSM in the upper 5 cm of the soil (Kimball et al., 2018; Xu et al., 2021). In this study, we employed the SMAP Level 4 passive product (Reichle et al., 2019), which features a spatial resolution of 9 km (Table 1).

MERRA-2, generated by NASA Global Modeling and Assimilation Office (GMAO), stands as the most recent atmospheric reanalysis covering the modern satellite era, offering global reanalysis data spanning from 1980 to the present (Gelaro et al., 2017). Reichle et al. (2017) validated MERRA-2 against in-situ measurements in North America, Europe, and Australia, demonstrating that its performance is slightly superior to that of ERA-Interim/Land. In this study, we utilised MERRA-2, which featured a spatial resolution of 56 km by 70 km (Table 1). CFSv2 was made operational at the National Centers for Environmental Prediction (NCEP) in March 2011 (Saha et al., 2014; Dirmeyer and Halder; 2016). The soil moisture dataset within the CFSv2 encompasses four layers (at 5, 25, 70, and 150 cm depths) and boasts a spatial resolution of approximately 22 km (Saha et al., 2011). It's worth mentioning that while SMAP directly measures SSM, MERRA-2 and CFSv2 rely on model-based reanalysis data generated through physical models, observations, and data assimilation. These three SSM datasets were resampled to 0.25º × 0.25º by bilinear interpolation (Xu et al., 2021) to ensure consistency across the datasets.

2.2.2. Factors Influencing SSM

We selected the candidate factors that influence SSM based on previous studies (e.g., Chen et al., 2004; Patel et al., 2009; Zhao et al., 2017; Raduła et al., 2018; Han et al., 2023) and data availability (Table 2, Figure S1). We used the catchment-averaged monthly data of the following parameters: precipitation, potential evapotranspiration, solar radiation, wind speed, normalized difference vegetation index (NDVI), and groundwater table depth. Additionally, time-invariant catchment attributes such as, distance from water bodies, clay fraction, organic matter fraction, elevation, and topography roughness index were included (Table 2). The mean monthly precipitation for the studied catchments was computed using the Radial Basis Function (RBF) interpolation method (Du, 2008) in Python. This was based on daily precipitation data collected from 422 synoptic stations between 2015 and 2023 from IMO (IMO, 2023). The mean monthly potential evapotranspiration for the studied catchments was determined using MODIS global evapotranspiration product (i.e., MOD16A2, Running et al., 2017) at a 500 m resolution.

Solar radiation and wind speed data were obtained from the ERA5-Land dataset (Muñoz-Sabater et al., 2021) at a spatial resolution of 11 kilometers, while an analysis of vegetation dynamics was conducted using the NDVI index and Sentinel-2 data (Copernicus, 2023) with 10 m resolution. The NDVI proves to be a suitable index for detecting vegetation changes, particularly within arid and semi-arid regions (Parizi et al., 2021). To generate a time series of groundwater table depth, we utilised monthly data from 11,003 observation wells collected by the Iran Water Resources Management Company (IWRMC, 2023). We derived the average monthly groundwater table depth for each catchment using the RBF interpolation method. The mean distance of each catchment from water bodies was determined by utilising the water bodies' data (NOAA, 2023) and the Euclidean Distance method in ArcGIS (ESRI, 2013). The SoilGrids250m dataset (Hengl et al., 2017) was employed to extract the clay and organic matter fractions within the upper vadose zone. Finally, we calculated the mean elevation and topography roughness index for the studied catchments using an ALOS DEM with a 30 m resolution (Tadono et al., 2016) and Focal Statistics tools in ArcGIS. All factors influencing SSM were extracted using the Google Earth Engine platform (Amani et al., 2020), Python and ArcGIS software. Additional details regarding the factors influencing SSM and their correlation are provided in the supplementary material.

2.3. Methods

2.3.1. Statistical Metrics

The performance of global SSM products was assessed using the Root Mean Squared Error (RMSE), Relative Bias (RBias), Kendall’s Tau (τ), and Kling-Gupta efficiency (KGE′) (Tang et al., 2020; Saemian et al., 2021; Fahrudin et al., 2020; Parizi et al., 2022):

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{(\mathrm{T}\mathrm{a}{\mathrm{r}}_{\mathrm{i}}-\mathrm{R}\mathrm{e}{\mathrm{f}}_{\mathrm{i}})}^2}{\mathrm{N}}}}$$

(1)

$$\mathrm{R}\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{N}}(\mathrm{T}\mathrm{a}{\mathrm{r}}_{\mathrm{i}}-\mathrm{R}\mathrm{e}{\mathrm{f}}_{\mathrm{i}})}{{\sum }_{\mathrm{i}=1}^{\mathrm{N}}\mathrm{R}\mathrm{e}{\mathrm{f}}_{\mathrm{i}}}}$$

(2)

$$\tau ={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\mathrm{T}\mathrm{a}{\mathrm{r}}_{\mathrm{i}}-\mathrm{T}\mathrm{a}{\mathrm{r}}_{\mathrm{j}})\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\mathrm{R}\mathrm{e}{\mathrm{f}}_{\mathrm{i}}-\mathrm{R}\mathrm{e}{\mathrm{f}}_{\mathrm{j}})}{\mathrm{n}(\mathrm{n}-1)}}$$

(3)

where N represents the number of samples, Ref is the reference values (i.e., in-situ data), and Tar is the target values (i.e., SMAP, MERRA-2, and CFSv2 data) for each record (i). Furthermore, we utilized the KGE′ statistic, initially introduced by Gupta et al. (2009) and subsequently modified by Kling et al. (2012). KGE balances the contributions of correlation, bias and variability terms as follows (Tang et al., 2020):

$$\mathrm{K}\mathrm{G}\mathrm{E}=\mathrm{l}-\sqrt{{(\mathrm{r}-\mathrm{l})}^2+{(\beta -\mathrm{l})}^2+{(\gamma -\mathrm{l})}^2}$$

(4)

$$\beta ={\displaystyle \frac{{\mu }_{\mathrm{T}\mathrm{a}\mathrm{r}}}{{\mu }_{\mathrm{R}\mathrm{e}\mathrm{f}}}}$$

(5)

$$\gamma ={\displaystyle \frac{{\mathrm{CV}}_{\mathrm{T}\mathrm{a}\mathrm{r}}}{{\mathrm{CV}}_{\mathrm{Ref}}}}={\displaystyle \frac{{\sigma }_{\mathrm{T}\mathrm{a}\mathrm{r}}/{\mu }_{\mathrm{T}\mathrm{a}\mathrm{r}}}{{\sigma }_{\mathrm{Ref}}/{\sigma }_{\mathrm{Ref}}}}$$

(6)

where r represent the correlation coefficient between Ref and Tar datasets, γ is the variability ratio, β is the bias ratio, μ is the mean SSM, CV is the coefficient of variation, and represents is the standard deviation. The better values have higher KGE (Saemian et al., 2021).

2.3.2. Random Forest

The Random Forest (RF) model was initially introduced by Leo Breiman in 2001 (Naghibi et al., 2016). The RF algorithm does not necessitate any alterations, conversions, or modifications to the input data, and it autonomously handles missing values (Breiman et al., 2012; Pouyan et al., 2021). An RF model comprises a multitude of decision trees designed to be as uncorrelated as possible (Kaiser et al., 2022). To create an uncorrelated collection of trees, the RF employs bagging and feature randomisation during the construction of the decision trees. This implies that each tree is trained on a random sample drawn from the training set with replacement, a technique known as bootstrapping (Amini et al., 2022). Additionally, each tree is limited to a random subset of the available features (Breiman, 2001).

One of the advantages of the RF method, as highlighted by Tyralis et al. (2019), is that it yields consistent predictions while diminishing variance without augmenting prediction bias. In this study, we used the RandomForestRegressor from the machine learning library sklearn (Pedregosa et al., 2011; Buitinck et al., 2013) for the Python programming language. For enhanced convergence speed and to mitigate the impact of local extremes on training, the input variables undergo normalization to a range of 0.1 to 0.9 before the training process. As some models encounter issues when normalized between 1 and 0, we chose to normalize the input variables using the adjusted min-max method within the range of 0.1 to 0.9 (Gökhan Aksu et al., 2019).

The hyperparameters were tuned using a 10-fold cross-validation approach on the training set. This involved dividing the training set into ten subsets, with one-tenth serving as the test sample in each fold, while the remaining data was employed for constructing the regression trees. In each fold, the test sample comprised randomly selected points from the overall training set that had not been included as a test sample in previous folds. The comparison of RMSE values for each hyperparameter combination across the 10-fold cross-validation provided an assessment of model performance (Carranza et al., 2021). To maintain a distinct validation set, a cross-validation scheme is a precaution against model overfitting (Lever et al., 2016; Kaiser et al., 2022). We employed a random selection process for the training dataset to allocate 80% of the data for training purposes, while the remaining 20% was reserved for testing.

2.3.3. SHAP

SHAP method represents a model-agnostic game-theoretic technique for interpreting machine learning models (Lundberg and Lee, 2017). Unlike conventional methods that only quantify the influence of input variables on the model output, SHAP can disclose whether each variable exerts a positive or negative impact on the model (Wang et al., 2022). In other words, SHAP can analyse an individual prediction by considering it as a composite result of the combined effects of each input variable on the output value (i.e., the predicted value). This approach allows users to gain insight into the magnitude of influence and how each input variable affects the output value (Lundberg and Lee, 2017, Broeck et al., 2022). Using a pre-trained machine learning model denoted as M and a set of input variables x={x₁, …, x_q}, SHAP employs an explanation model E to ascertain the individual influence of each variable on the behaviour of model M (Liu et al., 2022). SHAP express as:

$$\mathrm{E}={\varphi }_0+{\sum }_{\mathrm{i}=1}^{\mathrm{q}}{\varphi }_{\mathrm{i}}{\mathrm{t}}_{\mathrm{i}}$$

(7)

$${\varphi }_{\mathrm{i}}(\mathrm{M},\mathrm{x})={\sum }_{\mathrm{t}\subseteq \mathrm{x}}{\displaystyle \frac{\left|\mathsf{t}\right|!(\mathsf{q}-|\mathsf{t}|-1)!}{\mathsf{q}!}}\left[\mathrm{M}\right(\mathsf{t})-\mathrm{M}(\mathsf{t}⧵\mathrm{i}\left)\right]$$

(8)

where q represents the number of input variables, t is the variable simplification, ɸ_i ∈ R represents the contribution of each variable to the machine learning model, and \ is the difference-set notation for set operations (Wang et al., 2022; Liu et al., 2022). In this study, the Python SHAP library was employed to assess feature importance in conjunction with the RF model.

2.3.4. Cluster Analysis

A two-step cluster analysis was conducted on the outcomes of the SHAP model to enhance our spatial understanding of the effects of factors influencing SSM within the studied catchments. In other words, employing a two-step cluster analysis can offer spatial interpretations for the factors that affect SSM, utilizing the outcomes generated by RF and SHAP. We follow the methodology for the two-step cluster analysis outlined by Fahy et al. (2019). Only a summary of the technique is provided here, and rather, the reader is directed to Fahy et al. (2019) for a detailed description of the approach. The method involves two steps: (1) whole records are probed by distance to construct a classification tree, whereby records in the same tree node are most similar (Qin et al., 2019); (2) nodes are classified using the cohesion technique and clustering results are evaluated using the Bayesian information criterion (BIC) or the Akaike information criterion (AIC), which determine the structure of the final cluster (Chiu et al., 2001; Satish & Bharadhwaj, 2010).

3. Results and Discussion

3.1. Performances of SSM Products

Figure 2 shows the SSM validation results based on 42 in-situ stations from April 2015 to March 2023. The findings suggest that among the evaluated datasets, SMAP stands out with the highest median values for τ and KGE (0.740 and 0.690), and it also exhibits the lowest median values for RMSE and RBias (0.068 and 0.030). Following SMAP, the MERRA-2 product demonstrates median values of τ and KGE at 0.684 and 0.604, along with median RMSE and RBias values of 0.085 and 0.034, respectively. CFSv2 shows the lowest median values for both τ and KGE (0.550 and 0.500), and it also has the highest median values for RMSE and RBias (0.113 and 0.059) among the three products, as shown in Figure 2. These findings demonstrate that SMAP outperforms the other datasets in estimating SSM in Iran.

The performance of SMAP aligns with the findings of the study conducted by Chen et al. (2018), which reported that SMAP has a global average anomaly correlation of 0.76. Xu et al. (2021a) evaluated eight global root zone soil moisture products (0-1 m depth) across the globe. Their findings indicated that SMAP, MERRA-2, JRA-55, and ERA-5 consistently showed a stronger correlation with in-situ root zone soil moisture measurements compared to GLDAS, NCEP R1, and NCEP R2. Colliander et al. (2017) validated SMAP SSM using core validation sites. They reported that the SMAP radiometer-based SSM product meets its expected performance, achieving an unbiased root mean square error of 0.04 m³ /m³ for volumetric SSM. It's noteworthy that global evaluations have not incorporated the in-situ SSM data from Iran.

There have been limited studies that have focused explicitly on evaluating SSM in Iran. These studies either concentrate on a particular local area, such as the Lake Urmia Basin (Maleki et al., 2019; Saeedi et al., 2021) or cover a short period in Iran (i.e., 2015-2016, as demonstrated in Gheybi et al., 2019). Additionally, some studies were solely concerned with validating a single product (e.g., Jamei et al., 2020, Jamei et al., 2022, and Amini et al., 2023). For example, Gheybi et al. (2019) validated SSM products from SMAP, SMOS, and AMSR2, using 23 in-situ stations in Iran from 2015 to 2016. Their results pointed to SMAP as the best-performing satellite-based product. Also, Jamei et al. (2022) stated that SMAP has a strong capacity for SSM data retrieval in Iran.

3.2. Spatial-Temporal Pattern of SSM

We calculated SSM for 609 catchments in Iran using the SMAP dataset, which demonstrated its optimality (Figure 3). Figure 3 illustrates the mean catchments-averaged daily SSM across 609 studied catchments from April 2015 to March 2023, delineated by different seasons and based on the SMAP dataset. The mean SSM for different seasons reveals distinct patterns, with winter having the highest median at 0.175 m³/m³, followed by spring at 0.160 m³/m³, while autumn and summer exhibit lower median of 0.096 m³/m³ and 0.081 m³/m³, respectively. These findings indicate notable seasonal variations in SSM within the studied catchments. Fakharizadehshirazi et al. (2019) stated that many catchments in Iran lack natural moisture, especially during the summer, leading to a heightened demand for irrigation in agriculture during this season. According to Garcia-Estringana et al. (2013), differences in soil moisture levels between dry (summer) and wet (winter) conditions are more pronounced in the upper surface layers (0–20 cm) when compared to deeper layers. Figure 3 also demonstrates that catchments with SSM exceeding 0.20 m³/m³ are predominantly located in northern, northwestern, western, and southwestern Iran, primarily within regions characterized by cold-humid and moderate-humid climates. In contrast, catchments with SSM below 0.05 m³/m³ are primarily concentrated in central and southeastern Iran that have warm dry climates.

3.3. SSM in Different Land Covers

Numerous studies have shown that changes in SSM exhibit varying characteristics when subjected to different land cover conditions (e.g., English et al., 2005; Jin et al., 2018; Feng et al., 2018; Zhou et al., 2023). The boxplots of SSM within the six primary land cover categories in Iran, as determined by the land cover map of ESA with a resolution of 10 m (ESA-World Cover, 2020) for the period spanning April 1, 2015, to March 31, 2023, are displayed in Figure 4a. The results reveal significant variations in median SSM across different land cover. Specifically, we found that forests exhibited the highest SSM with a median value of 0.180 m³/m³, followed by grasslands at 0.148 m³/m³, croplands at 0.148 m³/m³, shrubland at 0.117 m³/m³, built-up areas at 0.107 m³/m³, and bare/sparse vegetation at 0.093 m³/m³. Janani et al. (2022) studied SSM in various land-use patterns along the lower Bhavani River in India. They concluded that SSM is higher in forested areas compared to fallow land and built-up areas. Figure 4a also indicates that grasslands exhibit the highest soil moisture diversity, while bare/sparse vegetation shows the least diversity across Iran.

Figure 4b indicates a time series of SSM for various land cover in Iran, spanning from April 1, 2015, to March 31, 2023. The results show that in the early months of 2019, especially in March and April, SSM reached higher levels compared to the same months in the preceding and subsequent years (Figure 4b). The analysis of the precipitation time series in Iran (as shown in Figure 4b) for the years 2015-2023 reveals that the increase in SSM during those particular months is due to the increase in precipitation. This finding aligns with previous research. Nikraftar et al. (2021b) demonstrated that a significant increase in precipitation during the early months of 2019 led to a rise in the water level of Lake Urmia in northwestern Iran. Khosravi et al. (2020), Sadeghi et al. (2021), and Parizi et al. (2022) have documented that due to the heavy and unprecedented precipitation occurring between mid-March and April 2019, widespread flooding events affected 25 out of the 31 provinces in Iran. These events resulted in more than 77 human fatalities and inflicted approximately US$ 2.2 billion in damages.

3.4. RF and SHAP

Figure 5 illustrates the comparison of the SSM data obtained by SMAP vs RF model in the testing phase for different seasons. The findings indicate that RF can yield SSM estimations with R² values of 0.89, 0.83, 0.70, and 0.75 for the winter, spring, summer, and autumn seasons, respectively. The decrease in R² values observed in the summer can be attributed to the significant shortage of SSM during this season compared to others, notably winter (as is shown in Figure 3). The good performance of the RF method in estimating SSM in this study aligns with the findings of previous research. For example, in the semi-arid region of West Khorasan-Razavi province in Iran, Adab et al. (2020) utilized several machine learning algorithms for SSM estimation. Their study concluded that the RF method provided the most precise results. Fathololoumi et al. (2020) conducted a study comparing spectral and spatial-based approaches to map local SSM variations in the Balikhli-Chay watershed in northwestern Iran. Their findings revealed that the RF approach outperformed others, demonstrating the highest level of performance in SSM modelling.

Machine learning techniques are often considered as black-box models, which limits their interpretability regarding the process of making predictions (Choubin et al., 2018). SHAP offers a way to understand the influence of each feature on the model's outputs (Ekmekcioglu et al., 2022). Recent research by Boueshagh et al., (2025) used RF and SHAP to investigate the impact of various environmental factors, including SSM, on spatial statistics of satellite-derived surface brightness temperatures to improve snowpack estimations in regional and global scales.

Figure 6 indicates summary plots of SHAP analysis for various seasons. These plots represent features along the vertical axis (y-axis) and SHAP values along the horizontal axis (x-axis). In these figures, each data point represents a SHAP value associated with its corresponding feature. Red points signify high feature values, while blue points represent low ones. When red points are positioned to the right, they signify an increase in a factor leading to an increase in SSM, whereas their presence on the left indicates that an increase in the factor results in decreased SSM. Conversely, if blue points appear to the right, a reduction in a factor increases the SSM, and when they're situated to the left, it signifies that a decrease in the factor leads to decreased SSM.

The findings reveal that the primary factors influencing SSM vary from season to season (Figure 6).

In the winter season, the key factors that exert the most influence on SSM are precipitation, distance from water bodies, solar radiation, clay fraction, potential evapotranspiration, and elevation, respectively (Figure 6a). Winter is typically associated with increased rainfall in many of Iran’s catchments (Tabari and Hosseinzadeh Talaee, 2011; Javari, 2016). Previous research has demonstrated a direct contribution of precipitation to SSM (e.g., He et al., 2012; Feng and Liu, 2015; Rascón-Ramos et al., 2021). The proximity of the studied catchments to water bodies, such as the Caspian Sea and the Persian Gulf, plays a vital role in determining SSM during the winter. For example, the distance from water bodies can influence relative humidity, which consequently impacts SSM (Carranza et al., 2021; Du et al., 2021). Solar radiation with inverse impact on SSM tends to be lower during the winter due to shorter days and reduced sunlight, which affects the rate of moisture evaporation from the soil as reported by Wenwu et al. (2018). So, lower solar radiation in winter can help maintain higher SSM.

As the fourth important factor influencing SSM in winter, clay fraction can enhance its water-holding capacity (Vachaud et al., 1985; Ojha et al., 2014; Han et al., 2022). Clay soils retain moisture more effectively, contributing to higher SSM levels. One of the other factors affecting SSM is potential evapotranspiration, which tends to be lower during the winter in the studied catchments due to cooler temperatures. This reduction in potential evapotranspiration can contribute to the SSM preservation. Carranza et al. (2021) highlighted potential evapotranspiration as a key determinant in estimating root zone soil moisture within the Raam catchment in the Netherlands. As the sixth parameter, elevation with inverse impact can play a role in determining SSM during winter. Catchments at different elevations may experience variations in temperature and slope, which can influence SSM retention (Pellet and Hauck et al., 2017; Xu et al., 2021b). The analysis of the published data on seven potential factors influencing the temporal stability of soil water content indicated that the influence of these factors appears to be interconnected rather than solely driven by a single dominant factor (Vanderlinden et al., 2012).

As the spring season commences and plants and trees in Iran start to grow, the significance of NDVI, distance from water bodies, clay fraction, and organic matter fraction becomes more pronounced compared to precipitation (Figure 6b). In other words, during spring, land cover and soil characteristics take on a greater importance compared to winter. On the other hand, increased plant growth and the decomposition of organic materials, driven by increased temperatures, can lead to an augmentation in soil organic matter content in spring. This, in turn, notably influences the soil's ability to retain water. Based on Han et al. (2023), an escalation in organic matter content results in an enhanced water-holding capacity due to the inherent affinity of organic matter for water. Fakharizadehshirazi et al. (2019) investigated long-term spatiotemporal variations in SSM and vegetation indices across Iran. Their findings concluded that NDVI exerts a noteworthy influence on the spatiotemporal variations of SSM. Ghasemloo et al. (2022) estimated agricultural farm SSM using spectral indices and demonstrated that NDVI and land surface temperature possess substantial potential for extracting valuable SSM information.

With the significant decrease in precipitation in summer, the dominant factors affecting SSM are proximity to water bodies, clay fraction, potential evapotranspiration, NDVI, organic matter fraction, and elevation (Figure 6c). These findings reveal that the influence of the clay fraction on SSM reaches its peak during this dry season, surpassing the impact seen in other seasons. Catchments with a substantial clay content can effectively retain SSM during this period. As autumn arrives and precipitation commences, there is a shift in the hierarchy of factors influencing SSM fluctuations, with precipitation emerging as the predominant factor impacting SSM during this season (Figure 6d). It's worth noting that additional factors, such as groundwater table depth, wind speed, and topography roughness index, have an inverse effect on SSM. These factors have a relatively minor impact on SSM in comparison to other factors. The low influence of the groundwater table depth on SSM can be attributed to the average depth of the groundwater table in the studied catchment, which stands at around 31 m (as indicated in Figure S1f). The impact of groundwater depth on SSM becomes noticeable primarily in catchments with shallower groundwater tables, especially in winter.

3.5. Cluster Analysis

The SHAP model yielded a large number of coefficients for the 609 studied catchments, presenting a challenge in terms of interpretation. Hence, there is a requirement for a methodology to categorize the non-stationarity results of factors influencing SSM across seasons in Iran. In this study, cluster analysis was used to process SHAP coefficients. Following the clustering process, there is a discernible resemblance in the factors affecting catchments grouped within the same category. When comparing these categories, the results inherently manifest substantial distinctions in the various factors across catchments. In SPSS software, the two-step clustering algorithm automatically assesses the suitability of segmenting catchments into multiple categories by considering lower BIC values. Once it automatically determines the optimal number of clusters, it becomes possible to gain more precise insights into the distinctions among these categories. The results of clustering enable us to comprehend the extent to which each factor influences SSM in every cluster. In this study, catchments were automatically classified into five types, as the clustering yielded the lowest BIC index at this point. Table 3 presents the average SHAP values of the six primary factors for every season and class. Figure 7 displays the spatial distribution of the clustering results for various catchments. Coefficients with significant absolute values carry a more pronounced influence on the SSM. These figures serve as valuable tools for comprehending each catchment type's spatial arrangement and identifying the factors that exert the most substantial influence on each catchment's SSM.

The clustering results for winter are shown in Figure 7a,e. The first catchment types are located mainly in northwestern and western Iran. In these catchments, the factors that most significantly affect SSM are, in respective order, precipitation, clay fraction, solar radiation, potential evapotranspiration, distance from water bodies, and elevation, as illustrated in Figure 7a,e and detailed in Table 3. Akbari Majdar et al. (2018) investigated the spatial and temporal variations in SSM with respect to topographic and meteorological factors in Ardabil province, which falls into the first catchments category of this study. Their research underscored a significant correlation between SSM and variables such as precipitation. The second type of catchments is primarily found in central, eastern, and northeastern Iran. The most influential factors affecting SSM in these catchments are precipitation, distance from water bodies, solar radiation, clay fraction, potential evapotranspiration, and elevation. The third catchment category is predominantly located in southern and northern Iran, in proximity to the Persian Gulf and Caspian Sea. Within this catchment type, the most pivotal factor affecting SSM is the distance from water bodies, with precipitation, clay fraction, solar radiation, elevation, and potential evapotranspiration following in significance. In the fourth catchment category, similar to the third catchment type, the primary factor influencing SSM is the distance from water bodies, albeit with an inverse effect. This is followed by solar radiation, precipitation, potential evapotranspiration, elevation, and clay fraction (Table 3). The fifth catchment category is distributed in southeastern Iran. In order of influence, the most significant factors on SSM are precipitation, clay fraction, solar radiation, distance from water bodies, potential evapotranspiration, and elevation.

The clustering results for the spring season are shown in Figure 7b,f. In this season, vegetation and organic matter fraction have a higher impact on SSM than in winter. The first catchment type is located mainly in north, northwestern and western Iran. Within these catchments, the factors exerting the most substantial influence on SSM are, in their respective order of impact: NDVI, organic matter fraction, clay fraction, precipitation, solar radiation, and distance from water bodies, as demonstrated in Figure 7 and expounded upon in Table 3. The second catchment type is mainly distributed in southern, southwestern, and north Iran. Factors that most significantly affect SSM are distance from water bodies, NDVI, clay fraction, organic matter fraction, solar radiation, and precipitation. In the spring of 2020, Bandak et al. (2023) collected 394 surface soil samples in the Gorgan province located northern Iran, a region which is located in the second catchment type of this study. The findings of their research demonstrated a strong correlation between NDVI and SSM. The third catchment category is primarily located in southern and central Iran. Within this catchment category, the most dominant factor affecting SSM is the NDVI, with organic matter fraction, precipitation, clay fraction, distance from water bodies, and solar radiation subsequently ranking in significance. In the fourth catchment category, the predominant factor influencing SSM is the distance from water bodies. This is followed by clay fraction, organic matter fraction, precipitation, solar radiation, and NDVI (Table 3). The fifth catchment category is distributed in southeastern and eastern Iran. The most significant factors on SSM are NDVI, clay fraction, precipitation, organic matter fraction, solar radiation, and distance from water bodies.

Figure. 7c.g and Table 3 indicate the clustering results for the summer season. During this season, the influence of distance from water bodies and clay fraction on SSM is more pronounced compared to other seasons. The first catchment type is located in northern, northwestern and northeastern Iran around the Caspian Sea. In these catchments, the most significant influences on SSM are associated with potential evapotranspiration, proximity to water bodies, clay fraction, NDVI, organic matter fraction, and elevation (Fig 7c,g and Table 3). The second catchment type is mainly located in southern and southeastern Iran around the Persian Gulf and the Gulf of Oman. In this class, distance from water bodies is the most important factor, followed by clay fraction, potential evapotranspiration, organic matter fraction, NDVI, and elevation. In the third and fourth catchment categories, the hierarchy of factors' influence closely mirrors that of the second catchment type, except for the notable difference that NDVI exerts a more significant influence than the organic matter fraction. The fifth catchment category is primarily located in southeastern and eastern Iran. The most influential factors are distance from water bodies, clay fraction, potential evapotranspiration, elevation, NDVI, and organic matter fraction (Table 3).

The clustering results for autumn are shown in Figure 7d,h and Table 3. As precipitation begins in Iran, it becomes the predominant factor affecting SSM in most classes. The first catchment category is primarily located in central and eastern Iran. Within these catchments, the factors exerting the most substantial influence on SSM are precipitation, distance from water bodies, potential evapotranspiration, elevation, solar radiation, and clay fraction, as demonstrated in Figure 7d,h and described in Table 3. The second catchment type is mainly distributed in southwestern and western Iran. The factors that most significantly affect SSM are precipitation, potential evapotranspiration, clay fraction, distance from water bodies, solar radiation, and elevation. The third catchment category is primarily located in southern (near the Persian Gulf) and northeastern Iran. In this catchment type, precipitation is the most influential factor on SSM, followed by the distance from water bodies, elevation, potential evapotranspiration, clay fraction, and solar radiation. In the fourth catchment category, the primary influencing factors on SSM are precipitation and clay fraction, potential evapotranspiration, distance from water bodies, solar radiation, and elevation (Table 3). Finally, the fifth catchment category is distributed in northern Iran near the Caspian Sea. The most significant factors on SSM are potential evapotranspiration, precipitation, distance from water bodies, solar radiation, clay fraction, and elevation.

4. Conclusion

In this study, we unravelled the impact of the spatial non-stationarity of the critical environmental factors influencing SSM. To this end, we have introduced a framework that combines the SHAP technique with a two-step clustering analysis to provide spatial interpretations for machine learning models such as RF. Considering poor in-situ data of SSM in Iran, we initially validated the global SSM datasets, which include SMAP, MERRA-2, and CFSv2 at a depth of 0-5 cm, against in-situ data. The main conclusions are as follows:

1) Results of dataset validation of SSM demonstrated that among the datasets, SMAP exhibited the highest median correlation and the lowest median root mean square error when compared to in-situ stations. Hence, it is recommended for applications such as hydrological modelling, water resources management, and drought monitoring in Iran, where SSM data is scarce.

2) Investigation of SSM in the different land cover in Iran revealed significant variations in SSM across land cover. So that, forests and bare/sparse vegetation regions exhibited the highest and lowest SSM with a median value of 0.180 m³/m³ and 0.093 m³/m³, respectively. These findings highlight the importance of understanding the spatial distribution of SSM in different land cover types, which can have implications for various environmental and ecological processes.

3) Findings indicated that the RF model can offer SSM estimates with R² values of 0.89, 0.83, 0.70, and 0.75 for the winter, spring, summer, and autumn seasons, respectively. This finding highlights the importance of seasonal investigation of SSM. With a shortage of SSM in the dry season, machine learning methods can't model SSM with very high accuracy.

4) The findings of the SHAP model and two-step cluster analysis indicated that winter and autumn SSM is primarily influenced by climatic factors. In contrast, spring and summer SSM is predominantly affected by vegetation and soil characteristics. These findings highlight the dynamic nature of SSM and how it is influenced by different environmental factors that change with the seasons. Understanding these seasonal variations is essential for effective SSM management and prediction in various regions.