Application of the artificial intelligence approach and remotely sensed imagery for soil moisture evaluation

Surface soil moisture (SM) is an important contributor to hydro-climatic processes, having critical roles to play in disaster prediction (e.g., flooding and drought), environmental monitoring (such as dust storms and erosion), and hydrological applications. Soil water content can affect the physio-biologic components of the biosphere and link the surface of the earth to the atmosphere via surface energy and moisture flux. SM is a source of water for the atmosphere through evapotranspiration from the land, which includes vegetable transpiration and bare soil evaporation. Moreover, SM conditions can affect the hydrological pattern of the land surface by controlling the infiltration capacity of the soil and portioning rainfall to runoff. The ecohydrology which focuses on the linkage between vegetation -water - climate relationships has been found to have a complicated dependency on the availability of the SM dynamics (Garcia-Estringana et al. 2013; Mulebeke et al. 2013). All these processes are highly characterized by the nonlinear behavior of SM and complex feedback mechanisms. Therefore, quantified conditions of SM are influential inputs for modeling agricultural, hydro-climatological, and meteorological attributes.

A set of components governs the dynamics of SM on the land surface with diverse time and spatial scales. Variations in both weather and climate are, therefore, influenced by the condition of SM. Reynolds (1970) classified the static (e.g., soil texture and topography) and dynamic (e.g., precipitation and vegetation) controlling elements of SM. An assessment of SM depends on the condition of the associated variables. Many of these elements are interrelated and range spatially and/or temporally, rendering the recognition of relationships among SM patterns and its riding variables complicated. Landscape elements, including topography, vegetation, and land use, are spatial and temporal governing elements of SM. SM's spatial variations are notably related to the features of the terrain (e.g., slope, elevation, and topographic wetness index). Therefore, the attributes of the terrain were utilized to estimate the parameters of SM patterns through regression, geospatial, and hydrological modeling in some previous research studies (e.g., see Western et al. 1999, 2004; Adab et al. 2020; Li et al. 2021). Additionally, the impacts of vegetation cover (e.g., type and distribution) on SM variation have been noted in various research studies. Also, the effects of spatial attributes on vegetation (normally interpreted from remote sensing images) have been utilized to generate SM patterns (Mohanty et al. 2000; Hupet & Vanclooster 2002).

Generally, long-term time series of SM can spatially detect changes in the water cycle associated with the weather or hydrological condition. In a large study area, the wide variety of networks and gauging SM remain restricted and furthermore, obtaining reliable approximations from field measurements is a challenging task due to excessive variability and the lack of correlations among the parameters. In several applications of SM, a wide variety of satellite-based products have the potential to help hydrologists to measure the SM condition in large regions. Since the 1970s, some remote-sensing strategies have been developed for analyzing and mapping SM by measuring specific regions of the electromagnetic spectrum ranging from optical to microwave areas (Musick & Pelletier 1988; Engman 1991; Wang & Qu 2009). As remote sensors cannot measure SM content directly, mathematical-based methods that can explain the relationship of measured signals and SM content need to be extracted. Microwave remote sensing strategies that include the Advanced Microwave Scanning Radiometer-Earth staring at system (AMSR-E) on-board the Aqua satellite for pc (from 2002), Soil moisture and ocean salinity satellite (SMOS from 2009), Multi-frequency Scanning Microwave Radiometer (MSMR from 1999), and Soil Moisture Active Passive (SMAP) (from January 2015) are currently operational, producing daily satellite records globally. Although these strategies provide many techniques to gauge SM on a large scale, they have almost low resolution (typically ∼25 km) and are no longer suitable for small catchments or discipline scales. Optical/thermal infrared remote sensing records are known as the Surface Temperature/Vegetation Index Method and offer finer resolution (∼1 km). Recently, Zhang & Zhou (2016) presented a new approach by which SM estimation could be made from optical/thermal remote sensing, and this method especially relies upon the affiliation among SM and the surface reflectance and the temperature or vegetation index. Retrieval strategies in this area, like thermal inertia, emphasize soil thermal traits or triangular dating techniques, indicating that the linkages among SM, the Normalized Difference Vegetation Index (NDVI), and the Land Surface Temperature (LST) of a given area are being utilized in different applications. However, their applications are limited due to the lack of sufficient spatial data including topography or low-density vegetation cover maps and data. The remotely sensed vegetation indices to estimate SM (e.g., NDVI, the Normalized Difference Water Index (NDWI), and the Normalized Multi-Band Drought Index (NMDI)) are proper alternatives; however, the distribution of SM cannot be anticipated by a single parameter and by figuring out parameter modifications among specific land-surface-aspect intensities. Extensive efforts have been made to estimate SM by satellite images by establishing a connection between remotely sensed LST and vegetation indices (e.g., Dari et al. 2021; Zhu et al. 2021). One of the practical advantages of remotely sensed imagery is that vegetation and LST parameters are available with a high spatial resolution (from 30 m to 1 km) via imagery, in addition to topography data. Predicting the SM condition by utilizing the structured landscape factors extracted from remotely sensed imagery instead of in situ measurement makes it viable to obtain quick and real-time tracking of the SM condition.

