Enhancing Deep Learning Soil Moisture Forecasting Models by Integrating Physics-based Models

Soil moisture (SM) plays an important role in climate and hydrological systems by balancing the interaction of water and energy exchange processes (Seneviratne et al., 2010; Crow et al., 2018). Thus, accurate predictions of SM could improve various crucial applications, e.g., drought monitoring and water resource management (Dorigo et al., 2017).

Physics-based (PB) models have long been used to forecast SM (Kanamitsu et al., 2003; Wood and Lettenmaier, 2006), and show excellent potential in simulating SM dynamics (Xia et al., 2014; Esit et al., 2021). Therefore, PB models have been widely used for early warning of drought (Luo et al., 2007) and for identifying SM–precipitation feedback (Santanello et al., 2019). However, current PB models still have some flaws, such as uncertainties in model parameters and insufficient representation of land-surface processes (Brooks et al., 2015). On the other hand, deep learning (DL) models are well known for their ability to learn nonlinear mapping relationships and show a remarkable capability in short-term SM forecasting (Fang et al., 2017; Fang and Shen, 2020; Li et al., 2022a). For example, Fang et al. (2017) utilized the long short-term memory (LSTM) model to extrapolate Soil Moisture Active Passive (SMAP) SM data spatiotemporally, based on meteorological forcing, and concluded that LSTM outperformed traditional machine learning models (e.g., Random Forest). In addition, Li et al. (2021) utilized the convolutional LSTM (ConvLSTM) model, which can capture the information of both temporal and spatial dimensions, to forecast SM and showed that the ConvLSTM model outperformed the LSTM model in short-term SM forecasting.

To take advantage of both PB and DL models, various studies have developed so-called hybrid methods to incorporate physical information into DL models. Hybrid models have many categories [see Slater et al. (2023) for a comprehensive study], and one popular and effective approach is to utilize the output of PB models to enhance DL models. There are two main methods for achieving this. (1) Averaging the outputs of the PB and DL models is the most commonly used hybrid method because of its simplicity. Fang et al. (2019) averaged the SM simulated by LSTM and the Noah land surface model (LSM) and confirmed that this model combination could effectively improve the long-term simulation of SM. However, simple averaging can transfer the system error of PB models into DL models, which is not the optimal way to exploit the value of both models. (2) Another straightforward way is to feed the output of PB models as features into DL models. For example, Daw et al. (2022) enhanced the ability of the DL model to predict lake temperature by using the output of the general lake model as input features for the neural network. Cui et al. (2022) proposed an LSTM-based encoder–decoder model that improved the accuracy of multi-step-ahead flood forecasting by adding the output of the Xinanjiang hydrological model as input features into the decoder. Similarly, Klocek et al. (2022) used the output of the HRRR (high-resolution rapid refresh) model as an exogenous input into the decoder of the ConvLSTM-based encoder–decoder model, and significantly improved the long-term predictability of precipitation.

However, previous hybrid models have predominantly depended on DL models to automatically extract significant features from PB model outputs without fully capturing the essential spatiotemporal features in different regions and forecast time scales. The attention mechanism is typically used to extract important spatiotemporal patterns from various features for SM predictions (Li et al., 2022a). However, to the best of our knowledge, the attention mechanism has not been used to incorporate physical information into DL models for improving SM forecasts. By relying on the attention mechanism, it is possible to adaptively learn important spatiotemporal patterns from PB model outputs to aid DL models in making predictions at each forecast time scale.

Moreover, while various hybrid methods have been widely used and demonstrated promising performance in hydrological forecasting (Slater et al., 2023), including flood (Speight et al., 2021), streamflow (Liu et al., 2021), and precipitation (Klocek et al., 2022), no attempts have been undertaken to develop a model that leverages the benefits of these hybrid methods specifically for SM predictions. Ensemble methods have been proven to combine the strengths of different DL models and enhance predictability regarding each DL model (Lee et al., 2022). Therefore, ensembling the forecasts of different hybrid models can unify their predictions in a single model and reduce the errors from every single model.

In the present study, we aimed to effectively incorporate the outputs of PB models into DL models for accurate multi-step-ahead (1–16 days) SM predictions. To achieve this, we first developed an attention-based hybrid model to leverage the benefits of PB models at each forecast time scale. Then, we further developed an ensemble model that combined the advantages of different hybrid models to achieve exceptional multi-step-ahead SM predictions. More specifically, we (1) verified whether attention mechanisms can adaptively take advantage of PB models at different forecast time scales and regions to improve SM predictions spatiotemporally; (2) verified whether ensemble methods can utilize the advantages of different hybrid models to create an optimal model among all forecast models; and (3) thoroughly analyzed and discussed the pros and cons of different hybrid SM forecasting methods for multi-timescale predictions.

