Surface and Subsurface Soil Moisture Estimation Using Fusion of SMAP, NLDAS-2, and SOLUS100 Data with Deep Learning

1. Introduction

Soil Moisture (SM) is a key state variable that controls the exchange of water, energy, and carbon fluxes between the land surface and the atmosphere [1,2]. SM has numerous applications in hydrology, agriculture, and meteorology, such as irrigation management, drought monitoring [3], flood prediction [4], and weather forecasting [5]. However, accurate measurement and monitoring of SM in different spatial and temporal scales is still a challenge. While in situ SM networks such as International Soil Moisture Network (ISMN) provides more than 2,800 active stations, their spatial coverage is insufficient due to uneven distribution of the stations which limits their ability to represent the spatial variability of SM caused by heterogeneity of land surface factors, especially those related to soil and vegetation types [6].

During the past two decades, optical and thermal satellite remote sensing observations have been widely used for SM estimation due to their higher spatial resolution compared to microwave satellite observations [7,8,9,10]. These methods typically employ triangle or trapezoid models that use parameters such as land surface temperature – vegetation index [7,11] or shortwave infrared reflectance – vegetation index [12,13] to estimate SM. While these methods can provide SM estimates with high spatial resolution (20 to 500 m pixel size), they have limitations such as low penetration through vegetation canopies and reduced or irregular temporal resolution due to unstable atmospheric conditions (e.g., cloud cover, aerosols, etc.) [14]. In contrast, passive and active microwave satellite observations are unaffected by atmospheric conditions and their longer wavelengths enable deeper penetration through vegetation and soil surface. Microwave signals are highly sensitive to the dielectric properties of soil, which are primarily influenced by large dielectric constant of water [15,16]. NASA’s Soil Moisture Active Passive (SMAP) satellite was launched in 2015 to specifically measure global SM at L-band (1.4 GHz) within the top 5 cm of soil under various land covers with a spatial resolution of 9-36 km and a temporal resolution of 2-3 days using both descending (6:00 AM) and ascending (6:00 PM) overpasses [17]. SMAP passive observations in the form of brightness temperature are typically combined with other soil and vegetation attributes including scattering albedo, surface roughness, vegetation optical depth (VOD), and soil dielectric properties to estimate SM using the inversion of the tau-omega radiative transfer equation model, knows as zero-order model [18,19]. The tau-omega model is solved with the single channel algorithm (SCA) or the dual channel algorithm (DCA). In the SCA, brightness temperature observations at either horizontal (H) or vertical (V) polarization are used to retrieve SM, while the DCA uses both H and V polarizations to retrieve SM and VOD simultaneously [20,21]. The DCA generally yields more accurate results than the SCA, as it uses both polarizations and is currently the baseline algorithm to produce official SMAP SM products [22,23]. However, some studies have shown that the SCA-V algorithm may perform better under specific land covers [24,25]. Although SMAP primary product (Level 3) provides surface SM (top 5 cm) with a target unbiased root mean squared error of 0.04 cm3 cm-3 over areas with low vegetation water content (less than 0.5 kg m-2), it does not provide daily estimation and faces challenges in accurately retrieving SM in regions with dense vegetation cover, complex topography, and snow-covered soils [26]. Moreover, a limitation of the tau-omega model is that it does not account for temporal dynamics of vegetation and growth stages which may affect the accuracy of SM retrievals. The SMAP Level 4 SM product provides an average SM estimate from the surface down to 100 cm depth with a spatial resolution of 9 km and temporal resolution of 3-hours. This product is generated by assimilating SMAP brightness temperature observations into the NASA catchment land surface model [27]. Although this SM product has very high temporal resolution, its vertical resolution (within the soil profile) is course that limits its usefulness for applications that require surface, near-surface or root zone soil moisture information.