As stated above, because of the complicated nature of the SM distribution process, a wide variety of ambiguous land covers (especially vegetation cover) and the effect of different functions, which include topography and terrain, make SM studies more complicated. The growing availability of remotely sensed imagery, linked to advances in modern computers, effectively contributes to low-priced broad-scale monitoring of hydro-environmental processes. In this regard, black box (lumped) methods may be superior to physical-based methods for a spatial-temporal estimation of SM. Unlike the data-driven methods, physical-based models require sequential weather elements and parameters of soil and vegetation types. Therefore, physical-based models are not so flexible for real-world practices because of the requirements of spatial parameters and high-quality field data. Neural networks is a kind of data-driven approach that describes nonlinear and complicated relationships to map inputs into a desirable output based on a sequence of complicated mathematical functions without any mechanistic information of the system. Several research studies have already described the application of neural network–based models to estimate SM over watersheds (Jiang & Cotton 2004; Pandey et al. 2010; Kornelsen & Coulibaly 2014).

Moreover, to handle massive data and records provided by remotely sensed satellites (for every pixel, there are numerous numbers or indices) that are to be incorporated into the model that uses remotely sensed land surface parameter data (e.g., SM), some data preprocessing algorithms are required. Clustering is one such statistical preprocessing technique of classifying data to homogeneous groups to make an operational assessment of the records through a regionalization of the study area. A cluster analysis can be employed for classifying multivariate data into subgroups to enhance the overall performance of the modeling. Clustering strategies recognize the shape of unlabeled data sets by dividing the members with minimum within-group-item dissimilarity and maximum between-group-item dissimilarity into homogeneous groups. There are several methods to cluster the data. A self-organizing map (SOM) is a neural-based model permitting unsupervised classification based on similarity among separable data sets. Therefore, the SOM approach for clustering an evaluation of the complicated spatiotemporal data sets and dominant patterns has gained the attention of hydrologists. For instance, Hsu & Li (2010) employed the SOM to recognize the most homogenous patterns of precipitation stations at the high-dimensional wavelet-transformed space. For an evaluation of the drought frequency, Chen et al. (2011) used the SOM-based clustering method by employing the standardized precipitation index (SPI) and the geographic characters of the gauges to group the rain gauges into different clusters. For groundwater data clustering, Nourani et al. (2015) used the SOM to recognize the patterns of groundwater level to be fed into a feed-forward neural network (FFNN) for one- and multi-step-ahead groundwater level forecasting. The evaluation, simulation, and analysis of the SOM in order to solve hydrologic problems have been reviewed by Kalteh et al. (2008) for the fields of meteorology and oceanography, respectively. In the field of ecohydrology, Olden et al. (2012) used the SOM to regionalize the ecohydrological function and provided an explanation, because the final result given by the SOM ought to replicate the ecohydrological traits of regionalization. Due to the aforementioned advantages, the SOM as a robust method to recognize the spatial patterns of remotely sensed SM conditions, was used in the present study.