The remainder of the paper is organized as follows. Section 2 describes the data used in the study and the preprocessing procedures. Section 3 describes the proposed hybrid models and existing practical hybrid models. The evaluation processes, based on conventional metrics, triple collocation, and extreme indices are also shown in this section. Section 4 presents the evaluation results with respect to in situ and gridded data, along with some further discussion of our findings. Lastly, section 5 summarizes our conclusions.

The training target surface SM (0–5 cm) was obtained from the SMAP L4 dataset (Reichle et al., 2017), which is commonly used for SM simulation based on DL models (Karthikeyan and Mishra, 2021; Cai et al., 2022; Kannan et al., 2022; Li et al., 2022b). The SMAP satellite measures surface SM globally using a passive L-band radiometer (Entekhabi et al., 2010), and it has been the most promising satellite for SM monitoring due to its higher capacity for penetrating vegetation at the L-band compared to that at the C- or X-band (Wigneron et al., 2018). The SMAP L4 data are obtained by merging SMAP observations with the estimates from the Catchment LSM through the ensemble Kalman filter (Entekhabi et al., 2010), and provides surface SM at 3-hourly and 9 km resolution globally. The SMAP L4 data were averaged into a daily time scale.

Meteorological forcing (i.e., precipitation, longwave radiation, specific humidity, surface pressure, downward shortwave radiation, surface temperature, and wind speed), and static physiographic attributes (i.e., land cover, elevation) were used as input features fed into DL models. Climate forcing was derived from the SMAP L4 dataset [note that the climate forcing in this dataset is not from the SMAP satellite but is generated from NASA’s GEOS (Goddard Earth Observing System) model]; the land-cover data were obtained from the United States Geological Survey (USGS); and the elevation data were from MERIT DEM (Yamazaki et al., 2017). The physiographic attributes were spatially mapped at 9 km resolution for consistency with the SMAP data. We used input features with a previous time of 1–7 days and predicted the SM for the next 1–16 days.

The PB model for the SM forecast was obtained from Global Forecasting System (GFS). The GFS forecast products are publicly distributed by NOAA (the National Oceanic and Atmospheric Administration). The core LSM in GFS is the Noah model, and the GFS surface SM (0–10cm) forecast is widely adopted by different national organizations (e.g., the USGS and the National Meteorological Centers in Asia). GFS provides up to 16-day forecasts at a 0.25° spatial resolution with forecast time steps of 3 h from 0 to 240 h, and 12 h from 240 to 384 h. GFS SM forecasts were temporally averaged into a daily time scale and bilinearly interpolated into a 9 km resolution.

In situ SM observations used for evaluation were obtained from the China Meteorological Administration (CMA). The observations contain 2966 stations over China at an hourly temporal resolution from 1 January 2010 to 31 December 2018. Unrealistic data (e.g., discontinuities and unchanged data) were first labeled (Dorigo et al., 2013), and then observational stations with more than 5% missing data were excluded from our analysis. These requirements led to a reduction in the number of sites from 2966 to 1164 for further analysis [Fig. S1 in the electronic supplementary material (ESM)]. The stations were densely distributed in southeastern China but sparsely distributed in western China. In total, the selected stations covered all five major climate zones (i.e., tropical, arid, temperate, cold, and polar).

Figure 1.  The detail of the (a) ConvLSTM-ED model and (b) inner structure of ConvLSTM.

Two gridded datasets used for triple collocation analysis (TCA) were obtained from ERA5-Land and SoMo.ml. We selected these two datasets because of their high quality and because the different types and sources of these two datasets led to error independence that could ensure the assumption of TCA (Stoffelen, 1998). The hourly, 9 km surface SM (0–7cm) obtained from ERA5-Land (Muñoz-Sabater et al., 2021) was averaged into a daily time scale for further analysis. Evaluation against in situ observations shows the value of ERA5-Land in the description of the hydrological cycle, particularly with SM (Beck et al., 2021; Muñoz-Sabater et al., 2021). SoMo.ml data, on the other hand, were derived through a DL model based on the training dataset from in situ observations (O and Orth, 2021) and provided daily surface SM (0–10 cm) at a 0.25° spatial resolution from 2000 to 2019. SoMo.ml data show excellent performance, particularly in terms of temporal dynamics, compared with existing SM data (e.g., ERA5).