Alternative databases such as the North American Land Data Assimilation System phase-2 (NLDAS-2), which is the output of process-based models, offers SM estimates with high temporal resolution (hourly) and 0.125° (~12.5 km) spatial resolution for multiple soil layers [28]. However, the accuracy of this product is influenced by regional factors such as topographic complexity (TC), VOD, and uncertainties in model parameterization caused by the use of simplified representations of physical processes, such as one-dimensional vertical water flow in soil [29], which often lead to underestimation of SM during wet seasons and overestimation during dry conditions [30]. Despite such limitations, several studies have shown that NLDAS-2 SM outperforms SMAP SM estimates under forested and unforested land covers [31]. NLDAS-2 high temporal resolution and its ability to provide SM estimates for multiple soil layers give it an advantage, especially for root zone and profile SM monitoring where SMAP satellite observations, which primarily focus on the surface layer, are limited.

In recent years, machine learning and deep learning (DL) techniques have shown superior performance compared to traditional statistical and physical models [32]. These methods can be considered as an alternative to the direct SCA and DCA algorithms for SM estimation at large scale (33-38). The main advantage of DL methods includes their strong learning ability to find highly non-linear spatiotemporal relationships and capability to handle multi-source multi-scale observations [39]. In contrast, radiative transfer (e.g., SCA and DCA) and land surface models struggle to capture the inherent heterogeneity in soil and vegetation and the non-linear relationships between SM and land surface parameters [40]. The NLDAS-2 SM data can be used as the benchmark to train DL models for SM estimation as DL uses many neurons and hidden layers to identify complex patterns in the data, making them a highly promising approach for SM estimation [41,42]. However, much of the existing research on DL for SM estimation has focused on small scale modeling, which limits the ability to fully capture the intricate spatial and temporal non-linear relationships. Recently, studies presented a deep neural network that fuses SMAP ancillary variables with reanalysis SM data from ERA5 to produce a high accurate SMAP SM product for a wide range of land-covers and climate regimes [39]. The hybrid DL models such as Convolutional Neural Network – Long Short-Term Memory (ConvLSTM) has shown significantly better performance in this regard [43]. ConvLSTM combines the strengths of convolutional neural networks [44] for spatial feature extraction and long short-term memory (LSTM) networks [45] for capturing temporal dependencies. This enables ConvLSTM to concurrently model spatial and temporal variations in SM more effectively. As a result, ConvLSTM has gained attraction for estimating and forecasting of various hydrological variables, such as precipitation [46,47], evapotranspiration [48,49], flood [50], and streamflow [51].

Incorporating soil information (e.g., soil physical properties like texture) into machine learning and DL models significantly enhances SM estimation by enabling the models to capture complex relationship between soil characteristics and SM [52]. Soil texture, which affects the pore size distribution and water-holding capacity of soil, plays an important role in determining retention and movement of water within the soil profile [53]. By integrating soil physical properties like texture maps, especially the new released NRCS digital soil maps, as additional features, ConvLSTM can identify spatiotemporal patterns and correlations that contribute to more accurate estimates of SM. Moreover, including spatially distributed soil physical properties facilitates a better understanding of variability across different landscapes, improving model generalization and performance under diverse soil types, land covers, and climatic conditions.

This study aims to employe the ConvLSTM model to generate daily SM products for surface and subsurface soil by integrating multi-source multi-scale ground and satellite observations across the contiguous United States (CONUS). The specific objectives of this paper are: (i) to apply ConvLSTM for creating new surface and subsurface SM products by integrating SMAP SM and ancillary data, NRCS soil physical properties maps, and NLDAS-2 SM data, and (ii) to evaluate the accuracy of the new SM products with in-situ SM networks for different soil textural classes and land covers.

2. Materials and Methods

The Materials and Methods should be described with sufficient details to allow others to replicate and build on the published results. Please note that the publication of your manuscript implicates that you must make all materials, data, computer code, and protocols associated with the publication available to readers. Please disclose at the submission stage any restrictions on the availability of materials or information. New methods and protocols should be described in detail while well-established methods can be briefly described and appropriately cited.

2.1. SMAP Data and SOLUS100 Maps