On the other hand, for the subsequent section of ANN-based modeling, the FFNN model with specific properties for nonlinear modeling was used in this study to predict the SM condition in the future. The capacity to recognize existing patterns from a given data set makes the FFNN to ferret out complicated relationships among the inputs and outputs. The conjunction of the FFNN with the mathematical tool of the clustering technique can be applied for the prediction of physical processes.

The main objective of the current study was to present a method to estimate the SM condition over the land surface. To this end, the time series of Moderate Resolution Imaging Spectroradiometer (MODIS) soil and vegetation indices and the digital elevation model (DEM) were used to obtain information about the SM condition over the study area. The method employs remotely sensed imagery to assist SM mapping, particularly for in situ data-poor areas. In order to fulfill this aim, two approaches in the field of artificial intelligence (AI) were applied. First, the SOM was used as a clustering technique to recognize homogeneous and comparable areas and essential pixels, after which for every cluster, a three-layer FFNN model for each centroid pixel of each cluster was designed as the method to discern the underlying nonlinear pattern within the data, spatially and temporally. Based on the literature survey and to the best of the author's knowledge, this is the first application of SOM-based clustering analysis to improve AI-based modeling of the SM to enjoy the advantages of both methods within a robust unique framework.

In the subsequent sections of the paper, remotely sensed data sources and their processing and the study area are presented. Thereafter, the proposed and used methods (i.e., SOM and FFNN) are reviewed and the technique for model development and evaluation criteria are explained. The outcomes acquired by the proposed technique are then presented and discussed, followed by the conclusions.

Located in the southwestern region of the USA (situated between latitudes of 32 °30′ and 34 °01′N and longitudes of 110 °18′ to 113 °20′W), in the south-central part of Arizona, Phoenix is a semi-arid desert region (Figure 1). It is situated at a mean elevation of 331 m, in the northeastern part of the Sonoran Desert. The topography of Phoenix is largely flat, except for highlands located around the city, including scattered and low mountains surrounding the valley: in the northeast, McDowell Mountain; in the west, White Tank Mountain; in the south/southwest, South Mountain and Sierra Estrella; and in the west, Camelback Mountain. North Mountain, Sunnyslope Mountain, and Piestewa Peak are located at the midpoint of the valley. The borders of Phoenix consist of large fields of irrigated cropland and Native American reservation lands. The Salt River flows westward through the city, but owing to massive diversions for irrigation purposes, the riverbed often goes dry.

Typical of the Sonoran Desert, the Phoenix climate (Koppen climate classification BWh) consists of subtropical desert. In the Phoenix metropolitan area and the surrounding agricultural land (in the center of Phoenix), the mean annual precipitation and evapotranspiration potentials are 180 and 2,286 mm, respectively. The summers of Phoenix are long and extremely hot and the winters are short and mild to warm. In USA, the average high temperatures in summer are the hottest among any main city. From late May to early October, there are 107 days in a year, on average, and on most of the days, the temperature is at least 38 °C. But even in peak summer when there is extreme heat, the highest average annual temperature is not reached.

In Sonoran, a hot desert, the average annual precipitation is 76–500 mm, which is dependent on the location, with considerable inter- and intra-annual changeability in timing and quantity. The air temperature in summer exceeds 40 °C and frequently reaches 48 °C. Severe thunderstorms of the summer monsoons occur due to the interaction of this high near-surface temperature with cool, moist air in the atmosphere. Subsequently, moisture on the soil surface and near-surface air evaporates and temperatures may drop by 10 °C or more, and this drop occurs following storms within a matter of minutes. In winter, the temperature is mild, with valley bottoms typically without frost, while the surrounding mountains may have compressed snow cover at high elevations and in the northern and eastern parts. During any season, daytime temperature swings of 15 °C or more are common, because the dry atmosphere and relatively low vegetation cover simplify the process of re-radiation of daylight heat into the atmosphere overnight.