The ConvLSTM-based encoder–decoder model (ConvLSTM-ED) was used to provide the next 1–16 days of SM predictions, and the hybrid models were then built to incorporate the SM forecasts of the GFS model into the ConvLSTM-ED model.

The encoder–decoder model consists of two neural networks designed for multi-step-ahead predictions. During the encoding process, critical information from the input features is compressed into an encoding vector. The decoding process then converts the encoding vector into the target forecasts. The encoder–decoder model can consider the temporal association among sequences (Cui et al., 2022), making it suitable and extensively used for multi-step-ahead SM forecasting (El Saadani et al., 2021; Li et al., 2021).

We used ConvLSTM in a recursive encoder–decoder structure (Cho et al., 2014; Cui et al., 2022) to predict SM in the next 1–16 days (Fig. 1a). In the decoder, the input of each time step is the output from the previous time step. Through the recursive process, the decoder can receive not only the important information extracted from the encoder but also the changing information from the previous time step. Both the encoder and decoder adopted ConvLSTM, and the principle of the ConvLSTM cell (Fig. 1b) is shown as follows:

where $ t $ represents the $ t $th step of the input timesteps; $ {\mathbf{X}}_{t} $ denotes the input features; $ {\mathbf{C}}_{t} $ is the cell state and $ {\mathbf{H}}_{t} $ is the hidden state; $ * $ denotes the convolution operator; and $ \cdot $ denotes the dot product. Notably, compared to the LSTM model, the ConvLSTM model adds the convolution operator into the state-to-state and input-to-state transitions to further encode spatial information. $ \sigma $ is the $ \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{d} $ function, and $ \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h} $ represents the hyperbolic tangent function. $ {\mathbf{W}}_{\mathrm{x}\sim} $, $ {\mathbf{W}}_{\mathrm{c}\sim} $ and $ {\mathbf{W}}_{\mathrm{h}\sim} $ are 2-dimensional convolution kernels, and $ {\mathbf{b}}_{\sim} $ is bias.

Two widely used hybrid models were adopted for benchmarking our proposed hybrid models. Firstly, SM was forecasted by simply averaging the SM forecasts for the next 1–16 days of both GFS and ConvLSTM-ED [namely, the average model, from Fang et al. (2019)]. Secondly, the GFS SM forecasts were fed into the decoder of ConvLSTM-ED to make the next 1–16 days of forecasts [namely, the condition model, from Cui et al. (2022) and Klocek et al. (2022)] (Fig. 2a).

Figure 2.  The model structures of (a) the condition model, (b) the attention model, and (c) the attention block. The dimension of each intermediate feature of the attention block is annotated. H, W and F are the sizes of the height, width and features dimensions, respectively. The abbreviations of intermediate features are the same as given in the main text. The colors in (c) indicate weights.

Based on the condition model, we developed the attention model (Fig. 2b). We used an attention block (Fig. 2c) to adaptively extract the important patterns from GFS SM forecasts at different forecast time scales. The extracted patterns were then used as the input for the decoder to create SM forecasts. The attention block was calculated through the following four steps. (1) For each forecast timestep t, the input tensors $ {\mathbf{U}}_{t} $ were first averaged over their spatial dimensions and reduced to a vector of features $ {\mathrm{z}}_{t} $ ($ {\mathbf{F}}_{\mathrm{s}\mathrm{q}} $). (2) A fully connected layer was applied to obtain the vector of learned weights $ {\mathrm{s}}_{\mathrm{t}} $ of the compressed tensor $ {\mathrm{z}}_{\mathrm{t}} $, representing the reassigned importance of the original input features ($ {\mathbf{F}}_{\mathrm{e}\mathrm{x}} $). (3) The learned weights were multiplied by the input features to obtain a final importance-reassigned tensor $ {\tilde{\mathbf{U}}}_{t} $ ($ {\mathbf{F}}_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}} $). (4) The final step is the output process, in which we send the input tensors $ {\mathbf{U}}_{t} $ and importance-reassigned tensor $ {\tilde{\mathbf{U}}}_{t} $ into a fully connected layer to obtain the final output ($ {{\mathbf{U}}^{\mathbf{*}}}_{t} $) of the attention block. For the ensemble model, we ensemble-averaged the forecasts of the average, condition and attention models to utilize the advantages of different hybrid methods. Notably, the ensemble model was constructed by simple averaging rather than using weighted ensemble methods (e.g., Bayesian model average). The aim of our study is to emphasize the benefits of integrating different hybrid methods (even using the simplest ensemble methods), rather than exploring the advantages of using more advanced ensemble methods.