The NASA’s Soil Moisture Active Passive (SMAP) mission was launched in 2015 to provide global SM estimates within the upper 5 cm soil. SMAP includes several SM products such as the Level 3 (L3) enhanced SM product with 9 km spatial and 2-3 days temporal resolution [19,54]. Here, we used SMAP L3 ancillary data including brightness temperature, albedo, effective SM, surface roughness, vegetation opacity, and vegetation water content as the ConvLSTM model inputs. These variables are currently used in the SMAP SM algorithms (e.g., DCA and SCA) to estimate surface SM. The SMAP L3 data (Version 6) were acquired from the National Snow and Ice Data Center (NSIDC) at EASE-Grid 2.0 Global projection (EPSG:6933) from January 2018 to December 2022.

Additionally, we used soil physical properties including texture (i.e., sand, silt and clay content) and bulk density from the USDA-NRCS Soil Landscapes of the United States (SOLUS100) maps. SOLUS100 is a new digital soil map generated from the fusion of soil profile data with field descriptions and soil survey maps via random forest modeling to provide high accurate physical and chemical soil properties maps at 100 m spatial resolution at seven standard soil depths (0, 5, 15, 30, 60, 100, and 150 cm) across the CONUS [55]. In this study, we used SOLUS100 maps for sand, silt, clay and bulk density properties for two depths, where the soil maps at 5 cm depth correspond to NLDAS-2 SM at 0-10 cm, and the maps at 30 cm depth correspond with NLDAS-2 SM at 10-40 cm.

2.2. Producing Benchmark Soil Moisture Data

To develop ConvLSTM models, we integrated two distinct SM products: SMAP SM (satellite product) and NLDAS-2 SM (process-based model) to leverage the strengths of both as pointed out above. Specifically, the SMAP Level 3 enhanced soil moisture product (SPL3SMP_E), which includes three SM products based on DCA, SCA-V, and SCA-H, and the NLDAS-2 Noah SM product at 0-10 cm and 10-40 cm depths [56], were utilized from January 1st, 2018, to December 30th, 2022. As mentioned above, there is no considerable study that compare the performance of SMAP SM products with other SM products such as NLDAS-2, as the performance is variable depending on land cover conditions [57]. Based on our validation results (Table 1), SMAP DCA SM demonstrated better correlations with ins-situ SM from both SCAN and USCRN networks, while the NLDAS-2 Noah SM showed better RMSE values across the CONUS. Therefore, the ConvLSTM model was designed to learn and produce new SM from the high correlation of SMAP DCA while consider low RMSEs from the NLDAS-2 Noah SM (depth 0-10 cm) for the surface depth (i.e., 5 cm). For the subsurface (i.e., 25 cm), the model relied solely on NLDAS-2 subsurface SM (10-40 cm) as SMAP Level 3 lacks SM data for subsurface soil. Recently, it has been demonstrated the applicability of integrating SMAP surface SM with reanalysis SM data from ERA5 to produce benchmark SM data for improving the training of deep neural networks over various land cover types across the CONUS [39].

To do so, the soil maps and NLDAS-2 data were resampled to a 9 km pixel size using bi-cubic (dis)aggregation method [58] to obtain a uniform spatial resolution based on SMAP L3 data. To preserve temporal continuity (i.e., 1 day), the SMAP data were gap-filled temporally using the Continuous Recurrent Unites (CRU) method prior to integration and subsequent processing using the ConvLSTM algorithm. The CRU is an encoder-decoder neural architecture designed to gaps in time series data by capturing temporal continuity between hidden states and optimally integrating noisy observations [59]. Figure 1 depicts the steps of data integration and deep learning modeling for surface and subsurface SM estimation with ConvLSTM model.

2.3. Model Input Scenarios