R stands for the reflectivity of each band. In MODIS sensors, 860 nm represents the Channel 2 band, 1640 nm represents the Channel 6 band, and 2130 nm represents the Channel 7 band (Wang & Qu 2007). This is an index exhibiting good sensitivity to drought severity, and, therefore, Zhang et al. (2009) used this index for Henan drought supervision and ascertain whether the NMDI can be used in areas with modest vegetation coverage. Mobasheri & Bidkhan (2013) studied the relationship between two SM indices: NMDI and PDI, and SM content (%) on the basis of in situ spectral measurements. They showed that the mentioned indices are well suited for monitoring SM, particularly for bare soil or thinly vegetated areas.

A robust absorption of near-infrared (860 nm) light in the vegetation canopy is used in this index, as well as the sensitivity of short-wave infrared 1,640 and 2,130 nm bands for the absorption of variance between soil and vegetation water contents. The SM condition includes three classes (on the basis of volumetric SM range): dry: 0–0.1; intermediate: 0.1–0.2 and wet: >0.2 (Miller et al. 2004); consequently, a pixel will be mapped as a dry soil condition, related to bare soil or areas with weak vegetation, if the NMDI is <0.7, intermediate condition if the NMDI is in the range of 0.6–0.7, and wet if the NMDI is <0.6. The DEM extracted from ASTER, a sensor mounted on Terra satellite, as part of NASA's EOS (Earth Observing System), was used in the current study with a 30 -m resolution retrieved via https://earthexplorer.usgs.gov/.

The raw data collected from MODIS images are in the form of sinusoidal projection and could be transferred to the UTM coordinate system via the MODIS Conversion Toolkit (MCTK). The images from the MODIS surface reflectance and DEM have 500- and 30-m spatial resolutions, respectively; hence, they were resampled to match the spatial degree of the MODIS LAI and LST products with a 1- km spatial resolution. Using vector data related to the study region boundaries, a water mask was applied, and finally the clouds were removed from the images. Each of the satellite images was cropped to the geographic extent of the study area. For each of MODIS product (e.g., LST, LAI, and NMDI) in order to create 8-day time series matrix for the study period (2001–2004), a three-dimensional matrix 184×(168 × 193) was formed, including 184 records (with 8-day interval) for 168 × 193 spatial cells. After fully removing unvalued data, the size of the input matrix for SOM classification included 553 (184 LST images+184 LAI image+184 NMDI image+1 DEM) rows and 21,022 (168×193 unvalued data) columns.

The 8-day MODIS products (LST, LAI, and NMDI) of the Phoenix for 2001–2004 were used for FFNN modeling of the centroid pixels of clusters: the first 3 years for the training phase and the rest for the verification phase. Extreme values of the inputs were separated for the training phase to make the learning process more reliable.

It is necessary to combine appropriate information about SM conditions in different pixels of the region and recognize the forceful pixels to predict SM in future. Therefore, the SOM as an ANN-based clustering approach, a common machine learning method for pattern recognition, is able to discover the structure of relationships among groups of unlabeled data without considering the physically inherent natural phenomenon. One simple description of the SOM in this study is that it is able to separate the images into meaningful regions with homogeneous behavior. The proposed model consists of two phases. After preprocessing satellite data, in phase 1, a two-step SOM was employed for classifying pixels into classes with close patterns. Generally, such a two-step SOM clustering is proposed to gain an understanding of homogeneous regions and estimate the number of clusters with regard to the topology of the plain. Then, for determining the predominant (center) pixel of each cluster, the Euclidean distance criterion was used. The predominant pixel was determined in order to extract the best sample to represent the pattern of the whole cluster.