The domain of this study was restricted to China (14.7°–53.5°N, 72.3°–135°E). We divided this region into 24 patches, and each patch contained 112 × 112 pixels. A local DL model was trained for each patch, and the training period was from 31 May 2015 to 31 December 2017 (80% and 20% of the dataset were used for training and validation, respectively), and the testing period was from 1 January 2018 to 31 December 2018. Before the training process of models, we normalized all input features to speed up the convergence. The backward propagation was based on the Adam gradient descent method (Kingma and Ba, 2017). The loss function was set as the mean squared error (MSE) between the observations and predictions.

Two types of hyperparameters were tuned, including the parameters that determine the structure of DL models (i.e., the number of hidden cells), and the parameters that optimize DL models during the training process (i.e., batch size, epochs, learning rate) (the detail is shown in Text S1 in the ESM). We tuned the hidden sizes of [8, 16, 32, 64] and the batch sizes of [16, 32, 64], and selected the “best” model structure with the minimal MSE on the validation dataset. After tuning, the number of hidden cells of the ConvLSTM models was set to 16, and the batch size was set to 32. The learning rate was set to 0.001 and was adaptively updated by multiplying the reducing factor (set to 0.1) when validation loss ceased improving. DL models were optimized for 50 epochs by selecting the “best” model with the lowest validation loss. Our study indicated that 50 epochs are sufficient for finding the optimal model as the validation loss stopped decreasing.

The DL models were trained on a single NVIDIA A100 GPU, and the training processes of each hybrid model cost about 120 h (for nearly 3 years of training data). While training requires a one-time effort, it significantly reduces the runtime during inference. For real-time applications, the model can generate 16-day SM forecasts over China in just a few minutes using one A100 GPU.

Two performance assessment criteria were used for in situ evaluation: the Pearson correlation coefficient (R) and unbiased root MSE (ubRMSE) (Text S2 in the ESM). R measures the correspondence between forecast and in situ time series variables in terms of temporal variability. ubRMSE is a scaled error criterion for describing the differences between observed and predicted time series. We evaluated the performance using the SM forecast time series of the nearest grid of each in situ station.

TCA (Gruber et al., 2016) was further used to evaluate the error of the forecast SM toward an unknown truth using collocated datasets (the detail is shown in Text S3 in the ESM). TCA can directly estimate the sensitivity of data using covariance and variance. Here, we used the TCA-based signal-to-noise ratio (SNR) to represent the error of SM forecast models. The TCA-based SNR is calculated using the following equation (Kim et al., 2020):

where $ {\mathrm{\sigma }}_{\mathrm{\Phi }}^{2} $ is the variance of the unknown truth of SM, and $ {\mathrm{\sigma }}_{{{\varepsilon }}_{\mathrm{i}}}^{2} $ is the variance of the random error of dataset i. The TCA-based SNR is symmetric around zero, and the zero value indicates that the SM signal ($ {{\beta }_{\mathrm{i}}}^{2}{\sigma }_{\mathrm{\Phi }}^{2} $) of forecast data is identical to the noise ($ {\sigma }_{{{\varepsilon }}_{\mathrm{i}}}^{2} $).

The soil wetness deficit index (SWDI) is used to estimate the severity of drought (Text S4). The SWDI is commonly used in agricultural drought assessment (Martínez-Fernández et al., 2015; Yang et al., 2017; Li et al., 2022b) and is calculated as

where $ \theta $ is the forecasted SM, and $ {\theta }_{\mathrm{F}\mathrm{C}} $ and $ {\theta }_{\mathrm{W}\mathrm{P}} $ are the field capacity and wilting point, respectively. In our study, the 5th and 95th quantile values of the historical SM time series were used to represent the wilting point and field capacity. We defined SWDI > 0 as no drought, −5 < SWDI < 0 as moderate drought, and SWDI < −5 as extreme drought, based on Zhang et al. (2021). The probability of detection (POD) metric (Ford and Quiring, 2019) was applied to answer how well the model correctly identifies drought events. POD is calculated as follows:

where “hits” represents the number of observed drought events that the SM forecast models correctly depicted (i.e., we observed extreme drought based on in situ observations and also predicted it in the forecast models) and “misses” represents the number of observed drought events that SM forecast models do not correctly depict (i.e., we observed extreme drought based on in situ observations but did not predict it, or predicted only moderate drought, in the forecast models). POD ranges from 0 to 1, with 1 representing a perfect score.