To estimate surface a subsurface SM, two scenarios were defined for each depth. While the SMAP ancillary data are used for both surface and subsurface and the two scenarios, for the surface (5 cm), the scenarios are based on whether soil physical attributes are included in the predictors or not. For the subsurface (25 cm), the scenarios are defined based on whether the generated surface SM is included as a predictor. These scenarios enable us to assess the effect and contribution of SMAP ancillary data and soil properties to estimating surface and subsurface SM. Table 2 summarizes the used scenarios and predictors for each soil depth.

2.4. ConvLSTM Model Architecture and Development

To estimate SM, we integrated convolutional neural network with long short-term memory network, named as ConvLSTM, to identify both spatial and temporal patterns in SM dynamics and high nonlinear relationships between SM and the predictors. The ConvLSTM is mathematically defined as follows [43].

$$i_t=\sigma (W_{xi}\mathrm{*}{x}_t+W_{hi}\mathrm{*}{h}_{t-1}+W_{ci}\circ {C_{t-1}+b}_i) (\mathrm{Input} \mathrm{Gate})$$

(1)

$$f_t=\sigma (W_{xf}\mathrm{*}{x}_t+W_{hf}\mathrm{*}{h}_{t-1}+W_{cf}\circ {C_{t-1}+b}_f) (\mathrm{Forget} \mathrm{Gate})$$

(2)

$$o_t=\sigma (W_{xo}\mathrm{*}{x}_t+W_{ho}\mathrm{*}{h}_{t-1}+W_{co}\circ {C_t+b}_o) (\mathrm{Output} \mathrm{Gate})$$

(3)

$$\widetilde{C_t}=tanh(W_{xc}\mathrm{*}{x}_t+W_{hc}\mathrm{*}{h}_{t-1}+b_c) (\mathrm{Candidate} \mathrm{Cell} \mathrm{State})$$

(4)

$$C_t=f_t\circ C_{t-1}+i_t\circ \widetilde{C_t} (\mathrm{Cell} \mathrm{State} \mathrm{Update})$$

(5)

$$H_t=o_t\circ tanh\left(C_t\right) (\mathrm{Hidden} \mathrm{State})$$

(6)

$$Y_t=W_y\mathrm{*}{h}_t+b_y (\mathrm{Output} \mathrm{Layer})$$

(7)

where $x_t$is the current input at time step $t$, $h_{t-1}$ is the previous hidden state (short-term memory),$C_{t-1}$ is the previous cell state (long-term memory), $f_t$, $i_t$, and $o_t$ are forget, input and output gates, respectively, $\widetilde{C_t}$ is the candidate cell sate, W are the network weights, $b$ are the biases, $\sigma$ is the sigmoid activation function, $tanh$ is the hyperbolic tangent activation function, $\mathrm{*}$ denotes the convolution operation, and $\circ$represents element-wise multiplication. At each time step $t$, the network outputs a value ($Y_t$) which is compared to the reference soil moisture data to compute the loss and update the network weights during training.

As shown in Figure 1, after unifying data (i.e., gap-filling, resampling, and reprojecting), each predictor was reshaped into thousands of smaller images to expedite the execution of the ConvLSTM and enhance the model’s capability to capture parameter variations [43]. Subsequently, the predictors were organized into sequences of three consecutive days (lookback: 3) to better explore the temporal relationships [60]. The dataset was split into training (years of 2018 to 2020), validation (year 2021), and testing (year 2022) sets which approximately correspond to 60% training, 20% validation, and 20% testing.

The developed ConvLSTM model consists of ten layers with varying number of filters (i.e., 16, 32, 48, and 64) in each layer, capturing progressively more complex spatiotemporal patterns as data passes through the layers [43]. The kernel size was set as 3, the activation function as ‘sigmoid’, padding as ‘same’, and optimizer as ‘Adam’. A dropout rate of 0.2 was used after each layer along with early stopping to prevent overfitting [61]. Batch normalization was applied after each layer to stabilize and accelerate the training process. Moreover, the model was equipped with two Attention Blocks after the second and fifth layers that enables the model to focus on important features in both spatial and temporal dimensions [62]. Each attention block performs global average pooling operation to summarize features across spatial dimensions, followed by two dense layers that first reduce the number of features and then restore them to the original size, generating attention weights. These weights are then multiplied by the original inputs to assign higher importance to more important features and enable the model to focus on key spatial patterns [64]. As mentioned above, to develop ConvLSTM model, we used the SMAP SM and NLDAS-2 SM products as benchmark data. For the surface depth, this was based on the NLDAS-2’s low RMSE values and the SMAP DCA SM’s high correlation with in-situ SM measurements. Hence, we defined the loss function to minimize both the RMSE and (1-R²) between the model outputs and SMAP and NLDAS-2 SM products, respectively, as formulated in the following equation.