The SOM is one of the high-performance tools for analyzing data with high dimensions. It orderly maps a distribution with high dimensions to a grid with regular low dimensions. Therefore, it can convert statistically complex and nonlinear relationships of data with high dimensions to simple geometric relationships with low dimensions. Moreover, the SOM is proven to effectively preserve the topological structure of the data during the conversion process (Kohonen 1990). The application of the SOM network in image clustering is to place the pixels with similar attributes in the same cluster by reducing the dimension and creating a map of commonly one or two dimensions. Therefore, the SOM achieves two goals of distinguishing similarities and dissimilarities and reducing the dimension of the data. Commonly, SOM architectures are double-layered, consisting of the input layer and the Kohonen layer. The described two-level SOM neural network (see Figure 2) is shown to be an efficient method to obtain an initial understanding of the complicated data sets. The classical SOM network can be extended by another one-dimensional Kohonen layer, where neurons are interconnected to the neurons of the previous Kohonen layer.

Steps of iterative learning of an SOM network are described as follows:

• The weights of each node are randomly assigned for initializing the network.

• A vector from the input vector (x) is selected randomly.

• All nodes are inspected to find the weights most similar to the input vector. The Euclidean distance is used as distance measurement (similarity and dissimilarity).

• The Best Matching Unit (BMU) is the winning node with the closest weights to the related input. The BMU and the nodes in the neighborhood are calculated and permitted to learn by changing weights through iterations. The iterations are processed to obtain the closest distance between the input and the weights.

• The training process continues until obtaining the convergence results.

• Finally, the clusters (homogeneous region) will be recognized on the basis of the structure of the SOM.

The Levenberg–Marquardt was selected as an optimization algorithm for a three-layer FF–BP network model in this study. It is worth mentioning that the efficiency of the developed FFNN-based models is associated with the appropriate selection of the neuron number in the input and hidden layers and the number of training iterations (epoch). Additionally, the model may display overfitting in the training phase due to a high epoch number and poor quality or quantity of data, thus limiting the model's performance by limiting the ability to predict new data outside the training set.

In Equations (3) and (4), M denotes the number of data, is the observed target value of the data and is the calculated value for intercepting the observed target, and is the average value of the observed data. In a high-performance model, the DC is closer to one, and the RMSE is closer to zero.

The satellite-based data (as introduced in Table 1) were used for SOM-based clustering of the area. A total of 553 normalized preprocessed satellite images in the form of (193 × 168) matrixes with a pixel resolution of 1,000 m were considered as inputs of the SOM in order to classify the region of the case study. Studies conducted on the NMDI illustrate that it does not respond to moisture variation of the soil, starting from LAI equal to 2, which means no soil background impact is found on the NMDI for any SM domain (Wang & Qu 2007). Therefore, for a vegetation canopy with an LAI equal or greater than 2, which means heavily vegetated areas, the NMDI becomes a complete index for predicting vegetation water content, compared with the SM index. In order to achieve the overall view of the situation of SM over the study region, pixels with LAI values higher than 2 were removed from the data.

Since selecting the appropriate number of neurons is based on a topology of inputs in the first step to form the two-step SOM, the size of the Kohonen layer could be determined as a 10-by-10 lattice by a trial and error process. Figure 4(a) illustrates the hits plan of the output layer for the applied data. The SOM output layer is explained by the hits plan in which every single neuron represents the number of arranged input vectors. Also, the size of the colored patch reveals the proportionate number of vectors. Figure 4(b) presents the neighbor weight distances in which the darker colors represent large distances and the lighter colors exhibit small distances. After forming a 10-by-10 Kohonen layer, the associated zones of the clusters on the map could be recognized on the basis of the number of zones allotted to the neurons (see Figure 4(a)). A one-dimensional SOM with neurons arranged within a 1-by-5 Kohonen layer was then organized in the second step of the clustering to group the pixels into five distinct classes, in which the number 5 for the final classes of numbers would be determined in the first step of clustering (two-dimensional map), visually.