We first evaluated GFS, ConvLSTM-ED, and different hybrid models against in situ observations (Figs. 3, 4). The ConvLSTM-ED model performed reasonably well in short-term predictions, with comparable predictability to SMAP L4 data (i.e., training target SM). However, its performance degraded dramatically as the forecast time scale increased (Fig. 3), particularly in southeastern China (Figs. 4g–i). The GFS model showed inferior performance compared to the ConvLSTM-ED model during all forecast time scales, especially for short-term forecasting (Fig. 3). Although the ConvLSTM-ED model outperformed the GFS model for different forecast time scales, it showed degraded performance in wet regions (e.g., southeastern China) compared to the GFS model (Figs. 4d–f, Figs. 4g–i). Notably, we cannot conclusively determine which model (GFS or ConvLSTM-ED) was superior based on the aforementioned results. Moreover, this determination was beyond the scope of our paper (the reasons are shown in Text S5).

Figure 3.  The mean (a) R and (b) ubRMSE of different forecast models at different forecast time scales. Dash lines denote the performance of SMAP L4 data evaluated by in-situ observations. The abbreviations of model names are the same as in section 3.

Figure 4.  The spatial distribution of performance (R) in 1-, 7- and 16-day forecasts of different models. We used the average model as the baseline hybrid model to evaluate the performances of the different hybrid models. Panels (a–c) show the performance of the average model, while the remaining rows show the differences between the R of the target model and the R of the average model. Red points indicate that the model improved the performance compared to the average model, while blue points show a declined performance.

The average model dramatically improved the ConvLSTM-ED model at nearly all forecast time scales, demonstrating the benefit of adding physical information into DL models. However, the performance of the average model still dropped dramatically for long-term forecasting (Fig. 3), particularly in southeastern China (Figs. 4a–c), consistent with the property of the ConvLSTM-ED model (i.e., decreasing performance in southeastern China as the forecast time scale increased). Moreover, the average model performed poorly in northern China, caused by introducing the bias of the GFS model in this region. These results indicate that the simple averaging method could not fully exploit the benefits of both the PB and DL models.

The condition model greatly improved the long-term predictability compared to the average model over nearly all stations (Fig. 3, Figs. 4j–l), demonstrating that adding GFS SM forecasts as the exogenous inputs into the decoder of the ConvLSTM-ED model could significantly improve long-term predictions. However, the performance of the condition model decreased significantly in short-term predictions compared to the average model (Fig. 3). In short-term predictions, we observed that the underperforming regions in the condition model (blue dots in Fig. S3b) were consistent with those in the GFS model (red dots in Fig. S3a in the ESM) when compared to the ConvLSTM-ED model. This result highlighted a problem with the condition model, i.e., it introduced the biases from the short-term predictions provided by the GFS model into the DL models. Notably, although the condition model could propagate short-term forecast errors, incorporating the PB SM evolution (i.e., sharpened predictions from GFS) in the decoder could still improve the performance of the ConvLSTM-ED model significantly. This emphasizes the significance of integrating physically consistent predictions into pure DL models.

The attention model further improved the short-term predictability and had equal long-term predictability when compared to the condition model (Fig. 3). Spatially, the attention model retained the high predictability of the ConvLSTM-ED model in northern China and effectively overcame its deficiencies in southeastern China (Figs. 5g–i). Moreover, the attention model outperformed the condition model in most regions (Figs. 5j–l). In particular, it showed significant improvement in short-term predictions where the condition model introduced the bias of the GFS model (Fig. 5). These findings indicate that the attention mechanism can adaptively learn to exploit the benefits of the ConvLSTM-ED and GFS models for different forecast time scales and soil conditions, thereby significantly improving the model performance spatiotemporally.

Figure 5.  The R of the (a–c) GFS, (d–f) ConvLSTM-ED, and (g–i) attention models, and (j–l) the improvement of the attention model compared to the condition model.

The ensemble model outperformed all other hybrid models (in terms of both R and ubRMSE) at all different time scales (Fig. 3), especially in long-term predictions; for example, the ensemble model improved the mean R by 65% (from 0.205 to 0.340) compared to the ConvLSTM-ED model for the 16-day predictions. Moreover, the ensemble model was able to further reduce the bias of the predictions in southeastern China (provided by the ConvLSTM-ED model) compared to the attention model (Figs. 4p–r), and the ensemble model outperformed the ConvLSTM-ED model over 79.5% of the in situ stations. These results underline the value of ensemble methods and emphasize the exceptional spatiotemporal predictability of the ensemble model.