$${\mathcal{L}}_{ConvLSTM}=min\left\{\alpha .RMSE+\beta .(1-R^2)\right\}$$

(8)

where $\alpha$ and $\beta$ are the weight parameters set to unity. For estimating subsurface SM, the loss function specifically focuses on minimizing RMSE values between the model outputs and NLDAS-2 SM, as there is no SMAP SM product available for the subsurface soil layer.

2.6. In Situ SM and Model Validation

The performance of the developed ConvLSTM models for surface and subsurface SM estimation was assessed using in-situ SM measurements from the Soil Climate Analysis Network (SCAN) [65] and the U.S. Climate Reference Network (USCRN) [66] for depths of 5-, 10-, and 20- cm for the testing set (year 2022). To further evaluate the model performance in different land-cover types and soil textures, we used MODIS satellite land-cover product (MCD12Q1) [67] and SOLUS100-based soil textural classes based on USDA classification scheme. The locations of in-situ stations along with the distribution of soil texture classes are displayed in Figure 2. The error metrics including root mean squared error (RMSE), mean bias error (Bias) and unbiased RMSE (ubRMSE) were used.

3. Results and Discussion

In this section, we first evaluate the performance of the ConvLSTM model for surface and subsurface soil moisture estimation based on two sets/scenarios of predictors and then examine the model performance for various soil textural classes and land cover types across the CONUS.

3.1. Surface and Subsurface SM Estimations with ConvLSTM

As explained above, two different scenarios including SMAP ancillary data with and without soil physical attributes from SOULS100 maps were considered for SM estimation in each soil depth. When comparing the overall performance of the ConvLSTM model for each scenario, we found that the second scenario, which includes SOLUS100 data, outperformed the first scenario (Table 3), revealing that soil physical properties including sand, silt, clay content and bulk density play an important role in spatiotemporal SM estimation even when the maps are aggregated to coarse pixels with 9-km size at large scale. The importance of soil physical properties in SM estimation and monitoring has been shown in other studies [52,68,69,70]. The mean R and ubRMSE values for the surface depth and for scenarios 1 and 2 were 0.54 and 0.053 cm³ cm^-3, indicating similar accuracy for the two scenarios. For the subsurface depth, scenario 2 slightly outperformed scenario 1 with mean R and ubRMSE values of 0.45 and 0.043 cm³ cm^-3 (scenario 1) and 0.51 and 0.041 cm³ cm^-3 (scenario 2) respectively. This could be due to including the estimated surface soil moisture as an additional input (see Table 2). The mean ubRMSE values meet the target accuracy for SMAP SM products (i.e., ubRMSE of 0.04 cm³ cm^-3) for vegetation water content less than 0.5 kg m^-2 [54] however the ConvLSTM SM estimates include advantages like daily SM for surface and subsurface soil depths. Though the ConvLSTM model was trained with both NLDAS-2 SM and SMAP L3 SM data, it still provides acceptable accuracy for estimating surface and subsurface SM.

3.2. Effects of Soil Texture on Soil Moisture Estimation