Since SM is known as a dynamic component of the underlying process happening on the surface of the earth, it is important to consider the effects of dynamic temporal variability in the classification parameters. In order to investigate the impact of time series as the input of SOM classification parameters for the study area, the time series and the average of time series for 4 years were imposed to the SOM separately. The silhouette coefficient as a cluster validity measure was utilized to assess the homogeneity of the clusters obtained by the proposed two-step SOM.

Figure 5 illustrates the performance of each SOM clustering for two types of selected input data. Horizontal and vertical axes show the number of clusters and silhouette value, respectively. Although an averaging of time series produces higher silhouette values for all numbers of clusters, and the average of the data reduces the complexity and the costs of calculation, averaging can conceal the temporal fluctuation within the data. Therefore, the whole time series imposed to the SOM was preferred for the classification. Since the maximum silhouette value was obtained when using clusters 4 and 5, these are the eligible cluster numbers and a strong-enough reason to choose 5 as the best cluster number in this study.

It is clear (as seen in Figure 5) that by increasing the number of clusters, the Silhouette coefficient decreases, which denotes a higher homogeneity of the clusters and efficiency of the clustering. On the other hand, more clusters mean more FFNN models are needed to model the SM over the region. It is seen in Figure 5 that there is a significant decrease in the Silhouette coefficient till the cluster number 5, and thereafter, by increasing the cluster number, only a lower homogeneity can be seen in the clusters. Therefore, and to have a tradeoff between the modeling cost and accuracy, the clustering size of 5 was considered in this study to classify the region.

The mentioned SOM clustering technique applied to the set of remotely sensed data clustered the region into five clusters. Figure 6 shows a spatial distribution of the clustered pixels and the SM condition in the study region with regard to the input data used for the clustering (LST, DEM, and LAI). Considering the DEM map of the area (Figure 6(b)), the diversity of the clusters in the higher elevations located in the northwest part of the map indicates that the rate of changes in parameter comparison with other areas is remarkable. Also, it is clear by considering Figures 6(a)–(d) that the spatial distribution of clusters is correlated more to the mean LAI and LST data than to the DEM, where such data include temporal variations and can handle the seasonality of the process, but the DEM does not have temporal variations. To explore the capability of the SOM, the main statistical features of each attribute are described and discussed in the following paragraphs to characterize the spatiotemporal variability of the process.

Figure 7 shows the distribution of elevation and average time series of the LST for each cluster sorted on the basis of increasing the mean value of each cluster. As expected, the elevation and LST have a negative relationship, which means clustering with high elevation typically produces lower average temperatures. Clusters 1 and 2 have a wide variation range of both parameters among other clusters. The downstream areas with low elevation were allocated to clusters 4 and 5, and gradually with increasing the elevation, the pixels located in other clusters, which illustrates the critical role of elevation in the clustering outcome (see Figure 7(a)).

LST histogram presented by Figure 8 shows the high non-homogeneity of LST values in the clusters. The histograms of the NMDI for the clustered areas demonstrate a regular distribution of this variable for each cluster (Figure 9). Most of the percentages of the NMDI values of the pixels of clusters 1, 2, and 3 are lower than 0.6 (based on the NMDI description in the data set section), so these clustered areas are mapped having a wet soil condition, and cluster 4 has a normal soil condition; cluster 5, which tends to the right of the diagram, denotes the dominant nature of the dry soil condition.

According to the range of LAI variation in Figure 10, clusters 1 and 2 are included in a wide range of LAI values, so this indicates that the non-homogeneity of these clusters is remarkable. But the remaining clusters have a low range of LAI values.

As is illustrated in Figure 7, clusters 4 and 5 are located in low elevation lands and high-temperature areas, which cover the Sonoran Desert. During dry conditions in these clusters, a low SM (see Figure 8(d) and (e)) can also strongly limit evapotranspiration, and, in turn, can reduce the potential variability caused by vegetation covers (see Figure 10(d) and (e)).

The downhill of the northeastern mountains is allocated to cluster 3. The high evapotranspiration rates in summer, as well as poor vegetation cover (Figure 10(c)), reduce the SM in this area. But the impact of topography on water redistribution, which is caused by rainfall at high altitudes, will increase the SM during rainy days and improve vegetation cover in wet conditions.