We further evaluated the predictability using gridded datasets by TCA, which could evaluate the performance toward an unknown truth. The spatial distribution of SNR is shown in Fig. 6. The results were similar to those of the in situ evaluation, and also demonstrated the superior performance of the proposed hybrid models (attention and ensemble models). For example, the condition model enhanced the long-term predictability but decreased the short-term predictability (Figs. 6j–l to Figs. 6p–r). In addition, the attention model was able to further correct the bias of the condition model in short-term predictions (Figs. 6j–l to Figs. 6m–o), whereas the ensemble model achieved the best performance among all hybrid models (Fig. 6). We further used another collocated dataset (i.e., SMOS data, in situ data from the CMA) and evaluated the predictability by TCA. The result was also consistent with that of the in situ evaluation (Text S6, Fig. S4 in the ESM).

Figure 6.  TCA-based SNR of different models. The triplets of the TCA are [*, ERA5-Land, SoMo.ml], where * denotes the forecast models. Panels (a–c) show the average model. The remaining rows show the difference between the SNR of the target model and the average model.

Figure 7 shows the SNR of SMAP L4 data (i.e., the training target of the DL models) and the ConvLSTM-ED and ensemble models. The ConvLSTM-ED model showed inferior performance compared with SMAP L4 data in most regions. Additionally, the regions where the ConvLSTM-ED model underperformed were consistent with those of SMAP L4 data. It was found that the performance of the ConvLSTM-ED model depended entirely upon the quality of the SMAP L4 data and was limited by the SMAP L4 data as a performance ceiling. Notably, the ensemble model outperformed the SMAP L4 data in most regions for short-term predictions (Fig. 7), particularly in drought-prone areas (e.g., the North China Plain), suggesting that the ensemble of different hybrid models could “break” the performance ceiling constrained by the training data in some areas. This is attributable to the introduction of physical information into the pure DL models. The in situ validation of the SMAP L4 data and the ensemble model further confirmed this result (Fig. S5 in the ESM). However, the long-term predictability of the ensemble model was still far inferior to the SMAP L4 data. Moreover, all forecast models still did not show satisfying performance (i.e., signal larger than noise, SNR > 0) in more than half of the regions (Fig. S6c in the ESM) in long-term predictions, indicating the challenge of long-term SM forecasting, which necessitates further investigation.

Figure 7.  TCA-based SNR of (a-c) ConvLSTM-ED, (d–f) ensemble model predictions and (h) the SMAP L4 datasets. The triplet of TCA is [*, ERA5-Land, SoMo.ml], where * denotes the forecast models and SMAP L4.

We further evaluated the drought predictability of the different forecast models. Figure 8 illustrates the kernel density curves of the SWDI of the in situ observations and different models. Surprisingly, the SWDI of the in situ observations contained two peaks located at the SWDI values of nearly −10 and −2. We further used more stringent quality-control processes for the in situ observations (Dorigo et al., 2013) and found the same two peak structures (Fig. S7 in the ESM). It is noted that these two peaks may be a unique property of the in situ observation datasets used in our study. We found that the GFS and average models tended to forecast the right peak of the in situ observations (i.e., they gave relatively stable predictions and were unaware of some extreme events such as SWDI < −10). On the contrary, the attention model tended to simulate the left peak (i.e., extreme drought events) better than the other hybrid models, showing the effectiveness of the attention mechanism for extreme drought forecasting. Furthermore, although the ensemble model provided the best general performance, it tended to forecast the mean SWDI of observations. This result emphasizes that ensemble methods could provide a more stable prediction by correcting the bias of each model but may also “remove” some extreme events, which may not be suitable for drought forecasting.

Figure 8.  The kernel density curve of the SWDI of the in situ observations from different forecast models (lines with different colors) at the (a) week-1 and (b) week-2 forecast.

We further evaluated the fractions of observed drought events that were correctly depicted by the forecast models. Table 1 summarizes the POD values of the different models over different Köppen–Geiger major climate zones. Generally, the attention model was able to accurately detect 60.6% and 56.8% of drought events at 1- and 2-week forecasts and achieved the best detection over arid, temperate, cold, and polar regions. Moreover, the ensemble-average operation always yielded an average prediction of drought events among ensemble members (see average and ensemble models), reinforcing the prior results. It is noted that the GFS model excelled in temperate regions but performed the worst over arid, cold, and polar regions among all the forecast models, indicating a poor representation of SM dynamics over these regions.