To better determine the effect of soil texture on the accuracy of SM estimation with the ConvLSTM model, we further classified the twelve soil textural classes (Figure 2) into three major classes: fine-textured (sandy clay, silty clay, and clay), medium-textured (loam, clay loam, silt loam, sandy clay loam, silty clay loam, silt) and coarse-textured (sand, loamy sand, and sandy loam) soils [70]. This grouping reduces the complexity of the soil texture classes for interpreting the results and helps identify general trends in SM estimation. The estimated SM was evaluated within each soil texture group against in-situ SM measurements from both SCAN and USCRN. The boxplots shown in Figure 3 display the error metrics of SM estimates for the surface and subsurface depths for the three groups of soil texture. As seen, based on the ubRMSE values, the accuracy of SM estimates for coarse-textured soils is slightly higher than that of medium-textured soils, followed by fine-textured soils. This could be partly attributed to the narrow range of SM content in coarse-textured soils due to their smaller porosity (i.e., saturation water content) and low water holding capacity. Similarly, the mean ubRMSE values for the surface depth are slightly lower than those of NLDAS-2 SM for fine- and coarse-textured soils, while no difference was observed for the subsurface SM estimates between the ConvLSTM and NLDAS-2. Based on the correlation coefficient (R) values, the model performed better for fine-textured soils, followed by medium-textured soils. The model performed slightly better for estimating subsurface SM compared to surface SM, which could be attributed to the lower dynamics of SM in deeper soil layers. By comparing scenarios 1 and 2, it is observed that the model performed slightly better in terms of ubRMSE for scenario 2 in medium- and fine-textured soils, while no significant improvement was observed between scenarios 1 and 2 for coarse-textured soils. Recently, SMAP level 3 SM products have been tested for estimating soil texture in China [71], which indicates the potential of SMAP SM data as an indirect method for characterizing soil physical properties at global scale.

3.3. Effects of Land Cover Types on Soil Moisture Estimation

Figure 4 displays the error metrics for surface and subsurface SM estimates calculated over various land-cover types against in-situ SM measurements from SCAN and USCRN networks for the test set. Overall, in terms of R and ubRMSE values, the model performed better for estimating subsurface SM than surface SM for most land-cover types, where the lower unRMSE values were obtained over savanna and woody savanna land-covers, followed by grasslands. Compared to the NLDAS-2 SM product, the ConvSLTM model performed better in terms of R values over shrublands, croplands, permanent wetland, savanna, and woody savannas. Regarding ubRMSE, the model outperformed NLDAS-2 SM product in permanent wetlands, savanna and woody savanna. In terms of ubRMSE values, the model estimates from scenario 2 outperformed scenario 1 estimates over permanent wetlands and savanna land-covers, while both scenarios performed similarly over croplands and shrublands. In grasslands and woody savannas, scenario 1 outperformed scenario 2. Over croplands, as a major land-cover type, and in terms of R and RMSE values, the ConvLSTM for both scenarios 1 and 2 performed better than that of NLDAS-2 SM product, though the model was trained based on NLDAS-2 SM products. In terms of mean bias, the ConvLSTM model shows the smallest values for surface and subsurface depths in most land-cover types, especially in croplands, shrublands and savanna. However, due to low learning ability of the model or low accuracy of NLDAS-2 SM products, the bias is large in other land covers like grasslands and permanent wetlands.

3.4. Spatiotemporal Validation of Soil Moisture Estimates

We further evaluated the accuracy of SM estimates from scenario 2 with in-situ SM measurements across the CONUS. Here, we are showing the result only for scenario 2 as it outperformed scenario 1 in most cases because of the effect of soil physical attributes and incorporating estimated surface SM as an additional predictor for subsurface SM estimation. Figure 5 illustrates the spatial patterns of the error metrics for each method, showing good agreement between the estimated and measured SM in terms of all error metrics for most of the stations, especially in western, northwest and southwest of the CONUS, where the model resulted in median R values of 0.5 and median ubRMSE values of 0.05 cm³ cm^-3 (surface depth) and 0.04 cm³ cm^-3 (subsurface depth). For example, in Utah state where the density of ground stations is high and because of topography the SMAP SM struggle to effectively capture variability [36], the model performed well. However, the bias values are large (more negative) in the southeast and mid-west regions, indicating that the model tends to underestimate surface and subsurface SM in these regions, while overestimating in the western CONUS areas. Based on the histograms of R and ubRMSE values, the model performed slightly better for estimating SM in the subsurface than the surface layer.