Clusters 1 and 2, which include the mountainous region and agricultural land surrounding the Phoenix metropolitan area, illustrates the wet soil condition. During the growth season in these areas, decreasing temperatures and maximum rainfall could improve the vegetation cover, and consequently, the evapotranspiration rate of the SM declines.

In the last step of the modeling (prediction step), an FFNN was trained for each of the representative pixels of each cluster to predict the NMDI one time-step ahead. The LAI and LST data sets of the present step and NMDI values of the previous step were used as the inputs for the training phase. These input orderlies, especially the NMDI lags, were selected by sensitivity analysis to consider the seasonality and autoregressive characteristics of the process. The results showed that the model efficiency was not strongly dependent on cases where larger lags were considered. Therefore, to obtain satisfying results, the minimum numbers of the input data were used to reduce the error of imposing redundant information and noise. The verification data were used to test the efficiency of the trained models. The obtained results are presented in Table 2. According to these results, it is clear that by imposing the NMDI value at the previous time step (NMDI(t − 1))) as input, the performance of the model improves, which indicates the autoregressive property of the process.

When each of the LAI or LST parameters was individually imposed (as the sole input) to the FFNNs, the modeling performance (in terms of verification DC) dipped, and the DC was less than 0.4 in.all clusters, which indicates poor modeling by these input sets. But when both LAI and LST were used as inputs of FFNNs simultaneously, the verification DC was at least 0.6. Adding the NMDI (t − 1) to the input layer resulted in a DC verification of at least 0.76, which shows the autoregressive property of the process.

It should be noted that the temporal modeling size of the study was 8 days, but when the modeling lag (horizon) was increased to have two- or three-time-step-ahead modeling (as a multi-step-ahead prediction scheme), the verification DC dropped below 0.5 in most cases, and because of this, the results of the multi-step-ahead modeling are not provided in the table. The reason for this poor multi-step-ahead modeling performance is that, after 8 days, climatic and soil conditions can undergo large changes, which make it more difficult to find a mathematical relationship between the present conditions and those that prevailed 16 (or 24) days back.

In the present research, a combined SOM clustering method and FFNN approach was applied to devise an AI-based black box methodology for the multivariate SM condition modeling of the study region located in the southwestern part of the USA. In the first step, the SOM was utilized to recognize the pixel data with high similarity and to determine the dominant pixel in each cluster to identify the best patterns as a cluster indicator. The selection of suitable pixels by the SOM resulted in the definition of regions with similar SM conditions. Also, it could be concluded that sufficient in situ SM measurements can be done in the future. By selecting such representative and dominant pixels, more tools and investments on these dominant points for in situ measurement may be employed to improve the data-recording performance. In this study, nonlinear relationship between LAI, LST and previous (lagged) values of NMDI were discovered by the FFNN for each of the representative pixel of each cluster to predict the SM condition. The obtained results indicated encouraging evidence when comparing the observed and predicted NMDI values.

As with any kind of black box (data-driven) method, the numerical results obtained through the current modeling are case sensitive for the study area data. But this methodology can be applied similarly to any other area, and it is expected that the same regression will be seen between the same data but with different order of strengths. As a suggestion for future studies, the proposed methodology can be further developed using other hydrological time series and variables as additional inputs for predicting the NMDI (e.g., precipitation or/and evapotranspiration). The combination of the SOM and FFNN models can be improved and updated by more recent remotely sensed images and field analysis. Due to the fact that the presented method of the current research mostly relies on data, which is a limitation of this study, using a high volume of data may increase the noise and produce abundant information, resulting in a time-consuming training phase with inefficient results. To overcome this limitation, an SOM-based clustering approach that employs wavelet-based entropies of time series as inputs of clustering is suggested to optimize the SOM input layer (see Nourani et al. 2015).

Data cannot be made publicly available; readers should contact the corresponding author for details.