It was found in this study that embedding physical information in the DL models through useful hybrid methods dramatically improved the SM predictability compared to using pure DL models, and this could be attributed to several possible reasons. Firstly, it is well known that pure DL models may produce unrealistic predictions because of a lack of physical consistency (e.g., mass and energy balance). For example, Fang et al. (2019) found that pure DL models provided highly fluctuating simulations non-physically. Thus, physical information provided by PB models that obeys physical laws can be used to correct the non-physical predictions of pure DL models. Secondly, pure DL models can benefit from the assimilation of high-quality observations in PB models (Fang and Shen, 2020); for example, pure DL models cannot predict the corresponding SM variation if a rainfall event is missing in the forcing data. However, data assimilation can remedy the forcing errors with high-quality observations, resulting in better temporal representation of SM dynamics. One benefit of using the GFS (including data assimilation) forecasts in our study was to help the pure DL models in correcting the bias induced by the forcing errors. Thirdly, Daw et al. (2022) pointed out that pure DL models rely heavily on the training data quality and can only depict the evolution of existing SM (Klocek et al., 2022). This may lead to significant biases over regions with poor-quality data (e.g., wet regions in SMAP L4 data). On the contrary, PB models can depict the dynamics of SM over different soil conditions (e.g., precipitation infiltrates more easily in regions with high soil porosity), and can provide stable and realistic simulations with high-quality rainfall forcing [e.g., wet regions; see Maggioni et al. (2012)]. In addition, the GFS model can simulate SM in different water states (e.g., solid, liquid) through soil dynamics, which pure DL models struggle to do accurately because of the poor quality of the training datasets during the freezing period. Thus, incorporating physical information into pure DL models might help to overcome the deficiencies derived from data (Daw et al., 2022).

Although we introduced PB features to improve the model, the proposed hybrid models still inherited the uncertainties from the supervised DL models, i.e., the uncertainty from the training data. In addition, another source of uncertainty came from the selection of hybrid schemes, as demonstrated in section 4. Furthermore, the quality of the PB models also contributed to the uncertainty. Parameterizations and inadequate representation of land processes can introduce uncertainties in hybrid models. However, when compared with the PB models, the hybrid models benefited from the fitting ability of the DL algorithm and the vast amount of data, which could partially correct systematic errors. Moreover, the introduction of PB features also alleviated the limitation of the training data when compared to the pure DL models. These findings suggest that hybrid models are a promising way of enhancing the prediction skill for meteorological and hydroclimatic variables (Slater et al., 2023).

The potential applications of SM forecasting models have been comprehensively discussed in Peng et al. (2021), and we highlighted two important application directions. Firstly, the proposed model could provide accurate initializations of land-surface conditions for numerical weather prediction (NWP) systems. Indeed, the integration of SM into several NWP models has been found to improve forecasts of atmospheric variables (Dharssi et al., 2011; Muñoz-Sabater et al., 2019; De Rosnay et al., 2020). Secondly, accurate predictions of SM could be utilized for monitoring, analyzing and providing early warnings of hydrometeorological disasters, including agricultural drought (Mishra et al., 2017) and floods (Li et al., 2018). Additionally, these predictions could inform decision-making processes, such as in watershed management (Heimhuber et al., 2017) and irrigation water management (Lawston et al., 2017).

In our study, we aimed to investigate the benefits of incorporating physical information into DL models, but exploring the interpretability of the proposed models is beyond the scope of the present paper. However, these complex hybrid models may have low interpretability and should be used with caution in practical applications. Explainable artificial intelligence (XAI) provides tools to aid in decision-making processes when applying DL models in real-world applications. Several studies have explored the interpretability of DL SM forecasting models using XAI tools. For example, Huang et al. (2023) adopted various post-hoc interpretation methods to assess the feature effects on SM predictions and showed that a comprehensive understanding of the relationship between input features and predicted SM could be achieved. Different interpretation methods used in their study, such as “shapely values” and “partial dependence plots”, could be used to investigate the contributions of different features (e.g., GFS forecasted values) to our proposed models, which deserves further exploration.