In Figure 6, all SM estimates from the ConvLSTM model based on scenario 2 for surface and subsurface depths are plotted and compared with the in-situ SM measurements. As shown, the performance for surface and subsurface depths is approximately similar, with correlation values of 0.72 and 0.75 and RMSE values of 0.09 and 0.10 cm³ cm^-3 for surface and subsurface depths, respectively. Visually, the model overestimates SM for dry conditions while underestimates SM for wet conditions.

To better understand the performance of the ConvLSTM model for temporal estimation of SM, we plotted the estimated surface and subsurface SM against in-situ measurements from three SCAN sites covered by cropland, grassland, and savanna land-covers across sandy loam (coarse-textured) and silt loam (medium-textured) soils. As shown in Figure 7, SM estimates from both scenarios well capture in-situ SM temporal dynamics, especially at low SM values, and clearly respond to precipitation events. For the selected three sites, the SM estimates and NLDAS-2 SM underestimate in-situ measurements, particularly in wet soil conditions. With respect to the cumulative density function (CDF) curves, the ConvLSTM model underestimates SM in humid sites while slightly overestimates SM in dry sites. In addition, scenario 2 slightly outperformed scenario 1 in both selected sites in arid and humid conditions.

4. Summary and Conclusions

Accurate information about surface and subsurface soil moisture (SM) with adequate temporal resolution (e.g., daily) at large scale is critical for a range of applications in hydrology, weather, and agriculture. Advances in deep learning techniques for data fusion and spatiotemporal modeling offer exceptional opportunity for integrating multi-source multi-scale data and capturing complex spatial and temporal relationships between SM and land surface parameters.

In this study, we developed a new method for estimating surface and subsurface SM using convolutional neural network and long short-term memory (ConvLSTM) deep learning techniques. ConvLSTM can capture both spatial and temporal nonlinearities, making the method ideal for modeling the intricate dynamics of SM. To do so, we implemented a method that integrates SMAP satellite SM and its ancillary variables with process-based NLDAS-2 SM products. This fusion leverages the strengths of both datasets, enabling the daily estimation of surface (5-cm) and subsurface (25-cm) SM. We employed ConvLSTM model in conjunction with two scenarios of inputs: (1) SMAP ancillary data, and (2) SMAP ancillary data with soil physical attributes from SOLUS100 digital maps. We further evaluated the effect of including surface SM estimates for estimating subsurface SM. The predictors included SMAP ancillary data like brightness temperature, albedo, soil temperature, surface roughness coefficient, vegetation opacity, and vegetation water content, in conjunction with SOLUS100-based soil physical properties. The model was trained using SMAP DCA and NLDAS-2 SM products, aiming to maximize correlation (R) and minimize RMSE, respectively. Extensive validation using in-situ SM data from SCAN and USCRN stations across the CONUS demonstrated that the ConvLSTM model achieved high accuracy for both surface and subsurface SM, closely aligning with observed data. The inclusion of soil properties (i.e., sand, silt, clay and bulk density) significantly enhanced model performance, particularly in coarse-textured soils and cropland areas, where the model exhibited low unbiased RMSE (ubRMSE) values. Additionally, the model showed superior performance in estimating subsurface SM compared to surface SM, especially when surface SM estimates are used as an additional predictor, across various land-covers, especially in areas covered by savanna and woody savanna.

The proposed approach overcomes the limitations of SMAP satellite 2-3 days overpasses, provides daily SM estimates for surface and subsurface, and reduces the extensive input requirements typical of process-based models. Consequently, this integration facilitates the production of comprehensive global SM estimates for both surface and subsurface depths. It highlights the effectiveness of integrating satellite-based and process-based SM retrieval methods within a ConvLSTM deep learning framework for robust SM estimation at large scale.