We end our discussion by pointing out some limitations of our study. Firstly, we did not provide the “best” hybrid schemes to achieve the “best” forecast (i.e., general performance and drought predictability) at different forecast time scales and spatial regions. For example, the ensemble model achieved the best general performance at all forecast time scales (section 4.1), but the ensemble method may “remove” some extreme drought events (section 4.3). Therefore, we highlight that the choice of different hybrid methods might depend on the different applications, such as the ensemble model is suited for long-term, stable predictions, which mainly focus on the average state of SM, while the attention model is suited for forecasting extreme drought events. Secondly, we integrated GFS with the ConvLSTM-ED models because of its efficiency and widespread use (Fan and van den Dool, 2011; Yin et al., 2019). However, the GFS and ConvLSTM-ED models were both not the “best” PB and DL models for SM prediction. Thus, the result of our study may not fully represent the properties of PB and DL models. Nonetheless, we showed improvements in the different hybrid methods based on these two widely used models. Thirdly, hybrid models have different framework, e.g., physically guided DL (Willard et al., 2022a), or differentiable programming (Feng et al., 2022). In this study, we only focused on using the PB model outputs and observational features in a hybrid modeling setup to generate strong-performing SM predictions. We did not introduce any physical laws and principles to guide the DL models. Several “deep” hybrid frameworks have been developed (Read et al., 2019; Liu et al., 2022), which can “force” DL models to forecast based on physical consistency, thereby possibly providing more realistic and stable predictions (Willard et al., 2022a). Moreover, pre-training DL models using PB model outputs and fine-tuning them in the target data (i.e., transfer learning) may also utilize the physical information.

In this paper, we first propose an attention hybrid model based on condition hybrid schemes and an attention mechanism to utilize the advantages of both PB and DL models over different forecast timescales and regions. An ensemble model is then further proposed by averaging the outputs of two existing practical hybrid models (average and condition) and the proposed attention model. To the best of our knowledge, this is the first study taking both the attention mechanism and ensemble methods to integrate PB and DL models for SM forecasting. We thoroughly assessed the predictability of the two proposed hybrid models (i.e., attention and ensemble) and two existing hybrid models (i.e., average, and condition) based on in situ and gridded data evaluation. Generally, the proposed hybrid models outperformed the two existing hybrid models, and could greatly improve the long-term predictability, and predictability of drought events, compared to pure DL models. The main conclusions were as follows:

(1) The proposed ensemble hybrid model achieved the best general performance among all hybrid models under different soil conditions over all forecast timescales (from 1 to 16 days), especially for long-term forecasting. Notably, the ensemble hybrid model improved 65% of the mean values of R and 6% of the mean values of ubRMSE for the 16-day forecast compared to the ConvLSTM-ED model, and the ensemble hybrid model outperformed the ConvLSTM-ED model over 79.5% of the validation stations.

(2) The proposed attention hybrid model achieved the best drought predictability among all hybrid models. This model could accurately detect 60.6% and 56.8% of drought events for 1- and 2-week forecasts, respectively, and had generally the best drought detection ability over arid, temperate, cold, and polar regions. The attention hybrid model was able to detect an additional 2.4% and 3% of drought events compared with the ConvLSTM-ED model for 1- and 2-week forecasts, respectively.

(3) Different hybrid schemes had their pros and cons, and our proposed model solved the problems encountered by existing hybrid methods to some extent. For example, the average model is simple and effective, but the performance significantly degraded as the forecast time scale increased. The condition model could significantly improve the long-term predictability of SM but involves the bias of the GFS model in short-term predictions. The proposed attention model solves the problem of the condition model and is suitable for forecasting extreme drought events. The proposed ensemble method performed best among all hybrid models, and is suitable for long-term, stable predictions, which mainly focus on the average state of SM.

Finally, we provide some future ideas for improving SM forecasting by relying on hybrid models that integrate PB and DL models together. Firstly, DL and PB models should be further improved for SM forecasting separately. For example, DL models should be trained with multi-source data based on multimodal learning to avoid reliance on single datasets. Secondly, it is best to use some physical laws to guide the design of the structure, the parameter initializations, and the loss functions of DL models, which could help provide physically consistent predictions. Thirdly, the attention mechanism shows an excellent ability to utilize the benefits of different input features. The attention mechanism used in this study only focuses on the channel features. Therefore, axial attention modules should be developed to adaptively extract essential features of PB and DL models over different axials, which may further enhance the ability to identify valuable information.

Acknowledgements. Lu LI was supported by the Natural Science Foundation of China (Grant Nos. 42088101 and 42205149); Zhongwang WEI was supported by the Natural Science Foundation of China (Grant No. 42075158); Wei SHANGGUAN was supported by the Natural Science Foundation of China (Grant No. 41975122); and Yonggen ZHANG was supported by the National Natural Science Foundation of Tianjin (Grant No. 20JCQNJC01660). All data, source codes and example codes are available at https://github.com/leelew/HybridHydro.

Electronic supplementary material: Supplementary material is available in the online version of this article at https://doi.org/10.1007/s00376-023-3181-8.