Skillful subseasonal soil moisture drought forecasts with deep learning-dynamic models

Domain: earth_science
Published: 2025-08-12
Editorial outcome: major_revision

Abstract
Deep neural networks that learn from climate reanalysis data have produced skillful weather forecasts within ten days. However, it is still a great challenge for dynamic models to predict soil moisture, droughts, and other extreme events with lead times beyond two weeks. Here, we combine a recursive deep learning model (namely RISE-UNet) and subseasonal forecasts from dynamic models and achieve skillful forecasts of root zone soil moisture up to four weeks in advance. Our hybrid model, combining RISE-UNet and dynamic model forecasts, outperforms reanalysis-driven RISE-UNet models, while both methods show significantly higher performance than the latest European Centre for Medium-Range Weather Forecasts (ECMWF) and Global Ensemble Forecast System (GEFS) dynamic models, the post processed ECMWF or GEFS subseasonal forecasts by RISE-UNet, or ensemble model output statistics. The hybrid model shows skill in predicting flash droughts, which is higher than ECMWF and GEFS models in most cases, as demonstrated for major events in the United States, China, and Australia. The forecast skill of the hybrid modeling approach from weeks three to four is mainly due to the inclusion of the first two-week dynamic model forecasts and antecedent root zone soil moisture reanalysis. Our results indicate that combining deep learning with dynamic model forecasts can substantially improve the skill of subseasonal predictions beyond two weeks, particularly for root zone soil moisture and flash drought events.

Introduction
Root-zone soil moisture (RZSM) is an essential climate variable which influences water and energy cycles and regulates thermodynamic exchanges between land and atmosphere, thereby influencing regional weather and climate conditions 1–3. Excessive RZSM drying, also known as drought, have great socio-environmental impacts related to food security, forestry, human health, environmental sustainability, or biological community structure 4–7. RZSM influences both land surface and atmospheric processes through mass and energy exchanges, and its variability can be influenced by regional seasonality, internal atmospheric variability from climate oscillations, or human-induced through anthropogenic greenhouse gas emissions or land cover change 8–10. Given its vital importance for controlling land-atmospheric influences water and energy cycles and regulates thermodynamic feedbacks, regulating water and energy processes, and determining extreme events such as drought and heatwaves, skillful RZSM forecasts at long leads are crucial for drought management, climate adaptation, mitigation, and decision making 11–13.

Previous studies have consistently identified limited subseasonal forecast skill beyond two weeks (>14 days) in current dynamic forecast systems and land surface models, including but not limited to: Noah-MP, Weather Research and Forecasting (WRF), Variable Infiltration Capacity (VIC), Global Ensemble Forecast System (GEFS), and subseasonal-to-seasonal European Center for Medium Range Weather (S2S ECMWF) 13–20. The subseasonal time scales ranges between ~14–60 days (~2–9 weeks), and low skill within models has been attributed to several factors, including unskillful precipitation forecasts, imperfect model parameterizations, initialization biases, and cumulative error growth associated with stochastic processes 13,21. To increase dynamic models (DM) forecast reliability and reduce model biases, post-processing bias-correction techniques are applied to DM forecasts with respect to observations (or reanalysis), which have improved subseasonal forecast skill of temperature, precipitation, wind speed, and organized motion associated with the Madden Julian Oscillation (MJO) 22–25. While bias correction can improve forecasts of mean conditions, substantial challenges remain in predicting extreme events, as bias-corrected forecasts often show limited skill beyond two weeks 26.

RZSM has a memory in which dry or wet anomalies, relative to the climatological mean, decays slowly for time periods between days to months and enhance land-atmosphere feedbacks through radiative and mass exchanges 27–29. Soil moisture memory (SMM) varies across geographic regions, seasons, and soil depths, and is influenced by dynamic factors such as solar radiation, vegetation, and precipitation, as well as static factors like slope and soil texture 30. While SMM can be a potential source of subseasonal predictability, previous studies have robustly identified that the prediction skill of dry extremes of RZSM is still limited beyond two weeks using either DMs, hybrid statistical-DMs, or machine learning methods 13,19,26,30. This low predictability can be attributed to large errors in precipitation forecasts, which affect root zone soil moisture (RZSM) through well-known processes such as infiltration and recharge, as well as to poor representation of precipitation in computing dynamic mass and energy fluxes. To improve subseasonal RZSM predictability beyond current land surface models that rely on error-prone precipitation forecasts, novel methods are needed to overcome this limitation.

Deep learning (DL) has the strength of learning complex systems by automatic extraction of spatial-temporal features 31. By directly learning from historical reanalysis data, DL has achieved notable successes in ten-day medium-range weather predictions 32,33, five-day extreme flood predictions 34, and seasonal predictions of El Nino South Oscillation 35 (ENSO), all showing improved performance compared to state-of-the-art DMs. Recent studies have also developed DL models for subseasonal predictions of atmospheric variables, while the method did not outperform state-of-the-art DMs, particularly at shorter leads (<10 days) 36. Due to SMM and improved medium-range weather prediction, DL has a great potential to skillfully predict subseasonal soil moisture and droughts beyond two weeks.

Few studies to date have explored the potential of deep learning (DL) for subseasonal prediction of root zone soil moisture (RZSM) and drought. In addition, DL-based post-processing for bias correction of spatial RZSM forecasts remains largely unexplored. Here, we present a novel hybrid approach that skillfully predicts subseasonal RZSM and droughts by combining a recursive DL model (namely RISE-UNet) with state-of-the-art dynamic models (DMs). We investigate whether DL models trained on publicly available data can provide accurate long-lead predictions at subseasonal timescales and whether hybrid DL-DM models outperform reanalysis-driven DL models, DL- or statistically post-processed DMs, and the original DMs themselves. To test this, we designed a set of experiments including: hybrid DL-DM models that use inputs from DM reforecasts and reanalysis; DL- or statistically post-processed DMs that apply bias correction to DM outputs; and DL models trained on reanalysis data alone. All models are trained to predict RZSM reanalysis as the target variable. Given the global availability of training data, our approach can be implemented in other regions to support skillful subseasonal forecasting of RZSM and drought. Moreover, due to the central role of soil moisture in subseasonal-to-seasonal predictability, this framework could be extended to improve forecasts of other climate extremes, such as heavy precipitation and heatwaves, potentially benefiting water resources drought response, agriculture, urban planning, and ecosystem management.

Results
Weekly soil moisture forecast skill and feature importance
In this study, our DL model, namely the residual inception squeeze-and-excitation (RISE) UNet++ (RISE-UNet, see Methodology), takes designed inputs and learns complex mapping functions between input variables and the RZSM forecast target (Fig. 1 and Supplementary Table S1). RISE-Unet accepts an arbitrary number of input predictors (channels) from both reanalysis at weekly lags and reforecast at forecast leads, and passes weighted convolution kernels across each channel to produce intermediate predictions, which are concatenated. Additional convolution kernels are passed along all intermediate predictions, and the image size is reduced through pooling operations, which allow different representations of the data to be computed within a more restricted field space, allowing the model to identify relevant features at different scales. These intermediate results are concatenated again, and the convolution/pooling process repeats until a minimal image size is reached (in our instance, a 48×96 pixel image is reduced to a 6×12 pixel image). Once the terminal size is reached, intermediate predictions are upscaled until the original image size is reached. For the RISE-Unet, three separate predictions are computed and the loss function applied over each one independently, with the goal of the model to reduce the aggregated loss of all three predictions.

We trained the model separately for each week, lead 1–5 and verified skill results against the physically-based GLEAM reanalysis (see Methodology). For clarification, skill results are obtained for each forecast week lead made at weekly intervals across the calendar year/season for the two-year period 2018-2019. Results consistently indicate that both DL-DM and DL models show superior performance compared to DM, DM-BC_DL, and ensemble model output statistics (EMOS) (Fig. 2A, B). For week lead 3, both DL-DM and DL models achieve the highest average anomaly correlation coefficient [ACC] of ~0.60 (an increase of ~66% over DM forecasts EMOS) and the lowest average CRPS of ~0.015 (>27% improvement) across the contiguous United States (CONUS). Hybrid DL-DM and reanalysis-driven DL models consistently outperform EMOS post processed DM forecasts for all leads in ACC skill but only show an improvement in the CRPS at leads 3-4. DL-DM and DL model exhibits similar skill across all leads which can be the direct result of the model(s) converging on similar weights during the training process, additionally the models within scenarios are averaged across many experiments as well as across CONUS which can diminish the variability. Spatial differences between the DL experiment 24 (EX24) and the DL-DM experiment 29 trained on GEFSv12 reforecast (EX29_G) identify that variability in skill occurs across the entire region, albeit differences may be minimal (Fig. 2C, D). EX29_G mostly shows higher ACC than EX24 except in some of the eastern regions (Fig. 2D). EX29_G also consistently shows better CRPS than EX24 overall of CONUS, suggesting higher ensemble forecast skill (Fig. 2D). While multi-model ensemble is commonly used to increase forecast skill through leveraging the strengths of each model forecast, we identified that this method did not considerably increase ACC skill when combining ECMWF and GEFSv12 (Supplementary Fig. S1).

To identify how input variables influence RZSM forecasts, we conducted a permutation feature importance test on each experimental run (see Methods). We found the most influential predictor for week 3 was GLEAM RZSM reanalysis input at weekly lag 1 (Fig. 3). The ACC values begin to decrease more rapidly after EXP25 (which contained no reanalysis atmospheric predictors). Subsequent experiments after EX25 in Fig. 3 included reanalysis atmospheric predictors and lower skill may be due to a combination of a.) diminished effects of previous atmospheric conditions at longer leads (low signal-to-noise ratio), b.) evolving land-surface conditions (precipitation and other hydrometeorological processes) during the three-week forecasts which are more impactful on RZSM variability or c.) inclusion of too many predictors which inhibits model learning through overfitting. Specifically, research has identified that increasing the number of inputs in a neural network can create a less robust model prone to overfitting and is unable to generalize to out-of-training samples 37.

SHAP values were computed for week lead 1 prediction to further estimate each input feature’s contribution to the final model prediction (see “Methods”). Using GEFSv12 as the forecast, it was observed that the most influential features impacting model outputs included GLEAM RZSM at week lag 1 as the most influential feature, followed by ERA5 t max at week lag 1, GEFSv12 reforecast at week lead 1, GLEAM RZSM at week lag 2, and geopotential height (z 200) from both ERA5 and GEFSv12 (Supplementary Fig. S2). The higher density of SHAP values for GLEAM RZSM week lag 1, which corresponds with negative SHAP values, indicates that wetter soil conditions will lead to a drying of soils in the predicted RZSM. While this may seem counterintuitive, it is probable that the model has learned that saturated conditions lead to future drying. GLEAM RZSM at week lag 1 also identifies that wetter soils can increase the predicted RZSM output, as seen with a slightly lower density of positive SHAP values. ERA5 t max at week lag 1 identifies that lower temperature values lead to an increase in saturation in predicted RZSM. Z 200 results from both ERA5 and GEFSv12 identify that below normal heights, associated with troughs, lead to an increase in RZSM predicted values (towards wetness) as observed with a slightly higher density of values with positive SHAP. These results were identified to be robust across five different random seeds, which details low variability between different SHAP runs (Supplementary Fig. S3).

The best performing model runs included only inputs from observation RZSM and/or recursive RZSM predictions (Fig. 3). The highest skill model at week 3 (EXP29) contains inputs of only lagged RZSM reanalysis and recursive RZSM predictions. While this set of inputs is highly similar to EX18, the recursive RZSM predictions from weeks 1-2 within EX29 come from models which were trained on only week 1 lag of reanalysis atmospheric conditions and short-lead DM atmospheric forecasts (weeks 1-2); therefore, it is likely that this hybrid model EX29 retains the important atmospheric signals embedded within the recursive RZSM predictions in which it is trained (Supplementary Table S1). For clarity, EXP29 will be named as either DL-DM_G (trained on GEFSv12 reforecast and other data) or DL-DM_E (trained on ECMWF reforecast and other data). Hereafter, we refer ECMWF reforecast as ECMWF. It should be noted that predictions and subsequent verification metrics from DL-DM_G and DL-DM_E are highly similar; therefore, we only present one model within each figure.

We evaluate the spatial distribution of skill between DL-DM_G, EMOS, and GEFSv12 (Fig. 4). For CONUS at week lead 3, a high ACC score is achieved with the DL-DM model in many regions (DL-DM_G; ACC > 0.7), with both DL-DM_G being significantly better than GEFSv12 and EMOS in the Central, Northwest, and Southeast US. DL-DM models can attain ACC > 0.6 over 58% of CONUS when compared to EMOS or GEFSv12, having only ~13%. The worse performing region(s) across all models are the Mississippi and Ohio Valley regions. An assessment of seasonal variability displays that months December-February (DJF) contribute most to the overall skill for week leads 3-4 (~0.80 ACC), followed by months June-August (JJA) with ACC between 0.60-0.65, September-November (SON) with ACC between 0.50-0.60, and the lowest average skill is in March-May (MAM) with ACC between 0.47–0.51 (Supplementary Fig. S4). The KGE metric demonstrates that both DL-DM_G and DL-DM_E are better calibrated when compared to raw reforecasts for each season with >70% of grid cells have a positive KGE value as compared ECMWF and GEFSv12 (Supplementary Fig. S5).

ACC skill based on the most recent GLEAM RZSM input across CONUS displays that wetter conditions lead to a more skillful forecast at weekly leads 3-4 in the Western and Northern High Plains regions, while drier conditions lead to more skillful forecasts in the Northwest, Northeast, and Southeast (Supplementary Fig. S6).

We further investigate the model generalizability over the globe and re-train the models for China and Australia (Fig. 5). Again, we find that the highest skill is obtained with DL-DM_G spatial ACC improvements of 23% (18%) compared to the raw reforecasts for China (Australia) respectively for percentage of grid cells with ACC > 0.6. For seasonal variability in China (Supplementary Fig. S7), the highest contribution to overall ACC skill for week leads 3-4 is from SON with average skill across leads (~0.62 ACC), proceeded by DJF (~0.54), MAM (~0.53), and JJA (~0.40). Improvements in skill relative to GEFSv12 display that the Tibetan plateau has no improvements in skill in JJA for weeks 3-4, while the eastern China agricultural region has improved skill during MAM and JJA. Seasonal skill improvements in Australia are primarily contributed from season JJA with averaged week lead 3-4 ACC of ~0.50, proceeded by MAM (~0.40), and similar results between DJF and SON (Supplemental Fig. S8).

We analyzed the continuous ranked probability score (CRPSS) to determine the extent to which the CRPS of RISE-UNet predictions were improved relative to raw forecasts (Supplementary Fig. S9). We identified that over CONUS, the CRPSS of DL-DM_G and DL-DM_E were improved over the entire region, with the Central, Northern High Plains, and Eastern coastal regions showing the highest improvements. In China, DL-DM models have CRPSS > 0 in many regions except for the Central and Northwest, indicating that raw forecasts have a better CRPS in these regions. In Australia, DL-DM models also have CRPSS > 0 in the Western and Eastern portions of the continent, but coastal regions indicate the CRPS of the raw dynamic forecast is better relative to RISE-UNet.

Predictability of dry extremes
For dry weekly RZSM anomalies below the 20th percentile, DL-DM_G is more skillful than raw reforecast based on true positive rate for weekly leads 4-5, false positive rate at leads 3 & 5, and Gilbert skill score at leads 3-4 (Supplementary Fig. S10). The DL-DM_G outperforms GEFSv12 and ECMWF forecasts in the Central and Western United States for the true positive rate, while ECMWF performed well in the Southeast United States, and GEFSv12 performed well only in a few areas in the West (Supplementary Fig. S11). While the models were trained using large samples from 2000–2015, it is worth noting that the samples size was small (n=104) during the testing period (i.e., dry weekly RZSM anomalies below the 20th percentile from 2018-2019). A longer DM reforecast archive can improve the robustness of the evaluations. In addition, when comparing the difference between the predicted and observed (reanalysis) percentiles, DL-DM_G generally shows a smaller percentile difference in the Western United States, albeit a slight dry bias (Supplemental Fig. S12). ECMWF has a dry bias in the Southeast United States and a wet bias in certain locations in the Central and West, and GEFSv12 has both positive and negative percentile differences scattered in the South and West. Raw reforecasts of ECMWF and GEFSv12 and DL-DM_G predictions have a dry bias in the Northwest. Overall, the results indicate that the RISE-UNet generally outperforms dynamic model forecasts for extremes, while the false positives can predominate in certain regions. Over CONUS, we see that the Western United States has much higher true positive rates for anomalies under the 20th percentile, but the Southeast had low performance, which may be caused by shorter SMM in this region 38.

Flash drought case studies
Flash drought is a devastating natural hazard and results in a rapid change in water stress, and can directly impact RZSM. We investigate a series of flash drought case studies beginning with Eastern China, with a flash drought that intensified rapidly during August 2019 and was associated with a strong cyclonic anomaly over the western North Pacific and northeasterly wind anomalies over the northwestern side 39,40. RZSM reanalysis data from GLEAM accurately captured a rapid drying of RZSM anomalies across Southeastern China (Fig. 6). DL-DM_G shows the lowest mean absolute error (MAE), highest spatial correlation (r), and highest structural similarity (SSIM) for all 100% of initializations when compared to GEFSv12. We also provide additional analysis of flash drought case studies over Australia and CONUS, which both occurred in 2019 (Figs. 7 and 8). Australia results indicate DL-DM_G improved 100% of forecasts for SSIM, r, and MAE when compared to GEFSv12. For the flash drought over the Southeast United States, DL-DM_G achieves the lowest MAE in 50% of initializations, the highest correlation (r) in 100%, and the highest SSIM in 67%. However, it still modestly underestimates the spatial extent and intensity of the drought, indicating room for further improvement in extreme event prediction.

Methods
Reanalysis data
We obtained daily soil moisture from December 1999-February 2020 for surface (0–10cm) and rooting-zone (10–100cm) and created a weighted average of surface and rooting-zone for 0–100cm hereafter referred to as root-zone soil moisture (RZSM) from the Global Land Evapotranspiration Amsterdam Model (GLEAM v3.8a) data product. The weighted average was computed as the sum of each variable, which was multiplied by a scaling factor associated with depth with RZSM = surface*10cm + rooting-zone*90cm / 100cm. GLEAM is a set of algorithms which estimates terrestrial evaporation, surface- and root-zone soil moisture, and hydraulic fluxes using satellite observations 41,60. Soil moisture is modeled using both Soil Moisture and Ocean Salinity and European Space Agency Climate Change Initiative satellites and assimilated with Noah Global Land Data Assimilation System reanalysis with data available globally at a spatial resolution of 0.25x0.25°. Data was bilinearly interpolated to a horizontal resolution of 0.5° and a seven-day rolling mean was applied to each grid cell. Seasonal anomalies for each day and grid cell were created by subtracting the seasonal mean computed between 2000–2015 from each observation date, excluding the validation or testing period. Seasons are defined according to conventional categories of DJF (December, January, February), MAM (March, April, May), JJA (June, July, August), SON (September, October, November). Seasonal distributions were created by grouping observations within the same calendar season. The seasonal mean was then subtracted to isolate anomalies.

We retrieved hourly data for 2m dewpoint temperature, 2m maximum temperature (T max), 2m minimum temperature (T min), precipitable water (pwat), surface pressure, and 200 hPa geopotential height (z 200) from the fifth generation of European ReAnalysis (ERA5) between October 1999-February 2020. ERA5 has been shown to be superior to Modern Era Retrospective Analysis version 2 over North America and Europe for convective processes, boundary layer height, temperature, and cloud cover 61–63. We computed the difference between T max and T min, hereafter referred to as Δ tmp, and computed 2m specific humidity (spfh) using dewpoint and surface pressure 64. The rationale for including Δ tmp is that radiation resolving processes within models are less well resolved within low resolution models, which can be impacted by poor cloud representation, which impacts radiation feedbacks 65. Δ tmp considers the diurnal temperature difference and a higher/lower temperature difference over land explicitly implies that there was either less/more cloud cover that day or a frontal storm system occurred which can either increase or decrease Δ tmp. For volumetric soil water, the top 3 layers were weighted and summed by the same procedure as GLEAM RZSM. We bilinearly interpolated all variables to a horizontal resolution of 0.5°, applied a seven-day rolling mean to each grid cell, and constructed seasonal anomalies identically to GLEAM.

For each study region, including the contiguous United States (CONUS), China, and Australia, data was subset according to the specific latitude and longitude bounding box coordinates listed in Supplementary Table S5.

Retrospective forecast (reforecast) from dynamic models
We used subseasonal forecasts of weekly leads 1 to 5 from the Global Ensemble Forecast version 12 (GEFSv12) and subseasonal-to-seasonal European Center for Medium Range Weather Forecasting (S2S ECMWF) models reforecast archive. Within this study, we refer to the S2S ECMWF as ECMWF. The GEFSv12 model was initialized using GEFSv12-reanalysis, and the forecasts were generated by coupling the GFSv15.1 AGCM with Noah LSM to produce a global 11-member ensemble each Wednesday with a lead time of 35 days 66. The first ensemble realization is the control forecast, while the remaining 10 realizations are stochastically perturbed. While the GEFSv12 reforecast data is available for 2020–2022 in a different repository, initializations during this period are produced with a different model version. For consistency, we did not incorporate that period into this study and only considered the period of 2000–2019. The first ten days of the GEFSv12-reforecast are provided on a spatial resolution of 0.25° and temporal steps of three hours, while the last twenty-five days are provided at a spatial resolution of 0.50° and temporal steps of six hours. We downloaded 1043 initializations between the years 2000–2019 for T max, T min, z 200, spfh, and RZSM. For T max and T min, we retrieved only the maximum and minimum values for each day, respectively and computed Δ tmp. For all other variables, we obtained a daily mean by averaging all 3- or 6-hourly leads for each day. For RZSM, we computed the weighted average of soil moisture layer depths 0–10cm, 10–40cm, 40–100cm and obtained the daily mean for the final depth of 0–100cm. RZSM weighting was constructed similarly to GLEAM, where the scaling factor was computed according to the layer depth in cm. For all ECMWF and GEFSv12 predictors, we retain the model using the entire 11-ensemble member forecast. For each lead time, seasonal distributions were created by grouping forecasts whose lead dates fell within the same calendar season. The seasonal mean was then subtracted to isolate anomalies, which are used as model inputs.

The operational ECMWF real-time extended range forecasts are a set of global 51-member ensembles produced twice weekly (Monday and Thursday) with a lead time of 45 days for years 2015-present. For our experiment, we utilized the ECMWF retrospective forecasts (reforecasts), which use the same model version and configuration of the ECMWF real-time forecast but are initialized with past known conditions from ERA5 and are available at 1.5° resolution with a smaller set of 11 ensemble members 67. These 11-member ensembles are produced for the previous 20 years for the same calendar date when the real-time forecast is produced and are hereafter referred to as ECMWF reforecast. ECMWF reforecast is initialized with assimilated observational data, and the predictions are made with the coupled Integrated Forecasting System (IFS) for atmospheric physics and the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL) for land surface dynamics. The first ensemble realization is the control forecast, while the remaining 10 realizations are stochastically perturbed. We subset the data to only include the exact same 1043 initialization dates as GEFSv12 and consider the first 35 days lead time. Within the ECMWF data portal, atmospheric variables provided are not identical to GEFSv12; therefore, we incorporated daily averaged 2m mean temperature (T mean), pwat, 2m dewpoint temperature (T dew), and 0–100cm RZSM. RZSM was conservatively interpolated to 0.5° resolution because bilinear interpolation produces null artifacts along coastlines due to interpolation with missing values associated with oceans. Atmospheric variables were bilinearly interpolated to 0.5° resolution, and we computed anomalies for all variables identically to GEFSv12.

ECMWF and GEFSv12 subseasonal reforecasts have different initialization dates, which can influence a direct comparison between models. To align ECMWF and GEFSv12 reforecasts, we matched lead times based on valid forecast dates, rather than attempting to align initialization dates. For example, an ECMWF reforecast initialized on Monday, Dec 2, 2019 (and its 20-year hindcast ensemble) would be compared against a GEFSv12 reforecast initialized on Dec 2, 2019, even if ECMWF did not produce a forecast on every intervening day. We verified that this approach preserves forecast lead equivalency (e.g., week-3 lead = days 15–21), and we tested sensitivity to potential day-of-week effects, which showed negligible differences in skill scores across the majority of grid cells across CONUS. The process is straightforward and creates a new ECMWF file with a 43-day forecast with the first lead date matching the first lead date from GEFSv12. While short-lead forecasts are sensitive to initialization timing due to the strong influence of initial conditions, at longer leads (e.g., week-4), the marginal degradation in skill from a two-day difference in initialization date is negligible. Therefore, aligning forecast periods across GEFSv12 and ECMWF, even with minor mismatches in start dates, remains a reasonable approximation for comparative skill analysis. While baseline subseasonal forecast skill from dynamic models alone tends to decline significantly beyond two weeks, our results show that hybrid approaches, combining early-lead dynamic forecasts and antecedent soil moisture states within deep learning models, can substantially improve predictability at weeks 3–5.

The selection of GEFSv12 and ECMWF as our two primary models is specifically due to these models providing RZSM data at lower depths (0-1m), which can be directly compared with reanalysis products. This depth has been shown in previous literature to encompass 11 different crops in temperature regions, which are highly relevant within agriculture and can be susceptible to droughts which occur within the subseasonal timescale 13,68. In addition, other models within the SubX project (RSMAS, GMAO, NRL, ECCC) do not provide soil moisture data at unified depths; specifically, these models have a variable rooting depth between 1.5–3m, which can vary between models 16. We did not identify specific models within the S2S project, aside from ECMWF, which may have similar soil moisture depths, although some of these models may be suitable for root-zone drought forecasting and monitoring 67.

Data pre-processing
Feature (variable) scaling is necessary for statistical and DL models because features with different scales decrease model flexibility by forcing the model to simply reduce the error from the feature with the highest scale difference relative to the predicted output. This reduction in error does not imply an increase in testing performance because that feature may not be highly associated with the final prediction. To assist with model setup and standardization procedure, each observational variable was reformatted from a 3-dimensional (time, latitude, longitude) time series to a 5-dimensional (initialized date, model, lead, latitude, longitude) time series that matches the same lead-time format as GEFSv12 and ECMWF initializations. We split the data into training (years 2000–2015), validation (2016-2017), and testing (2018-2019). For each variable, seasonal anomalies were scaled using min-max normalization, where the minimum and maximum values were obtained over the entire spatial domain from only the training period, with the following equation x' ijk = [x ijk - min(x i)] / [max(x i) - min(x i)]. With x ijk and x' ijk being the unscaled and scaled feature value, respectively, for lead i, latitude j, and longitude k. All values will be within the range 0-1. For RZSM, ocean and water bodies have no values and are masked with a value of zero prior to training the model.

To prepare input for the RISE-UNet models (described in the next section), we stack the observations and reforecast predictors (see Supplementary Table S1) according to the initialized date and the ensemble realization to produce a final input shape during training of (9185 samples, 48 latitude grid cells, 96 longitude grid cells, N channels). For a more in-depth understanding of the creation of lagged observations, see Supplementary Table S6. The 9185 samples were obtained by training and initializing into a lower-dimensional latent feature space (835 initializations x 11 ensemble realizations = 9185 total samples). N channels is dependent upon the experiment number found in Supplementary Table S1. The verification and testing data were stacked similarly with sizes (1144 samples, 48 latitude grid cells, 96 longitude grid cells, N channels).

Deep learning architecture and experiment setup
Model architecture. A climate variable at each time step can be viewed as a two-dimensional image. Convolutional neural networks (CNNs) are powerful for learning complex functional relationships across spatial dimensions and covariate features among different images. The UNet architecture based on CNN 69 is a deep neural network model containing a symmetric encoding-decoding structure that learns both low- and high-level feature abstractions between covariates and across spatial regions. Previous studies have used UNet architectures to bias correct and post-process both short- and long-lead NWP ensemble forecasts of temperature, geopotential height, and precipitation, with results indicating better performance than raw forecasts for various leads and regions 36,46,70. A new variant of the UNet is the UNet++ 3, which contains a series of nested and dense skip connections connecting each encoder-decoder block with previous blocks across the same spatial dimensions [Fig. 1B] 71. This nested architecture has benefits over traditional UNet, which can suffer from information loss as downsampling occurs and is overcome through multiple nested skip connections. We modified the original UNet++ 3 architecture by adding inception, residual, and squeeze-and-excitation 72 (SE) blocks to increase the model’s spatial relation understanding and feature space representation. The new model architecture is named as residual inception squeeze-and-excitation (RISE) UNet++, hereafter referred to as RISE-UNet (Fig. 1). For each Inception operator block, we include four parallel convolutional filters using kernel sizes of 3×3, 5×5, 7×7, and max pooling layer 5×5 and filter sizes of [32, 64, 128, 256] for the 1st through 4th encoder backbone respectively. After each convolution, batch normalization (BN) is applied, followed by the rectified linear unit (ReLU) activation, and we pass another 1×1 convolution kernel with the same number of previous filters, followed by BN. We concatenate the four parallel convolution blocks, add the original input as a residual connection to reduce diminishing gradients and finally apply a SE block with sigmoid activation function to the summed residual/inception block in order to weight inter-channel features. Dropout layers are applied within each convolutional block to act as a regularization technique to prevent overfitting and improve generalization. Between 10–25% of the two-dimensional layers are randomly dropped during the training process, which creates stochastic predictions during model inference over the testing period. We train the model individually (and recursively for selected experiments) for weekly leads 1–5. Additional information on DL model training and hyperparameters can be found in Supplemental Table S3.

Cost function and ensemble generation process. We used the experimental Continuous Ranked Probability Score 73 (CRPS exp) as a loss function, which computes the mean absolute error (MAE) and subtracts a scaled standard deviation (f*σ). Where, CRPS exp = MAE - (f*σ) (1). MAE = 1 / j*k * Σ n i=1 |^y i - y i| (2). σ = sqrt( 1 / j*k * 1 / 11-1 * Σ 11 r=1 (^y r - ^y i) 2) (3). The scaling factor f is 0.08. σ is the standard deviation, j is the number of latitude grid cells, k is the number of longitude grid cells, ^y is the predicted value, y is the observation, i is the 11-member ensemble prediction, and r is the predicted realization. Higher standard deviations of the forecast will decrease the overall loss values, and we use a batch size of 66 to account for six initializations with eleven realizations and CRPS exp can be computed on each eleven-member realization to evaluate the standard deviation. Although the loss function is not an explicit calculation of CRPS, the CRPS exp computes the MAE on each of the eleven realizations independently and also computes the standard deviation across the forecast initialization. For clarification, the MAE and standard deviation calculation is computed on each batch, and the process begins with partitioning each batch size of 66 into 6 groups of 11-members. On each group, the MAE is calculated on each ensemble prediction and the standard deviation is computed across the 11-member ensemble. This set of functions allows the model to adjust model weights, which both decrease the error and increase the spread of the forecast, thereby retaining stochasticity similar to that of physically-based dynamic models.

Weather and climate ensemble forecast systems employ various perturbation methods to account for uncertainty in initial conditions and allow for the quantification of uncertainties due to errors in initial conditions, model formulation, and physical parameterizations 74–77. More specifically, when using a singular realization contrived from assimilated observations and processed through a physically-based land and atmospheric model, uncertainty can be imparted from instrument calibration error or model structural errors associated with parameterizations of complex processes. Small errors within initial conditions have been observed to rapidly amplify in weather forecasting and restrict the intrinsic predictability limit of forecast systems 78; therefore, perturbation methods allow for inclusion of a wider range of plausible initial conditions to evaluate uncertainty in ensemble forecasts through sensitivity experiments and probabilistic evaluations. The RISE-UNet DL architecture accepts an ensemble of stochastically generated initial conditions and passes these inputs through a series of operators (convolution kernels, squeeze-and-excitation, and max-pooling), non-linear activation functions (ReLU and sigmoid), and dropout layers to produce stochastic ensemble forecasts for each lead week. The CRPS exp loss function is computed over each set of 11 ensemble members generated from the RISE-UNet, which allows for gradient-based model weight adjustments based on the aggregated loss.

Experimental setup. To determine the optimal model and understand what impacts the forecast skill, we constructed a series of experiments (Supplemental Table S1) in which the RISE-UNet DL model only takes historical reanalysis (DL), only dynamic model (DM) reforecast (DM-BC_DL), or combined reanalysis and dynamic model (DM) reforecasts (DL-DM) as inputs to predict weekly RZSM anomalies for leads 1–5. For the DL and DL-DM scenarios, a forecast is made by recursively applying RISE-UNet to the previous prediction to produce a prediction at successive lead weeks 2–5 (N). For further clarification, weekly lag/lead inputs have been computed by the following pre-processing pipeline: apply a 7-day rolling mean, compute seasonal anomaly, and perform min-max standardization. For example, Week lag 1 of reanalysis represents the antecedent information from the day before the forecast initialization date through the previous week. Specifically, the scenario DM-BC_DL includes experiments EX0 and EX13, which bias-correct DM forecasts using reanalysis, where EX0 uses all previous reforecast weeks as input predictors, while EX13 only includes the current week (single predictor). The scenario DL includes only lagged reanalysis as inputs, with EX14-EX17 containing atmospheric variables and RZSM, while EX22-EX25 includes only RZSM. DL-DM is the combined model including both lagged reanalysis and DM reforecasts, which includes EX1-EX4, EX5-EX8, EX9-EX12, EX27-EX28, and EX29. EX1-EX4 includes atmosphere and RZSM from reanalysis, as well as recursively predicted RZSM inputs from the previous time step. Recursively predicted inputs are the predictions of the previous lead weeks from the DL-DM model, which are used to predict the next lead week. EX9-EX12 is identical to EX1-EX4 with the exception of including DM RZSM for week (N) as an input. EX5-EX8 includes only RZSM reanalysis and DM RZSM for week (N) as inputs, and EX18-EX21 further include recursively predicted RZSM. EX27-EX28 are identical to EX9-EX10, with the exception of including DM reforecast atmospheric inputs up to lead week two. The setup of experiments in this manner allows for a robust assessment of each set of inputs on forecast skill 79.

Feature importance evaluation
Permutation test. To determine the impact of each input variable on model skill in the testing dataset (2018-2019), we conducted a permutation test on each input individually for each experimental run. The order of operations for the permutation test for each experimental run goes as follows: First, a single input variable was selected (e.g., RZSM observation week lag 1) and random Gaussian noise was added, with mean and standard deviation from training data input statistics. Second, the input data was transformed via min-max normalization based on the new minimum and maximum. Third, with only a single input transformed and all other inputs identical, the combined inputs were fed into the RISE-UNet pre-trained model, and 11 ensemble predictions were made for each initialization in the testing period. Finally, the ensemble mean was computed and MAE evaluated.

SHapley additive exPlanations (SHAP). To further interpret the predictions of our deep learning model, we employed SHAP, which is a unified framework that attributes a prediction to input features based on cooperative game theory 80,81. SHAP assigns each feature an importance value for a particular prediction by estimating the contribution of each feature to the model output. SHAP values are grounded in cooperative game theory and satisfy the property of additive feature attribution, f(x) = φ 0 + Σ M i=1 φ i (4) where f(x) is the model prediction for input x, φ 0 is the expected model output, and φ i denotes the SHAP value corresponding to the input feature i. To efficiently compute SHAP values, we used the Gradient Explainer, which estimates SHAP values based on the model’s gradients with respect to its inputs, integrated along a straight-line path from a baseline input to the actual input and computed as, φ i = (x i - x' i) ∫ 1 α=0 ∂f(x' i + α(x i - x' i)) / ∂x i dα (5) where x' is a baseline input and ∂f / ∂x i is the gradient of the model output with respect to the input feature x i. This formulation captures how the model’s prediction changes as the input transitions from the baseline x to the actual value x i. Interpreting these results involves examining both the magnitude and sign of each φ i with larger absolute values denoting greater influence on the final prediction, while the sign indicates whether the feature pushes the prediction higher (positive φ i) or lower (negative φ i) relative to the baseline prediction φ 0. Therefore, SHAP Gradient Explainer shows the relative influence of each predictor on the final output and identifies how small or large changes in the feature influence the final output, mimicking the methods of a sensitivity test. For computation, we randomly selected 500 background samples from the training dataset to evaluate SHAP values over the testing dataset.

Ensemble model outputs statistics (EMOS)
As a comparison, DM is also bias-corrected using ensemble model output statistics (DM-BC_EMOS). EMOS is a parametric statistical post-processing technique which incorporates NWP ensemble predictions and fits a predictive model to minimize the continuous ranked probability score (CRPS) 82. For every observation, EMOS conducts multiple linear regression to compute a predictive probability density function (PDF). Let X 1, …, X n indicate N ensemble forecasts for a univariate continuous, positive-defined predicted variable Y. Then Y = a + b 1X 1 + ... + b nX n + ε where the a, b 1, ..., b n are the regression coefficients and ε is the error term which averages to zero. The variance of ε is modeled as a linear function of the ensemble spread, Var(ε) = c + dS 2 where S 2 is the ensemble variance, and c and d are non-negative coefficients. These terms yield the Gaussian predictive distribution N(a + b 1X 1 + ... + b nX n, c + dS 2) whose mean forms a bias-corrected weighted average of the ensemble forecasts, and whose variance depends linearly on the ensemble variance. With the predictive PDF, ensemble predictions can be made of any size.

For each weekly lead and grid cell, an EMOS model was trained on the 11-member GEFSv12 anomaly RZSM ensemble. We evaluated four different EMOS modeling schemes including training on data from: a) only dates that were 12 weeks before the predicted date with data from only the training period (2000–2015) [192 initialization dates], b) only dates +/− 12 weeks from the predicted date in the training period [384 initialization dates], c) only dates that were 12 weeks before the predicted date with data from only the training and validation period (2000-2017) [216 initialization dates], d) only dates +/− 12 weeks from the predicted date in the training and validation period [432 initialization dates]. The purpose of these four different schemes allows for a more robust metric comparison from models trained on different sample sizes (deep learning training = 835 initialization dates vs. EMOS training = ~192-432 initialization dates), and allows for an EMOS inter-model comparison. Training EMOS on all 835 initializations produced highly unskillful predictions (results not displayed). Models were evaluated over the testing period (2018-2019). More information on each EMOS experiment can be found in Supplementary Table S4).

Evaluation metrics
For all evaluation metrics, forecast skill results are averaged over forecasts made at weekly intervals across the entire calendar year for the two-year testing period 2018-2019. We evaluate probabilistic forecast skill with the CRPS and continuous ranked probability skill score (CRPSS), and deterministic skill with the Anomaly Correlation Coefficient (ACC). The CRPS measures the performance ensemble forecasts relative to observations by computing the forecast predictive cumulative distribution function (CDF) calculating the integral of the squared difference between the CDF and the observation, which is considered a step function. CRPS is concerned with the accuracy and spread of the forecasts, rather than an individual point estimate, with lower CRPS values indicate greater skill and a value of 0 indicates a perfect prediction model. CRPS is defined as follows, CRPS(F, y) = ∫ -∞ ∞ (F(x) - 1fx>y) 2 dx where F(x) is the predicted CDF, 1 is the Heaviside step function, where values of observation (x) that are greater than prediction (y) are given a value of 1 and y < x is given a value of 0. The CRPS was computed for each grid cell and averaged overall initialization dates within the testing period (2018-2019).

The CRPSS compares the CRPS of the ensemble forecast to that of a reference forecast to determine whether the ensemble spread is a useful measure of forecasting uncertainty and this metric was computed with the Python package 83. The CRPSS is formally defined as CRPSS = 1 - (CRPS RISE-UNet / CRPS raw) where values greater than 0 indicate that the RISE-UNet forecast is more skillful than the raw forecast, and values less than 0 indicate RISE-UNet is less skillful, with the range of CRPSS being (-1, 1].

The ACC describes the correlation between a forecast anomaly relative to forecast climatology and a verifying observation anomaly relative to observation climatology. ACC indicates the strength of association for both direction and magnitude of forecasts relative to observations, with values falling between −1.0 and 1.0. ACC = Σ M m=1 (f m - f 1) (o m - o 1) / [Σ M m=1 (f m - f 1) 2 Σ M m=1 (o m - o 1) 2] 1/2 where M is the total number of time samples, f is the forecast, f 1 is the forecast climatology, o is the observation, and o 1 is the observation climatology. ACC approaching 1.0 indicate strong positive agreement where forecasts have value, ACC values > 0.6 indicates synoptic prediction skill for large weather patterns, ACC values approaching 0.5 indicate the forecast errors are similar to errors from a forecast of the climatological average, and ACC values around 0 indicate the forecast has no value 84. We calculated the ACC using the mean of all ensemble members within the testing period 2018-2019 and averaged the results by region.

Categorical (binary) accuracy was assessed for anomalies falling below the 20th percentile using the True positive rate (TPR), False positive rate (FPR), and Gilbert skill score (GSS). Forecasts were evaluated using a contingency table in which each forecast lead time is verified against observations from the same date where the random variable true positive (TP) indicates that forecast and observation both classified anomalies being below the 20th percentile, false positive (FP) when the forecast predicted anomalies below the 20th percentile but the observation did not, false negative (FN) when the forecast predicted anomalies above the 20th percentile when observations were above the 20th percentile 13. The true positive rate (TPR), also known as sensitivity, assesses the percentage of correct forecast predictions compared to all predicted events, including events that did not occur in observations 85. TPR = Σ TP jk / Σ (TP jk + FN jk) where j, k are the respective grid cell and TPR values range from 0–100% where 100% indicates that all anomalies below the 20th percentile were correctly classified. We calculated TPR for each grid cell and averaged the results by region. The false positive rate (FPR) 85 assesses the percentage of predictions that were incorrectly classified as positive and is the ratio of incorrectly forecasted events below the 20th percentile to the total number of actual negative events below the 20th percentile FPR = Σ FP jk / Σ (FP jk + TN jk). The Gilbert skill score (GSS) 85, also known as the Equitable threat score, measure the fraction of events that were correctly predicted and is adjusted for hits associated with random chance and equitability allows scores to be compared more fairly across different regimes GSS = (Σ TP jk - Σ TP random) / (Σ TP jk + Σ FN jk + Σ FP jk - Σ TP random). TP random is the number of TP hits by random chance and is defined as TP random = (Σ TP jk + Σ FP jk) * (Σ TP jk + Σ FN jk) / N where N is the total number of observation events. Values for the GSS fall within the range −1/3 to +1 with values of 0 indicating no skill and values of 1 indicating perfect skill.

The Structural Similarity Index (SSIM) is commonly used to measure the similarity between a set of images and provides information on three image components, including luminance, contrast, and structure 86. Luminance (μ x) is the average of all pixel values between nearby pixels, and structure (s x) is the unit-normalized image values μ x = 1/N Σ i=1 N x i, σ x = sqrt(1/N-1 Σ i=1 N (x i - μ x) 2), s x = (x - μ x) / σ x. For each of these components, comparison functions are applied to the two images and SSIM evaluated. SSIM values lie within the range -1.0 and 1.0, with values of 1.0 indicating identically similar images, values of 0 indicating different images, and values approaching −1.0 indicating complete dissimilarity. Since SSIM is based on visible structures in the image, it measures the perceptual difference between two similar images, but not which one has higher accuracy.

The Kling-Gupta efficiency (KGE) 87 is an alternative to the Nash-Sutcliffe Efficiency (NSE) and is designed to address some of its limitations, particularly in balancing the bias, variability, and correlation of the simulated data with observed data. KGE = 1 - sqrt((r-1) 2 + (β-1) 2 + (γ-1) 2). It ranges from -∞ to 1 and is computed as the sum of three components, including correlation (r), bias ratio (β), and variability ratio (γ). Pearson linear correlation (r) measure the linear correlation between simulated and observed samples, β represents the ratio of the means of the simulated and observed data β = μ simulated / μ observed, and γ measure the ratio of the standard deviations of the simulated and observed data γ = σ simulated / σ observed. KGE values approaching 1 indicate near-perfect agreement between the simulated and observed data and values approaching -∞ indicate very poor model performance.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript presents an evaluation of dynamical model (DM) and deep learning (DL) hybrid subseasonal forecasts (DL-DM) for predicting Root Zone Soil Moisture (RZSM) up to five weeks ahead. The authors demonstrate that DL-DM forecasts outperform reanalysis-driven DL models, and both approaches significantly surpass DM alone and ensemble model output statistics. Additionally, DL-DM exhibits higher skill in predicting flash droughts compared to DM alone. This study offers valuable insights into the integration of dynamical models and deep learning, highlighting potential improvements in soil moisture and SM-based flash drought predictions. However, I have several concerns regarding data preprocessing, model selection, and case study selection that need to be addressed.\n1. The manuscript compares the RZSM forecasting skills of DL-DM and DM approaches. Could the authors clarify why they opted to use DM-simulated RZSM data instead of soil moisture (SM) data derived from established hydrological models such as VIC or Noah-MP?\n2. Have the authors conducted any evaluation or validation of the reanalysis (observation) data used in the study? Assessing the quality and reliability of the reanalysis datasets is crucial for interpreting the forecast skills and ensuring the robustness of the results.\n3. The study utilizes GEFSv12 and ECMWF models. Could the authors elaborate on the criteria or rationale behind selecting these two specific models over others? Have the authors applied any bias correction techniques to the reanalysis and reforecast data from GEFSv12 and ECMWF? Implementing bias correction is recommended to mitigate systematic errors in the forecast data, thereby enhancing the accuracy and reliability of the predictions.\n4. In Figure 2, DL-DM and DL models exhibit nearly identical skills in both Anomaly Correlation Coefficient (ACC) and Continuous Ranked Probability Score (CRPS). Could the authors explain why there is no significant difference in performance between DL-DM and DL in this context? It appears that the number of runs varies between different models or methods. What was the rationale behind selecting these specific numbers of runs for each model? Additionally, could the differing numbers of runs influence the comparative results and conclusions drawn from this analysis?\n5. Figures 4&5 : To provide a clearer understanding of the improvements achieved by the DL-DM approach, it would be beneficial to include additional figures that depict the differences in skill metrics between DL and DL-DM models. For instance, presenting a figure that shows the skill of subplot (C) minus the skill of subplot (E) would effectively highlight the enhancements brought about by the hybrid approach.\n6. The manuscript examines flash droughts using only a single case study from China. To enhance the robustness and generalizability of the findings, it is recommended that the authors include multiple flash drought events from diverse geographical regions and climatic conditions. Incorporating a broader range of case studies would provide more comprehensive evidence of the method’s effectiveness and applicability across different environmental contexts.\nMinor comments:\n1. In the introduction (Lines 73-75), the authors mention various methods used in the study. It would enhance clarity if the authors could provide more detailed explanations of the differences between these methods, particularly distinguishing between the hybrid DL-DM approach and DL-processed DMs. This additional information would help readers better understand the unique contributions of each method.\n2. Fig 3. I would recommend adding color legends for the runs.\n3. Line 99: or would it be because of the poor atmospheric prediction skills?\n4. lines 101-105, I'm a little confused about this sentence, could the authors make it clearer?\n5. On lines 114-115, the authors note that the Mississippi and Ohio Valley regions exhibit poorer performance. Are there specific factors or characteristics of these regions that contribute to this reduced performance?\n6. In the Results section, it is unclear what reference data was used for evaluating the forecasts. Although descriptions are provided in the Discussion, it would be beneficial to include a brief mention of the reference data either in the Results section or the Introduction.", "reviewer_2": "This manuscript examines using deep learning to predict soil moisture and drought conditions beyond the typical two-week forecast horizon. The proposed approach combines historical reanalysis data and dynamic model forecasts to predict weekly root-zone soil moisture anomalies up to five weeks ahead, benchmarked against dynamical forecasting systems and deep learning post-processed predictions.\nMajor Concerns:\n1. Fundamental Issues with Problem Formulation\n- The study does not address the operational reality that current soil moisture status is not immediately available, due to inherent delays in gathering observations and running analysis models for both soil moisture and atmospheric states\n- The methodology for generating ensemble predictions lacks clarity, particularly regarding how ensemble spread is achieved without apparent stochastic elements, which is crucial for CRPS optimization\n2. Limited Scientific Advancement\n- The key findings largely confirm known behaviors: soil moisture's relatively long memory compared to atmospheric processes, and the degradation of atmospheric prediction skill beyond two weeks\n- The study fails to provide meaningful process-based analysis of how early precipitation forecasts contribute to sustained soil moisture predictability\n- The proposed model lacks mechanistic interpretability that could advance our understanding of soil moisture dynamics\n3. Scientific Merit\nThis work appears to be a straightforward machine learning application study without sufficient novelty or significance for Nature Communications publication standards.\nTechnical Issues:\nAbstract and Introduction:\n- Imprecise terminology (e.g., \"Deep learning with historical reanalysis data\" should be \"Deep neural networks that learn from climate reanalysis data\")\n- Inconsistent logic in comparing deep learning + dynamical model performance versus deep learning + reanalysis\n- Confusion between organizational names and model references (e.g., ECMWF)\n- Inadequate definition of key terms and metrics\n- Operational feasibility concerns regarding the use of lagged reanalysis data\n- Mischaracterization of error growth mechanisms in dynamical systems\nMethodology:\n- Insufficient description of model architecture and problem setup prior to results presentation\n- Unclear definitions of key variables and temporal relationships\n- Imprecise statistical framework for extreme event analysis\nThe manuscript requires substantial revision to address these fundamental issues before it could be considered for publication.", "reviewer_3": "Key results\nThe paper demonstrates substantially improved soil moisture prediction skill on subseasonal time scales (1-5 weeks ahead) over the contiguous US, Australia and China, obtained using novel deep learning (DL) compared to the raw output from state-of-art dynamical prediction models (DMs), as well DMs calibrated with conventional linear statistics. The improved skill is found to be largely attributable to soil moisture memory. Purely empirical DL forecasts are found to perform fairly similarly to a hybrid method that combines DL with DMs, although the hybrid approach obtains the highest skill.\nValidity\nThe technical quality appears high and generally convincing. The main findings are shown to be robust across three different geographical regions. In the main MS evaluation is largely in terms of deterministic skill, using anomaly correlation, and over a short 2-year evaluation period, though other skill metrics are available in the SI.\nOne apparent technical flaw needs to be addressed: The GEFSv12 and ECMWF dynamical model reforecasts (years 2000-2019) use different configurations which means that the same forecast start dates are not the same between both both (L274). The main findings don’t depend on this, but it needs to be addressed because it (rightly or wrongly) indicates an incomplete understanding of subseasonal forecasting that otherwise casts doubt over the study.\nSignificance\nWhile the findings suggest a substantial improvement in soil moisture prediction skill obtained using DL, there is not enough comparison against simple benchmarks (persistence and climatological forecasts), or with the results of previous work, for the significance of the reported DL skill increases to be properly assessed. The very short length of the testing period (only 2 years) is a limitation that needs to be acknowledged.\nThe ACC values obtained here should be compared to previous studies in the published literature on seasonal and subseasonal forecasting of soil moisture. How do they compare in general with previous studies on both the seasonal and subseasonal time scale? The current results will not mean much unless they are set in context of previous work.\nRelatedly, a general weakness is the lack of physical interpretation of the results. In particular, the results are not stratified by calendar season, which is important for interpretation, and is standard practice. The testing period is only 2 years, but I think such a seasonal breakdown ought to be included in the SI.\nThe main skill assessment is in terms of deterministic metrics which makes the significance of the main findings hard to judge. It is widely accepted that forecasts need to be expressed probabilistically in order for them to be usable in societal decision making.\nData and methodology\nOverall technical quality appears high, with the caveats mentioned in my previous remarks. However, the labelling of the many experiments is highly confusing, and the selection of results presented in the figures often obscures simple interpretations. Why are different lead times (week 3 and 4) shown in Figs 3-6? Since weeks 3+4 is established practice in subseasonal forecasting, it would be simplest to choose that, or at least present all figures for either week 3 or week 4.\nSecond, the key question to answer is how much this new DL approach can improve on established practice with dynamical models, which is well represented by DM-EMOS. The ACC and MAE for each raw DM (ie without any bias correction) and the DM-EMOS should be included in Fig 3. Assuming the DM-EMOS outperforms DM-raw, then these should be the benchmarks for the case studies in Figs 6, S5, S6, not the raw DM.\nThe captions in Figs. 6, S5, S6 need to be include the dates of there target periods, not just the initialization dates. It should be clarified that GLEAM refers to the validation data.\nAnalytical approach\nWhy was the CRPS used, as opposed to CRPSS (L395)? The latter is a much more common metric in meteorological and climate studies because it is a Skill Score that expresses the forecast accuracy relative to a benchmark, typically the climatological forecast. Unless the authors can justify the choice of CRPS, and provide a reference where CRPS has been used in the meteorological field, then CRPSS should be used instead of CRPS.\nIt needs to be clarified that the results are not stratified by calendar season (L475). While year-round results may still be informative, there may be a contribution from mean seasonality that can inflate the skill estimates, depending on the details of how the seasonal cycle has been subtracted. This needs to acknowledged and addressed.\nMulti-model ensembling of DMs boost skill and is established practice in seasonal and subseasonal forecasting. Why is that not done here?\nSuggested improvements\nL48: Give “beyond week 2” in days or just say “beyond two weeks” and clearly define the “subseasonal” forecast range in days and weeks, here at the outset\nL66 & 67: explain “shorter” hirer in terms of weeks, and “medium range” for the Nature audience\nL84/L88 : Clearly define “weekly” A table with daily equivalents of the weeks may be useful. Is week 1 comprised days 1-7 for example. These definitions vary in subseasonal studies.\nL83: Explain near here that the skill results are averaged over forecasts made at weekly intervals across the entire calendar year, for the two-year period 2018-19. Those are important for there interpretation of the results and should be highlighted here in addition to the methods section.\nL108: Spatial *distribution” of skill\nL189: The statement “Even after DL postprocessing, dynamic models GEFSv12 and ECMWF show no skill for predicting RZSM beyond week 2” is at odds with Figs 2, 4 & 5. An ACC threshold of 0.5 is only a rule of thumb for skillful forecasts.\nL237: Please describe how the “seasonal mean” is defined, with respect to a particular weekly forecast target period. Are these rolling 3-month averages?\nL247: Clarify that Delta_tmp is the diurnal temperature difference. Note that passage of a frontal system could reduce this, as well as increase it.\nL287: Please state the lat-lon domains used, and clarify the lead time definition of the forecast weeks, eg using a table.\nL299: Clarify the meaning of “channels”\nL656: This Santoalla reference link is broken.\nACC Figures: The color scales need to be standardized across the figures, using a simple color scale, eg -0.4 to +0.9\nClarity and context, References\nThe manuscript is well written, generally clear and set well in context of past work and references previous literature appropriately.\nThe study of Lesinger et al (2024) is highly relevant and ought to be compared against more closely. The authors do cite it.\nL312: UNet has previously been used by Horat and Lerch (2024, MWR) for S2S forecasting, which ought to be cited and commented on, including the key difference between their model and the RISE-UNet.\nI did not see a link to the codes used in the MS. Including one would be very beneficial to researchers."}, "author_response": "Response to the Reviewers for the Manuscript Entitled “Skillful Subseasonal Soil Moisture Drought Forecasts with Deep Learning-Dynamic Models” (Ref.: NCOMMS-24-67423A)\nWe genuinely thank the editor and Reviewers for reviewing our manuscript “Skillful Subseasonal Soil Moisture Drought Forecasts with Deep Learning-Dynamic Models”. We also appreciate your time and effort for providing additional constructive feedback to further improve the manuscript and increase readability. We have carefully considered and incorporated all your suggestions into the manuscript. Please see below, in blue, for a point-by-point response to your comments. Bolded blue line numbers represent the new line numbers in the revised manuscript (clean version) when “No Markup” is selected in the “Review” tab within Word.\nReviewer #1:\nReviewer #1 (Remarks to the Author):\nThe manuscript presents an evaluation of dynamical model (DM) and deep learning (DL) hybrid subseasonal forecasts (DL-DM) for predicting Root Zone Soil Moisture (RZSM) up to five weeks ahead. The authors demonstrate that DL-DM forecasts outperform reanalysis-driven DL models, and both approaches significantly surpass DM alone and ensemble model output statistics. Additionally, DL-DM exhibits higher skill in predicting flash droughts compared to DM alone. This study offers valuable insights into the integration of dynamical models and deep learning, highlighting potential improvements in soil moisture and SM-based flash drought predictions. However, I have several concerns regarding data preprocessing, model selection, and case study selection that need to be addressed.\n1. The manuscript compares the RZSM forecasting skills of DL-DM and DM approaches. Could the authors clarify why they opted to use DM-simulated RZSM data instead of soil moisture (SM) data derived from established hydrological models such as VIC or Noah-MP?\nThanks for your comments. In this study, we used DM-simulated RZSM data instead of soil moisture (SM) data derived from established hydrological models because directly bias correcting RZSM forecasts from dynamic models may be more efficient and accurate than using land surface hydrological models forced by bias corrected inputs from dynamic models. A previous study has shown that dynamic models such as GEFSv12 can show low subseasonal soil moisture forecast bias in some specific regions and seasons, suggesting its potential value (Lesinger et al. 2024). By contrast, soil moisture forecasts produced by land surface hydrological models such as VIC or Noah-MP have proven ineffective beyond two weeks in predicting subseasonal drought even after dynamic downscaling, bias correction, and being forced with the ensemble mean of 6 subseasonal models within the SubX project (Su et al. 2023). This was due to high precipitation errors after week 2 which negatively influence soil moisture predictability. We agree that comparing DM-simulated RZSM data and soil moisture (SM) data derived from established hydrological models could deserve a separate study. To address your comments, we have made following changes to the manuscript.\n[L 48-52] “Previous studies have consistently identified limited subseasonal forecast skill beyond two weeks (>14 days) in current dynamic forecast systems and land surface models: including but not limited to Noah-MP, Weather Research and Forecasting (WRF), Variable Infiltration Capacity (VIC), Global Ensemble Forecast System (GEFS), and subseasonal-to-seasonal European Center for Medium Range Weather (S2S ECMWF).”\n2. Have the authors conducted any evaluation or validation of the reanalysis (observation) data used in the study? Assessing the quality and reliability of the reanalysis datasets is crucial for interpreting the forecast skills and ensuring the robustness of the results.\nThank you for your comments. We have now provideda discussion on the accuracy of the GLEAM data in the revised manuscript.\n[L 240-247] “We acknowledge that GLEAM may have its own systematic biases with respect to in-situ observations and improving the quality of the ground-truth reanalysis can be informative in future research to accurately predict rapid changes in RZSM and drought. Although we have not directly evaluated GLEAM soil moisture within our study, previous research has conducted this evaluation with in-situ observations finding a global average correlation of 0.64 (Martens et al. 2017). Research in the Tibetan plateau identified an R value of ~0.8 between in-situ measurements and GLEAM which was higher than GLDAS (Yang et al. 2022), and further evidence of a strong correlation between GLDAS and GLEAM has been identified regarding the direction of soil moisture trends (Guan et al. 2023).”\nReferences:\nMartens, B. et al. GLEAM v3: Satellite-based land evaporation and root-zone soil moisture. Geosci Model Dev 10, 1903–1925 (2017).\nYang, S., Zeng, J., Fan, W. & Cui, Y. Evaluating Root-Zone Soil Moisture Products from GLEAM, GLDAS, and ERA5 Based on In Situ Observations and Triple Collocation Method over the Tibetan Plateau. J Hydrometeorol 23, 1861–1878 (2022).\nGuan, Y. et al. Spatio-temporal variations in global surface soil moisture based on multiple datasets: Intercomparison and climate drivers. J Hydrol (Amst) 625, 130095 (2023).\n3. The study utilizes GEFSv12 and ECMWF models. Could the authors elaborate on the criteria or rationale behind selecting these two specific models over others?\nThank you for your comments. We have now provided an additional paragraph in the data & methodology section to discuss our chosen models.\n[L 439-448] “The selection of GEFSv12 and ECMWF as our two primary models is specifically due to these models providing RZSM data at lower depths (0-1m) which can be directly compared with reanalysis products of the same depth. This depth has been shown in previous literature to encompass 11 different crops in temperature regions which are highly relevant within agriculture and may be susceptible to droughts occurring within the subseasonal timescale (Fan et al. 2016; Lesinger et al. 2024). Additionally, other models within the SubX project (RSMAS, GMAO, NRL, ECCC) do not provide soil moisture data at unified depths and therefore cannot be applied; specifically, these models have a variable rooting depth between 1.5-3m which can vary between models (Pegion et al. 2019). We did not identify specific models within the S2S project, aside from ECMWF, which may have similar soil moisture depths due to data processing and time constraints, although some of these models may be suitable for root-zone drought forecasting and monitoring (Vitart et al. 2017).\nReferences:\nFan, J., B. McConkey, H. Wang, and H. Janzen, 2016: Root distribution by depth for temperate agricultural crops. Field Crops Res., 189, 68–74,\nhttps://doi.org/10.1016/j.fcr.2016.02.013.\nLesinger, K., Tian, D. & Wang, H. Subseasonal Forecast Skill of Evaporative Demand, Soil Moisture, and Flash Drought Onset from Two Dynamic Models over the Contiguous United States. J Hydrometeorol 25, 965–990 (2024).\nPegion, K., and Coauthors, 2019: The Subseasonal Experiment (SubX): A Multimodel Subseasonal Prediction Experiment. Bull. Amer. Meteor. Soc., 100, 2043–2060, https://doi.org/10.1175/BAMS-D-18-0270.1.\nVitart, F., and Coauthors, 2017: The Subseasonal to Seasonal (S2S) Prediction Project Database. Bull. Amer. Meteor. Soc., 98, 163–173, https://doi.org/10.1175/BAMS-D-16-0017.1.\nHave the authors applied any bias correction techniques to the reanalysis and reforecast data from GEFSv12 and ECMWF? Implementing bias correction is recommended to mitigate systematic errors in the forecast data, thereby enhancing the accuracy and reliability of the predictions.\nWe applied EMOS (the scenario DM-BC_EMOS) and DL (the scenario DM-BC_DL) to bias correct raw forecasts from GEFSv12 or ECMWF using reanalysis. The results indicate the reanalysis-driven DL model and hybrid DL-DM model consistently show better performance compared to the EMOS or DL bias corrected forecasts (See Figures 2 and 3).\n4. In Figure 2, DL-DM and DL models exhibit nearly identical skills in both Anomaly Correlation Coefficient (ACC) and Continuous Ranked Probability Score (CRPS). Could the authors explain why there is no significant difference in performance between DL-DM and DL in this context?\nThank you for your comments. We have now included additional panels within Figure 2 to display variability between different experiments from scenarios DL and DL-DM. Additionally, we have included a brief discussion of these results.\n[L 126-133] “DL-DM and DL models exhibit similar skill across all leads which can be the direct result of the model converging on similar weights during the training process, additionally the models within scenarios (DL and DL-DM) are averaged across many experiments as well as across CONUS which can diminish the variability. Spatial differences between a single DL (EX24_G) and DL-DM (EX29_G) model identify that variability in skill occurs across the entire region, albeit differences may be minimal (Fig 2; C,D). EX29_G mostly shows higher ACC than EX24 except in some of the eastern regions (Fig. 2; D). EX29_G also consistently shows better CRPS than EX24_G over all of CONUS, suggesting higher ensemble forecast skill (Fig. 2; D). ”\nIt appears that the number of runs varies between different models or methods. What was the rationale behind selecting these specific numbers of runs for each model?\nThe number of runs is determined by the plausible inputs that have been used for each model scenario (as described in the Deep learning architecture and experiment setup in the Methods section). For DL model scenario, we include combinations of soil moisture and atmospheric reanalysis inputs at different leads. For DL-DM models scenario, it additionally includes atmospheric and soil moisture reforecast from dynamic models as inputs, which resulted in a larger number of model runs. For scenario DM-BC_DL, we directly bias-correct DM forecasts using soil moisture reanalysis.\nAdditionally, could the differing numbers of runs influence the comparative results and conclusions drawn from this analysis?\nThank you for pointing this out. The performance of all model runs is shown in Figure 3. The primary purpose of creating multiple runs was to understand the effects of different inputs (reanalysis or reforecast) on the forecast performance, allowing us to more robustly compare different model scenarios. Differing numbers of runs does not influence the results and conclusion. Figure 3 shows EX28, as a DL-DM hybrid model, has the best performance compared to all the other model runs, which is consistent with the conclusion we draw from the averaged results in Figure 2, showing the DL-DM scenario provides the best performance.\n5. Figures 4&5 : To provide a clearer understanding of the improvements achieved by the DL-DM approach, it would be beneficial to include additional figures that depict the differences in skill metrics between DL and DL-DM models. For instance, presenting a figure that shows the skill of subplot (C) minus the skill of subplot (E) would effectively highlight the enhancements brought about by the hybrid approach.\nWe agree that the inclusion of a differencing plot would be more helpful to readers in understanding differences between skill. We have subsequently changed Fig. 4 and Fig. 5 to include the difference of skill.\n6. The manuscript examines flash droughts using only a single case study from China. To enhance the robustness and generalizability of the findings, it is recommended that the authors include multiple flash drought events from diverse geographical regions and climatic conditions. Incorporating a broader range of case studies would provide more comprehensive evidence of the method’s effectiveness and applicability across different environmental contexts.\nIn the original manuscript, we included three cases studies in US, China, and Australia (one in main text and two in the Supplemental). To address your comments, we have now moved the two case studies from the supplemental material into the main text (Figs. 7-8).\nMinor comments:\n1. In the introduction (Lines 73-75), the authors mention various methods used in the study. It would enhance clarity if the authors could provide more detailed explanations of the differences between these methods, particularly distinguishing between the hybrid DL-DM approach and DL-processed DMs. This additional information would help readers better understand the unique contributions of each method.\nThanks for pointing this out. In the revised manuscript, we provided clarification on the differences between different models in the Introduction.\n[L 88-96] “In this paper, we explore the extent to which DL, trained on publicly available data, can provide accurate long-lead soil moisture and drought predictions at the subseasonal timescale and address whether the hybrid DL-DM models can outperform the reanalysis-driven DL models, DL or statistical postprocessed (bias-corrected) DMs, and the current state-of-the-art DMs for predicting subseasonal RZSM and drought. To address this question, we designed experimental runs, including hybrid DL-DM models, which take relevant variables from different DMs reforecast and reanalysis data as inputs, DL or statistical postprocessed DMs, which directly bias correct DM forecasts, and reanalysis driven DL models, which take reanalysis as inputs. All models take RZSM reanalysis as the output target. ”\n2. Fig 3. I would recommend adding color legends for the runs.\nWe have now included a color legend within Figure 3.\n3. Line 99: or would it be because of the poor atmospheric prediction skills?\nThe decrease in performance cannot be associated with poor reforecast atmospheric predictions skill because only reanalysis atmospheric predictors were used in this model runs (See Supplemental Table S1 or Fig. 3). More specifically, reforecast atmospheric predictors were only used for only EX27-EX29 model runs for models predicting Weeks 1 and 2 RZSM. Our analysis and discussion of Week 3 experiments (Fig. 3) does not include any atmospheric predictions from GEFSv12. We have now added a few clarifying thoughts about potential for low skill.\n[L 139-146] “Subsequent experiments after EX25 in Fig. 3 included reanalysis atmospheric predictors and lower skill may be due to a combination of a.) diminished effects of previous atmospheric conditions at longer leads (low signal-to-noise ratio), b.) evolving land-surface conditions (precipitation and other hydrometeorological processes) during the three-week forecasts which are more impactful on RZSM variability or c.) inclusion of too many predictors which inhibits model learning through overfitting. Specifically, research has identified that increasing the number of inputs in a neural network can create a too complex model prone to overfitting and is unable to generalize to out-of-training samples38.”\n4. lines 101-105, I'm a little confused about this sentence, could the authors make it clearer?\nWe apologize for the lack of clarity and agree that the wording can be improved upon. We have now included a better description to clarify our intent.\n[L 147-153] “The best performing model runs included only inputs from observation RZSM and/or recursive RZSM predictions. The highest skill model at week 3 (EXP29), contains inputs of only lagged RZSM reanalysis and recursive RZSM predictions. While this set of inputs is highly similar to EX18, the recursive RZSM predictions from weeks 1-2 within EX29 come from models which were trained on only week 1 lag of reanalysis atmospheric conditions and short-lead DM atmospheric conditions (weeks 1-2); therefore, it is likely that this hybrid model EX29 retains the important atmospheric signals within the recursive RZSM predictions in which it is trained (Supplemental Table S1).”\n5. On lines 114-115, the authors note that the Mississippi and Ohio Valley regions exhibit poorer performance. Are there specific factors or characteristics of these regions that contribute to this reduced performance?\nThank you for bringing up this point. We have now provided an additional investigation into potential causes of low prediction skill of RZSM within this region through an autocorrelation analysis and discuss results in the discussion section. Additionally, we have placed it into context with previous research which also has identified this issue.\n[L 290-298] “We observed that DL-DM predictions in some regions have low skill such as in the Mississippi River, Ohio Valley, and Southwest regions within CONUS which has been identified in previous post-processing research for temperature post-processing. Horat and Lerch (2024) implemented a series of experiments with two different UNet architectures and found these same regions have low skill relative to climatology, regardless of the architecture. Here we argue that for RISE-UNet, low RZSM autocorrelation can be a primary reason for both low skill within physically based dynamical forecasts as well as post-processed outputs (Supplemental Fig. S10). In the Mississippi, Ohio Valley, and Southwest regions, anomaly autocorrelation (𝜌 ) k for lags (k) greater than 21 days have noticeably smaller values which indicate that previous anomalies are less likely to persist in the future, thus inhibiting predictability based solely on memory.”\nReference:\nHorat, N. & Lerch, S. Deep Learning for Postprocessing Global Probabilistic Forecasts on Subseasonal Time Scales. Mon Weather Rev 152, 667–687 (2024).\n6. In the Results section, it is unclear what reference data was used for evaluating the forecasts. Although descriptions are provided in the Discussion, it would be beneficial to include a brief mention of the reference data either in the Results section or the Introduction.\nWe agree that presenting information regarding the ground-truth dataset earlier in the manuscript would be beneficial. We have now added information regarding what we are verifying skill against within the results section.\n[L 117-118] “We trained the model separately for each week lead 1-5 and verified skill results against the physically-based GLEAM reanalysis (see Methodology).”\nReviewer #2:\nReviewer #2 (Remarks to the Author):\nThis manuscript examines using deep learning to predict soil moisture and drought conditions beyond the typical two-week forecast horizon. The proposed approach combines historical reanalysis data and dynamic model forecasts to predict weekly root-zone soil moisture anomalies up to five weeks ahead, benchmarked against dynamical forecasting systems and deep learning post-processed predictions.\nMajor Concerns:\n1. Fundamental Issues with Problem Formulation\n- The study does not address the operational reality that current soil moisture status is not immediately available, due to inherent delays in gathering observations and running analysis models for both soil moisture and atmospheric states\nThank you for this identifying potential roadblocks between research and operations. We agree that soil moisture reanalysis is currently not immediately available due to delays in gathering observations and running analysis. Our best performing model takes antecedent soil moisture and dynamic model forecasts. Both ECMWF and GEFSv12 dynamic models are already in an operational setting (through real-time forecasts). The soil moisture data (GLEAM) were not updated in real time. It deserves a separate study to bridge the gap between research and operations by including additional data sources such as soil moisture data from National Aeronautics and Space Administration Short-term Prediction and Transition Center (NASA SPoRT). We included a detailed paragraph within the discussion section.\n[L 329-340] “ECMWF and GEFSv12 subseasonal forecasts are not reliant on assimilated RZSM data during initialization which allows them to work efficiently within an operational setting. This can be a disadvantage due to the importance of accurate land surface initializations which are required for skillful forecasts of temperature and precipitation (Deangelis et al. 2020). As future models are developed and data assimilation procedures implemented into land surface model initializations, we can expect that the quality of the dynamical model forecast will increase and further post-processing using the RISE-UNet can subsequently correct model biases. Alternatively, future research can incorporate hydrological model which include a more extensive set of land surface parameterizations which account for land surface fluxes such as Noah-MP (Niu et al. 2011; Yang et al. 2011). Additional constraints which may impact operational performance include using lagged RZSM from reanalysis which may not be directly available on a weekly basis due to data assimilation and quality control procedures. This can be circumvented by incorporating near real-time soil moisture data such as from National Aeronautics and Space Administration Short-term Prediction and Transition Center”.\nReferences:\nDeangelis, A. M. et al. Prediction skill of the 2012 U.S. great plains flash drought in subseasonal experiment (SubX) models. J Clim 33, 6229–6253 (2020).\nNiu, G. Y. et al. The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with local-scale measurements. Journal of Geophysical Research Atmospheres 116, (2011).\nYang, Z. L. et al. The community Noah land surface model with multiparameterization options (Noah-MP): 2. Evaluation over global river basins. Journal of Geophysical Research: Atmospheres 116, 12110 (2011).\n- The methodology for generating ensemble predictions lacks clarity, particularly regarding how ensemble spread is achieved without apparent stochastic elements, which is crucial for CRPS optimization\nThank you for your comments. We included an explanation in the Methodology section regarding how the ensemble spread was achieved and how CRPS loss function was calculated.\nexp\n[L 512-517] “For clarification, the MAE and standard deviation calculation is computed on each batch and the process begins with partitioning each batch size of 66 into 6 groups of 11-members. On each group, the MAE is calculated on each ensemble predictions and the standard deviation is computed across the 11-member ensemble. This set of functions allows the model to adjust model weights which both decrease the error and increase the spread of the forecast. This set of functions allows the model to adjust model weights which both decrease the error and increase the spread of the forecast, thereby retaining stochasticity similar to that of physically-based dynamical models”\n2. Limited Scientific Advancement\n- The key findings largely confirm known behaviors: soil moisture's relatively long memory compared to atmospheric processes, and the degradation of atmospheric prediction skill beyond two weeks\nThis study provides innovations for subseasonal soil moisture forecasts based on deep learning with the understanding of soil moisture memory and atmospheric prediction skill. No previous study has yet to create models that can skillfully forecast extreme soil moisture anomalies after two weeks (sources are provided within the introduction Lorenz et al. 2021; Su et al. 2023). The existing understanding of atmospheric prediction skill deteriorate after two weeks and soil moisture memory guided our experimental design (Supplemental Table S1). Existing literatures on deep learning for weather and climate predictions focus on developing purely reanalysis driven deep learning model for 14-day weather (e.g., Lam et al. 2023; Bi et al. 2023; Price et al. 2025) or seasonal forecasts (e.g., Ham et al. 2019), while neglected the potential benefits of reforecast from dynamic models and subseasonal forecasts beyond 14-day but within a season. This study addressed these key research gaps by designing and evaluating unique sets of experiments that compare purely reanalysis driven deep learning models (DL), reanalysis and dynamic model reforecast driven deep learning models (DL-DM), deep learning postprocessed/bias corrected dynamic models (DM-BC_DL), statistical postprocessed dynamic models (DM-BC_EMOS), and raw forecasts from state-of-the-art dynamic models (GEFSv12 and ECMWF) for subseasonal soil moisture forecasts. A major finding suggested that the DL-DM hybrid model, which combine antecedent soil moisture and dynamic atmospheric forecast within two weeks, shows the best performance, providing significant insights for developing more skillful environmental predictions beyond solely relying on reanalysis data or postprocessing dynamic models. The new methods and findings from this study demonstrate a major advancement towards skillful subseasonal forecasts beyond two weeks.\nReferences:\nLam, R., Sanchez-Gonzalez, A., Willson, M., Wirnsberger, P., Fortunato, M., Alet, F., Ravuri, S., Ewalds, T., Eaton-Rosen, Z., Hu, W. and Merose, A., 2023. Learning skillful medium-range global weather forecasting. Science, 382(6677), pp.1416-1421.\nBi, K., Xie, L., Zhang, H., Chen, X., Gu, X. and Tian, Q., 2023. Accurate medium-range global weather forecasting with 3D neural networks. Nature, 619(7970), pp.533-538.\nPrice, I., Sanchez-Gonzalez, A., Alet, F., Andersson, T.R., El-Kadi, A., Masters, D., Ewalds, T., Stott, J., Mohamed, S., Battaglia, P. and Lam, R., 2025. Probabilistic weather forecasting with machine learning. Nature, 637(8044), pp.84-90.\nHam, Y.G., Kim, J.H. and Luo, J.J., 2019. Deep learning for multi-year ENSO forecasts. Nature, 573(7775), pp.568-572.\n- The study fails to provide meaningful process-based analysis of how early precipitation forecasts contribute to sustained soil moisture predictability\nThank you for your comments. This study provides new evidence showing early precipitation forecasts contribute to soil moisture predictability (week leads 1-2). It deserves a separate study to perform a process-based analysis of how early precipitation forecasts contribute to sustained soil moisture predictability. To address your comment, we provided a brief description in the introduction section to detail the relationship between precipitation forecasts and soil moisture predictability.\n[L 69-74] “...This low predictability can be attributed to high model errors within precipitation forecasts which drive changes within RZSM through well-known processes of infiltration and recharge, and a poor representation of precipitation impacting the computation of dynamical mass and energy fluxes. Low precipitation errors within forecasts would improve RZSM predictability by appropriately mapping water availability and allow the land surface model to appropriately compute fluxes related to evapotranspiration, runoff, and water storage.”\n- The proposed model lacks mechanistic interpretability that could advance our understanding of soil moisture dynamics\nWe agree that deep learning models often function as black boxes, limiting mechanistic interpretability. The aim of this study is to advance subseasonal forecasts, which shares a common goal with other deep learning focused studies recently published in Nature (e.g., Nearing et al. 2024; Price et al. 2025). This study provides methodological innovations and the findings are useful for accelerating long-range (subseasonal) forecasts beyond 14-day. We used the permutation test to increase the model interpretability. However, given the scope of this work, it deserves a separate study to improve its mechanistic interpretability for advancing our understanding soil moisture dynamics.\nReferences:\nNearing, G., Cohen, D., Dube, V. et al. Global prediction of extreme floods in ungauged watersheds. Nature 627, 559–563 (2024). https://doi.org/10.1038/s41586-024-07145-1\nPrice, I., Sanchez-Gonzalez, A., Alet, F. et al. Probabilistic weather forecasting with machine learning. Nature 637, 84–90 (2025). https://doi.org/10.1038/s41586-024-08252-9\n3. Scientific Merit\nThis work appears to be a straightforward machine learning application study without sufficient novelty or significance for Nature Communications publication standards.\nThis study presents a novel approach for making skillful subseasonal soil moisture forecasts based on the current understanding of soil moisture memory and atmospheric prediction skill. The novelty and significance of this study lie in the following aspects. First, the current approach cannot provide skillful forecast of extreme soil moisture anomalies beyond two weeks. The model presented in this study demonstrates significant advances in improving subseasonal soil moisture forecasts. Second, recent prominent studies on deep learning-based weather, climate, or hydrologic predictions focus on within 10-day (e.g., Bi et al. 2023; Lam et al. 2023; Nearing et al. 2024; Price et al. 2025) or seasonal forecasts (e.g., Ham et al. 2019). This study addresses a critical gap in existing literature concerning deep learning for subseasonal forecasts of soil moisture anomalies, which suggests a major advance in subseasonal forecasts. Third, recent prominent studies on deep learning for weather or hydrological predictions (e.g., Bi et al. 2023; Lam et al. 2023; Price et al. 2025) focus on developing reanalysis data driven deep learning forecast models, while neglected the potential benefits of reforecast from dynamic models in developing deep learning models. This study designed and evaluated novel experiments that compare purely reanalysis driven deep learning models (DL), reanalysis and dynamic model reforecast driven deep learning models (DL-DM), deep learning postprocessed/bias corrected dynamic models (DM-BC_DL), statistical postprocessed dynamic models (DM-BC_EMOS), and raw forecasts from state-of-the-art dynamic models (GEFSv12 and ECMWF) for subseasonal soil moisture anomaly and drought forecasts. A major finding suggested that combining antecedent soil moisture and dynamic atmospheric forecast within two weeks give the best performance, providing useful insights for developing more skillful environmental prediction systems beyond solely relying on reanalysis data or postprocessing dynamic models. In summary, the new methods, findings, and insights provided from this study demonstrate sufficient novelty or significance.\nReferences:\nBi, K., Xie, L., Zhang, H., Chen, X., Gu, X. and Tian, Q., 2023. Accurate medium-range global weather forecasting with 3D neural networks. Nature, 619(7970), pp.533-538.\nLam, R., Sanchez-Gonzalez, A., Willson, M., Wirnsberger, P., Fortunato, M., Alet, F., Ravuri, S., Ewalds, T., Eaton-Rosen, Z., Hu, W. and Merose, A., 2023. Learning skillful medium-range global weather forecasting. Science, 382(6677), pp.1416-1421.\nNearing, G., Cohen, D., Dube, V. et al. Global prediction of extreme floods in ungauged watersheds. Nature 627, 559–563 (2024). https://doi.org/10.1038/s41586-024-07145-1\nPrice, I., Sanchez-Gonzalez, A., Alet, F. et al. Probabilistic weather forecasting with machine learning. Nature 637, 84–90 (2025). https://doi.org/10.1038/s41586-024-08252-9\nHam, Y.G., Kim, J.H. and Luo, J.J., 2019. Deep learning for multi-year ENSO forecasts. Nature, 573(7775), pp.568-572.\nTechnical Issues:\nAbstract and Introduction:\n- Imprecise terminology (e.g., \"Deep learning with historical reanalysis data\" should be \"Deep neural networks that learn from climate reanalysis data\")\nThank you for your suggestion, we have made this alteration within the manuscript within the abstract.\n[L 6-7] “Deep neural networks that learn from climate reanalysis data have produced skillful weather forecasts within ten days.”\n- Inconsistent logic in comparing deep learning + dynamical model performance versus deep learning + reanalysis\nThank you for your comments. However, we are unsure which sections you mentioned in the text are unclear and cannot make any changes.\n- Confusion between organizational names and model references (e.g., ECMWF).\nWe have confirmed from several sources that the subseasonal-to-seasonal model created by ECMWF is conventionally called ECMWF. Ref 1, Ref 2, Ref 3, Ref 4.\nTo alleviate potential confusion, we rename the subseasonal-to-seasonal ECMWF (S2S) model as “ECMWF” within our manuscript in two places.\n[L 155] “...trained on subseasonal-to-seasonal ECMWF reforecast (now referred to as ECMWF) with GLEAM reanalysis.”\n[L 392-393] “...and subseasonal-to-seasonal European Center for Medium Range Weather Forecasting (S2S ECMWF) models reforecast archive. Within this study, we refer to the S2S ECMWF as ECMWF.”\n- Inadequate definition of key terms and metrics\nThank you for your comments. However, we are unable to make any changes due to a lack of feedback regarding which key terms or metrics are unclear.\n- Operational feasibility concerns regarding the use of lagged reanalysis data\nThank you for addressing this topic. We have now included a thorough discussion on the discrepancies between research and operations (with respect to lagged reanalysis data) and we provide potential directions for future research to overcome this issue.\n[L 336-340] “Additional constraints which may impact operational performance include using lagged RZSM from reanalysis which may not be directly available on a weekly basis due to data assimilation and quality control procedures. This can be circumvented in the future by incorporating near real-time soil moisture data such as from National Aeronautics and Space Administration Short-term Prediction and Transition Center.”\n- Mischaracterization of error growth mechanisms in dynamical systems\nDue to our regional analysis across the globe, we discussed common error growth mechanisms in dynamical systems, while a detailed analysis of error growth mechanisms in dynamical systems is out of the scope of this work.\nMethodology:\n- Insufficient description of model architecture and problem setup prior to results presentation\nWe apologize for not sufficiently discussing the model architecture prior to presenting results. We have now included a more thorough description of the model architecture and the order of events during training prior to discussing results.\n[L 106-116] “RISE-Unet accepts an arbitrary number of input predictors (channels) from both reanalysis at weekly lags and reforecast at forecast leads, and passes weighted convolution kernels across each channel to produce intermediate predictions which are concatenated. Additional convolution kernels are passed along all intermediate predictions, and the image size is reduced through pooling operations which allow different representations of the data to be computed within a more restricted field space, allowing the model to identify relevant features at different scales. These intermediate results are concatenated again and the convolution/pooling process can repeats until a minimal image size is reached (in our instance a 48x96 pixel image can beis reduced to a 6x12 pixel image). Once the terminal size is reached, intermediate predictions are upscaled until the original image size is reached. For the RISE-Unet, three separate predictions are computed and the loss function applied over each one independently with the goal of the model to reduce the aggregated loss of all three predictions.”\nWith regards to the problem setup, we believe that we have done a proficient job of detailing the limitations in predictability from all available previous soil moisture research as it relates to the prediction of dry anomalies. We have analyzed and reported the major findings for bias-corrected dynamical reforecasts, raw dynamical reforecasts, and statistical models. To further re-emphasize the problem we have added this new statement\n[L 58-61] “While bias-correction can improve forecast quality for mean conditions, there is still a sufficient need to identify additional methods which improve the predictability of extreme events due to low skill in bias-corrected forecasts after two weeks (Su et al. 2023).”\nReference:\nSu, L., Cao, Q., Shukla, S., Pan, M. & Lettenmaier, D. P. Evaluation of Subseasonal Drought Forecast Skill over the Coastal Western United States. J Hydrometeorol 24, 709–726 (2023).\n- Unclear definitions of key variables and temporal relationships\nWe apologize for any uncertainty of key definitions and temporal relationships. We have clarified the key variables and temporal relationships in the revised manuscript.\n- Imprecise statistical framework for extreme event analysis\nWe believe that our current set of tests and statistical frameworks are sufficient for the analysis of extreme events. Previous studies have incorporated metrics such as the true positive rate/false positive rate (Supplemental Figs. 8-9) for understanding model skill with respect to extreme events, and these metrics are widely accepted in the scientific community (Bathiany et al. 2020, Nearing et al. 2024, McKinnon et al. 2016, Su et al. 2023). We have analyzed our model 3 case studies (Figs. 6-8) and verified that the forecast quality is generally better than dynamical forecasts for mean absolute error, correlation, and structural similarity which provide a diverse set of metrics for forecast quality assessment.\nReferences:\nBathiany, S., J. Hidding, and M. Scheffer, 2020: Edge Detection Reveals Abrupt and Extreme Climate Events. J. Climate, 33, 6399–6421, https://doi.org/10.1175/JCLI-D-19-0449.1.\nMcKinnon, K., Rhines, A., Tingley, M. et al. Long-lead predictions of eastern United States hot days from Pacific sea surface temperatures. Nature Geosci 9, 389–394 (2016).\nhttps://doi.org/10.1038/ngeo2687\nNearing, G., Cohen, D., Dube, V. et al. Global prediction of extreme floods in ungauged watersheds. Nature 627, 559–563 (2024). https://doi.org/10.1038/s41586-024-07145-1\nSu, L., Cao, Q., Shukla, S., Pan, M. & Lettenmaier, D. P. Evaluation of Subseasonal Drought Forecast Skill over the Coastal Western United States. J Hydrometeorol 24, 709–726 (2023).\nThe manuscript requires substantial revision to address these fundamental issues before it could be considered for publication.\nThank you for your feedback. We have made substantial revisions to the original manuscript and have addressed all the issues as specified above.\nReviewer #3:\nReviewer #3 (Remarks to the Author):\nKey results\nThe paper demonstrates substantially improved soil moisture prediction skill on subseasonal time scales (1-5 weeks ahead) over the contiguous US, Australia and China, obtained using novel deep learning (DL) compared to the raw output from state-of-art dynamical prediction models (DMs), as well DMs calibrated with conventional linear statistics. The improved skill is found to be largely attributable to soil moisture memory. Purely empirical DL forecasts are found to perform fairly similarly to a hybrid method that combines DL with DMs, although the hybrid approach obtains the highest skill.\nValidity\nThe technical quality appears high and generally convincing. The main findings are shown to be robust across three different geographical regions. In the main MS evaluation is largely in terms of deterministic skill, using anomaly correlation, and over a short 2-year evaluation period, though other skill metrics are available in the SI.\nOne apparent technical flaw needs to be addressed: The GEFSv12 and ECMWF dynamical model reforecasts (years 2000-2019) use different configurations which means that the same forecast start dates are not the same between both both (L274). The main findings don’t depend on this, but it needs to be addressed because it (rightly or wrongly) indicates an incomplete understanding of subseasonal forecasting that otherwise casts doubt over the study.\nThank you for addressing this uncertainty and we apologize for not initially including a more thorough discussion of our data pre-processing. We had already accounted for this inconsistency between ECMWF and GEFSv12 initializations. We have now provided a description of our methodology for assessing subseasonal forecast with different initializations.\n[L 425-438] “ECMWF and GEFSv12 subseasonal forecasts have different initialization dates which can influence a direct comparison between models. To resolve this issue, we created new ECMWF initialization files which contained the same forecast dates as the nearest GEFSv12 file initialization date. This process involved selecting ECMWF 45-day forecasts which were initialized on Monday and subsetting these forecasts to be consistent with the 35-day forecast dates of GEFSv12 which were initialized on the following Wednesday. For example, ECMWF forecast initialized on Monday December 2, 2019 were subset to match the same 35-day forecast as GEFSv12 which was initialized on December 4, 2019. The process is straightforward and creates a new ECMWF file with a 43-day forecast with the first lead date matching the first lead date from GEFSv12. While the impact on ECMWF skill would be negligible for short-lead forecasts due to the strong influence of initial conditions on forecast skill, longer lead forecasts (e.g., week-4 forecasts) would be more impacted due to two additional days of accumulated error growth. Research has identified that long-lead subseasonal forecasts are highly unskillful regardless of predictor or model; therefore, this methodology to compare models with different initialization dates can be viewed as the best approximator with negligible uncertainties”\nSignificance\nWhile the findings suggest a substantial improvement in soil moisture prediction skill obtained using DL, there is not enough comparison against simple benchmarks (persistence and climatological forecasts), or with the results of previous work, for the significance of the reported DL skill increases to be properly assessed. The very short length of the testing period (only 2 years) is a limitation that needs to be acknowledged.\nWe have now included a separate analysis of CRPSS with the reference forecast being that of the forecast from GEFSv12 or ECMWF (Supplemental Fig. S6). Additionally, we have included in the manuscript the limitation of short length of testing period.\n[L 193-196] “While the models were trained using large samples from 2000-2015, it is worth noting that the sample size was small (n=104) during the testing period (i.e., dry weekly RZSM anomalies below the 20th percentile from 2018-2019). A longer DM reforecast archive can improve robustness of the evaluations.”\n[L 265-270] “We acknowledge that our sample size is small for the testing set years 2018-2019 (n=104), and for seasonal evaluation the sample size reduces to approximately n=31 which can decrease the predictive potential of the model (Hernandez et al. 2006).”\nReference:\nHernandez, P. A., Graham, C. H., Master, L. L. & Albert, D. L. The effect of sample size and species characteristics on performance of different species distribution modeling methods. Ecography 29, 773–785 (2006).\nThe ACC values obtained here should be compared to previous studies in the published literature on seasonal and subseasonal forecasting of soil moisture. How do they compare in general with previous studies on both the seasonal and subseasonal time scale? The current results will not mean much unless they are set in context of previous work.\nWe have now provided a comparison between previous subseasonal soil moisture forecast studies.\n[L 312-321] “Across CONUS, averaged ACC values at lead times >3 weeks are higher than previous research which assessed multiple subseasonal-to-seasonal forecasts (Duan et al. 2025) and ACC skill in week lead 3 in the Southeast is improved compared to Lesinger et al. (2024). RISE-UNet predictions show an increase in TPR for dry anomalies in the South, High Plains, and Western United States (Supplemental Figs. S9) which has been identified previously to be low within GEFSv12 (Lesinger et al. 2024). For lead times <2 weeks, RISE-UNet has much improved deterministic skill across eastern China when compared to research models within the Subseasonal-to-Seasonal database (Zhu et al. 2019). Predictions are also improved for the western and southeastern Australia regions as well as southeast and west China when compared to the Norwegian Climate Prediction Model (Nair et al. 2024), although a direct assessment is limited due to their research computing monthly average RZSM predictions.”\nReferences:\nDuan, Y. et al. Enhancing sub-seasonal soil moisture forecasts through land initialization. NPJ Clim Atmos Sci 8, 1–13 (2025).\nLesinger, K., Tian, D. & Wang, H. Subseasonal Forecast Skill of Evaporative Demand, Soil Moisture, and Flash Drought Onset from Two Dynamic Models over the Contiguous United States. J Hydrometeorol 25, 965–990 (2024).\nNair, A. S., Counillon, F. & Keenlyside, N. Improving subseasonal forecast skill in the Norwegian Climate Prediction Model using soil moisture data assimilation. Geoscientific Model Development Discussions 2024, 1–29 (2024).\nZhu, H. et al. Evaluation of the subseasonal forecast skill of surface soil moisture in the S2S database. Atmospheric and Oceanic Science Letters 12, 467–474 (2019).\nRelatedly, a general weakness is the lack of physical interpretation of the results. In particular, the results are not stratified by calendar season, which is important for interpretation, and is standard practice. The testing period is only 2 years, but I think such a seasonal breakdown ought to be included in the SI.\nThank you for your comments. We have now completed the seasonal assessment for each region (Supplemental Figs. S2, S4, S5) and included the following texts in the revised manuscript.\n[L 162-166] “An assessment of seasonal variability displays that months December-February (DJF) contribute most to the overall skill for week leads 3-4 (~0.80 ACC), followed by months June-August (JJA) with ACC between 0.60-0.65, September-November (SON) with ACC between 0.50-0.60, and lowest average skill is in March-May (MAM) with ACC between 0.47-0.51 (Supplemental Fig. S2).”\n[L 173-179] “For seasonal variability in China (Supplemental Fig. S4), the highest contribution to overall ACC skill for week leads 3-4 is from SON with average skill across leads (~0.62 ACC), proceeded by DJF (~0.54), MAM (~0.53), and JJA (~0.40). Improvements in skill relative to GEFSv12 display that the Tibetan plateau has no improvements in skill in JJA for weeks 3-4, while the eastern China agricultural region has improved skill during MAM and JJA. Seasonal skill improvements in Australia are primarily contributed from season JJA with averaged week lead 3-4 ACC of ~0.50, proceeded by MAM (~0.40), and similar results between DJF and SON (~0.37) (Supplemental Fig. S5).”\nThe main skill assessment is in terms of deterministic metrics which makes the significance of the main findings hard to judge. It is widely accepted that forecasts need to be expressed probabilistically in order for them to be usable in societal decision making.\nIn Figure 2, we presented the CRPS for evaluating the probabilistic forecasts which has been routinely used in weather forecast studies (Price et al. 2025, Weyn et al. 2021, Li et al. 2021). In the revised manuscript, we also included additional probabilistic forecast evaluations to show CRPSS in the US, China, and Australia (Figure S6).\nReferences:\nPrice, I., Sanchez-Gonzalez, A., Alet, F. et al. Probabilistic weather forecasting with machine learning. Nature 637, 84–90 (2025). https://doi.org/10.1038/s41586-024-08252-9\nLi, Y., D. Tian, and H. Medina, 2021: Multimodel Subseasonal Precipitation Forecasts over the Contiguous United States: Skill Assessment and Statistical Postprocessing. J. Hydrometeor., 22, 2581–2600, https://doi.org/10.1175/JHM-D-21-0029.1.\nWeyn, J. A., Durran, D. R., Caruana, R., & Cresswell-Clay, N. (2021). Sub-seasonal forecasting with a large ensemble of deep-learning weather prediction models. Journal of Advances in Modeling Earth Systems, 13, e2021MS002502. https://doi.org/10.1029/2021MS002502\nData and methodology\nOverall technical quality appears high, with the caveats mentioned in my previous remarks.\nHowever, the labelling of the many experiments is highly confusing, and the selection of results presented in the figures often obscures simple interpretations.\nWhy are different lead times (week 3 and 4) shown in Figs 3-6? Since weeks 3+4 is established practice in subseasonal forecasting, it would be simplest to choose that, or at least present all figures for either week 3 or week 4.\nThank you for your suggestion. For clarity, we present all figures for week 3.\nSecond, the key question to answer is how much this new DL approach can improve on established practice with dynamical models, which is well represented by DM-EMOS. The ACC and MAE for each raw DM (ie without any bias correction) and the DM-EMOS should be included in Fig 3. Assuming the DM-EMOS outperforms DM-raw, then these should be the benchmarks for the case studies in Figs 6, S5, S6, not the raw DM.\nWe have also included both raw DM and EMOS in Fig. 3. EMOS is only marginally better than raw GEFSv12. We have now provided a comparison between raw GEFSv12 and EMOS in Fig. 4. We only trained EMOS over CONUS region; therefore, for case study analysis we only retain GEFSv12 for comparison since this data is available for all regions.\nThe captions in Figs. 6, S5, S6 need to be include the dates of there target periods, not just the initialization dates. It should be clarified that GLEAM refers to the validation data.\nThank you for addressing that we have not included the verification dates. We have now included all of the new verification dates along with the initialization date for all case studies.\nAnalytical approach\nWhy was the CRPS used, as opposed to CRPSS (L395)? The latter is a much more common metric in meteorological and climate studies because it is a Skill Score that expresses the forecast accuracy relative to a benchmark, typically the climatological forecast. Unless the authors can justify the choice of CRPS, and provide a reference where CRPS has been used in the meteorological field, then CRPSS should be used instead of CRPS.\nPrevious studies identifying skill of subseasonal forecasts have utilized the CRPS extensively and we believe that this method is useful for diagnosing the efficacy of models for probabilistic predictions (Price et al. 2025, Weyn et al. 2021, Li et al. 2021). To evaluate models’ relative performance, we now have also included an analysis of the CRPSS over each study region where we evaluated over the climatology forecast (Supplemental Fig. S6).\nReferences:\nPrice, I., Sanchez-Gonzalez, A., Alet, F. et al. Probabilistic weather forecasting with machine learning. Nature 637, 84–90 (2025). https://doi.org/10.1038/s41586-024-08252-9\nLi, Y., D. Tian, and H. Medina, 2021: Multimodel Subseasonal Precipitation Forecasts over the Contiguous United States: Skill Assessment and Statistical Postprocessing. J. Hydrometeor., 22, 2581–2600, https://doi.org/10.1175/JHM-D-21-0029.1.\nWeyn, J. A., Durran, D. R., Caruana, R., & Cresswell-Clay, N. (2021). Sub-seasonal forecasting with a large ensemble of deep-learning weather prediction models. Journal of Advances in Modeling Earth Systems, 13, e2021MS002502. https://doi.org/10.1029/2021MS002502\nIt needs to be clarified that the results are not stratified by calendar season (L475). While year-round results may still be informative, there may be a contribution from mean seasonality that can inflate the skill estimates, depending on the details of how the seasonal cycle has been subtracted. This needs to acknowledged and addressed.\nWe have included stratified the results by calendar season (Supplemental Figs. S2, S4, S5). Since the evaluation results are based on anomalies and the seasonal cycles have been subtracted, the year-round results would not inflate the skill estimates. Additionally, we have added more explicit information on how seasonal anomalies were constructed.\n[L 369-372] “Seasons are defined according to conventional categories of DJF (December, January, February), MAM (March, April, May), JJA (June, July, August), SON (September, October, November). Seasonal distributions were created by grouping observations within the same calendar season. The seasonal mean was then subtracted to isolate anomalies.”\n[L 409-411] “For each lead time, seasonal distributions were created by grouping forecasts whose lead dates fell within the same calendar season. The seasonal mean was then subtracted to isolate anomalies which are used as model inputs.”\nMulti-model ensembling of DMs boost skill and is established practice in seasonal and subseasonal forecasting. Why is that not done here?\nWe acknowledge that multi-model ensemble of DMs may boost skill and is an established practice in S2S forecasting. We have now included a separate analysis in Supplemental Fig. S1.\n[L 133-135] “Multi-model ensembling is commonly used to increase forecast skill through leveraging the strengths of each model forecast, but we identified that this method did not considerably increase ACC skill when combining ECMWF and GEFSv12 (Supplemental Fig. S1).”\nSuggested improvements\nL48: Give “beyond week 2” in days or just say “beyond two weeks” and clearly define the “subseasonal” forecast range in days and weeks, here at the outset\nWe have now changed our wording to “beyond two weeks (>14 days)” and clarified the definition of the subseasonal time range.\n[L 48-49] “Previous studies have consistently identified limited subseasonal forecast skill beyond two weeks (>14 days) in...”\nL66 & 67: explain “shorter” hirer in terms of weeks, and “medium range” for the Nature audience\nWe have now clarified our meaning of “shorter” and “medium range”.\n[L 81-84] “Recent studies have also developed DL models for subseasonal predictions of atmospheric variables, while the method did not outperform state-of-the-art DMs particularly at shorter leads (<10 days)Click or tap here to enter text.37. Due to SMM and high predictability of weather at the medium rangeClick or tap here to enter text. 37(3-15 days), DL has a great potential to skillfully predict subseasonal soil moisture and droughts.”\nL84/L88: Clearly define “weekly” A table with daily equivalents of the weeks may be useful. Is week 1 comprised days 1-7 for example. These definitions vary in subseasonal studies.\nThank you for your suggestion. We have now added a new table which describes what each week and corresponding days encompasses (Supplemental Table S6).\nL83: Explain near here that the skill results are averaged over forecasts made at weekly intervals across the entire calendar year, for the two-year period 2018-19. Those are important for there interpretation of the results and should be highlighted here in addition to the methods section.\nThank you for your suggestion. We have now added clarification of the interpretation for skill results within the Results section and the Methods section.\n[L 118-120] “For clarification, skill results are obtained for each forecast week lead made at weekly intervals across the calendar year/season for the two-year period 2018-2019.”\n[L 575-576] “For all evaluation metrics, forecast skill results are averaged over forecasts made at weekly intervals across the entire calendar year for the two-year testing period 2018-2019.”\nL108: Spatial *distribution” of skill\nWe have now added “distribution” to the sentence.\n[L 158] “We evaluate the spatial distribution of skill between DL-DM_G, EMOS, and GEFSv12 (Fig. 4).”\nL189: The statement “Even after DL postprocessing, dynamic models GEFSv12 and ECMWF show no skill for predicting RZSM beyond week 2” is at odds with Figs 2, 4 & 5. An ACC threshold of 0.5 is only a rule of thumb for skillful forecasts.\nThank you for bringing up this point of note. We have now changed this sentence to more objectively compare the data.\n[L 262-265] “Dynamic models GEFSv12 and ECMWF show limited ACC skill for predicting RZSM beyond 2 weeks when looking at regional averages (Fig. 2), but seasonal assessments identified spatiotemporal improvements in skill can be significant (Supplemental Figs. S2,S4,S5).”\nL237: Please describe how the “seasonal mean” is defined, with respect to a particular weekly forecast target period. Are these rolling 3-month averages?\nWe have now included a description of the seasonal mean calculation within the Methods section.\n[L 369-372] “Seasons are defined according to conventional categories of DJF (December, January, February), MAM (March, April, May), JJA (June, July, August), SON (September, October, November). Seasonal distributions were created by grouping observations within the same calendar season. The seasonal mean was then subtracted to isolate anomalies.”\n[L 409-411] “For each lead time, seasonal distributions were created by grouping forecasts whose lead dates fell within the same calendar season. The seasonal mean was then subtracted to isolate anomalies which are used as model inputs.”\nL247: Clarify that Delta_tmp is the diurnal temperature difference. Note that passage of a frontal system could reduce this, as well as increase it.\nThank you for your suggestion. We have now included that Delta_tmp is the diurnal temperature difference and briefly mentioned frontal system effects.\n[L 381-383] “∆tmp considers the diurnal temperature difference and a higher/lower temperature difference over land explicitly implies that there was either less/more cloud coverage that day or a frontal storm system occurred which can either increase or decrease ∆tmp.”\nL287: Please state the lat-lon domains used, and clarify the lead time definition of the forecast weeks, eg using a table.\nThank you for your suggestion. We have now stated the lat-lon domains within Supplemental Table S5. Additionally we have provided more clarification for the lead week and lead day definition of forecast weeks (Supplemental Tables S6). We have also included these references within the Methods section.\nL299: Clarify the meaning of “channels”\nWe have now clarified the meaning of channels within the Methdology.\n[L 468-469] “For clarification, each channel is a single predictor variable for a single lead or lag (e.g., GLEAM RZSM inputs from six lag weeks will be six different channels).”\nL656: This Santoalla reference link is broken.\nWe apologize that the reference page has been subsequently deleted. We have now updated this reference link from ECMWF (ECMWF. Anomaly Correlation Coefficient - Forecast User Guide - ECMWF Confluence.\nhttps://confluence.ecmwf.int/display/FUG/Section+6.2.2+Anomaly+Correlation+Coefficient.).\nACC Figures: The color scales need to be standardized across the figures, using a simple color scale, eg -0.4 to +0.9\nWe have now standardized ACC figures using a simple color scale, -0.4 to +0.9\nClarity and context, References\nThe manuscript is well written, generally clear and set well in context of past work and references previous literature appropriately.\nThe study of Lesinger et al (2024) is highly relevant and ought to be compared against more closely. The authors do cite it.\nWe have now additional comparisons with Lesinger et al. (2024).\n[L 314-316] “RISE-UNet predictions show an increase in TPR for dry anomalies in the South, High Plains, and Western United States (Supplemental Fig. S8) which has been identified previously to be low within GEFSv12 (Lesinger et al. 2024).”\n[L 324-326] “Based on spatial variations in seasonal skill (Lesinger et al. 2024) which are also identified in our research, additional deep learning methods may need to account for different predictors being more impactful at different times of year.”\nReferences:\nDuan, Y. et al. Enhancing sub-seasonal soil moisture forecasts through land initialization. NPJ Clim Atmos Sci 8, 1–13 (2025).\nLesinger, K., Tian, D. & Wang, H. Subseasonal Forecast Skill of Evaporative Demand, Soil Moisture, and Flash Drought Onset from Two Dynamic Models over the Contiguous United States. J Hydrometeorol 25, 965–990 (2024).\nL312: UNet has previously been used by Horat and Lerch (2024, MWR) for S2S forecasting, which ought to be cited and commented on, including the key difference between their model and the RISE-UNet.\nThank you for your suggestion. We have subsequently added a new paragraph in the Discussion which details methods used by Horat and Lerch (2024), results obtained from their research, and a comparison with our work.\n[L 290-309] “DL-DM predictions in some regions have low skill such as in the Mississippi River, Ohio Valley, and Southwest regions within CONUS which has been identified in previous post-processing research for temperature post-processing. Horat and Lerch (2024) implemented a series of experiments with two different UNet architectures and found these same regions have low skill relative to climatology, regardless of the architecture. Here we argue that for RISE-UNet, low RZSM autocorrelation can be a primary reason for both low skill within DMs as well as post-processed DMs (Supplemental Fig. S10). In the Mississippi, Ohio Valley, and Southwest regions, anomaly autocorrelation (ρ_k) for lags (k) greater than 21 days have noticeably smaller values which indicate that previous anomalies are less likely to persist in the future, thus inhibiting predictability based solely on memory. While our research has combined both atmospheric and land surface conditions which have been proposed to increase forecast skill48, we find that combining these features did not adequately increase skill within these regions. A natural extension to this research would be the incorporation of additional predictors which influence precipitation and RZSM variability including sea-surface temperature or high-frequency modes of variability (e.g., North Atlantic Oscillation, Madden Julian Oscillation). The global UNet models utilized by Horat and Lerch (2024) contained traditional skip connections between layers which is different than our implementation of a regional RISE-UNet architecture which concatenates intermediate predictions between layers. It is currently not known what differences in skill for post-processed reforecasts between a traditional UNet vs. RISE-UNet, and we leave this to future research to comprehensively examine ensemble predictions trained with CRPSexp loss to identify if the decreased training cost associated with traditional UNet outweighs potential improvements over RISE-UNet.”\nI did not see a link to the codes used in the MS. Including one would be very beneficial to researchers.\nWe apologize for not having included the link to the codes in the original manuscript. All the codes are shared through OSF, a DOI minting repository, and can be obtained through this link https://osf.io/6y4kh/ (DOI: 10.17605/OSF.IO/6Y4KH)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The author has addressed all of my comments. I have no more questions.", "reviewer_2": "Thanks for the point-by-point respond. Yet, the following questions remain to be clarified.\n1. The methodology for generating ensemble predictions lacks clarity. The uncertainty may come from a) meteorological forcings, which are in-deterministic beyond the weather range, b) errors in representing land processes. However, the authors only apply a perturbation to the neural nets' weights for generating ensembles, which is not justified or well clarified.\n2. The study still fails to provide detailed, process-based analysis of how early meteorological forecasts, current soil moisture status contribute to sustained soil moisture predictability. While these aspects are declared, we may need further analysis showing how deep learning help to enhance our understandings, and possibly lead to transparent, process-based understandings.\nOverall, I still consider the work here to be trivial application of deep learning for a prediction task. The novelty and depth could not meet the criteria of the journal.", "reviewer_3": "I have looked over the authors' replies to my comments, along with the revised version of the manuscript.\nI am not convinced that this apparent technical flaw has been addressed properly:\nMy original comment:\nOne apparent technical flaw needs to be addressed: The GEFSv12 and ECMWF dynamical model reforecasts (years 2000-2019) use different configurations which means that the same forecast start dates are not the same between both both (L274). The main findings don’t depend on this, but it needs to be addressed because it (rightly or wrongly) indicates an incomplete understanding of subseasonal forecasting that otherwise casts doubt over the study.\nWhile agreeing with my comment, the revised MS (L412-438) still confuses real-time forecasts and reforecasts. Until the authors can demonstrate that they understand the distinction, and can provide a clear description of the data analysis procedure used to construct the combined GEFSv12 and ECMWF multimodal reforecast combination, I cannot recommend publication.\nIn particular:\n1. L412-413: \"ECMWF extended range reforecasts are initialized from ERA5-reanalysis and produce a global 11-member ensemble twice weekly (Monday and Thursday) with a lead time of 45 days between 1995-present.\"\nThis is incorrect. The *real-time forecasts* are initialized *2015-present* twice weekly (Monday and Thursday). Each forecast in that range has a set of reforecasts associated with it, initialized on the same *calendar date* as the forecasts, but for the previous 20 years, initialized from ERA. BUT these will not generally be on Mondays and Thursdays in those past years!\n2. L427-430: The example given is impossible: \"For example, ECMWF forecast initialized on Monday December 2, 2019 were\nsubset to match the same 35-day forecast as GEFSv12 which was initialized on December 4, 2019.\"\nThis would imply that the ECMWF real-time *forecast* was used, not the reforecasts. If that were the case, then the range of possible years would only be 2015-present (not 1995-present or 2000-present), and the full 51-member real-time ensemble could have been used, not the 11-member reforecast set which the paper says was used.\n3. L433-438: \"While the impact on ECMWF skill would be\nnegligible for short-lead forecasts due to the strong influence of initial conditions on forecast skill,\nlonger lead forecasts (e.g., week-4 forecasts) would be more impacted due to two additional days of\naccumulated error growth. Research has identified that long-lead subseasonal forecasts are highly\nunskillful regardless of predictor or model; therefore, this methodology to compare models with\ndifferent initialization dates can be viewed as the best approximator with negligible uncertainties.\"\nThis is incorrect. The week-4 forecasts would *not* be impacted much by an additional 2 days of lead time; the accumulated error growth due to two additional days at such long lead (21 days ahead vs 23 days ahead) would be negligible.\nFurther, the statement \"long-lead subseasonal forecasts are highly unskillful regardless of predictor or model\" seems to contradict the findings of the study."}, "author_response": "Response to the Reviewers for the Manuscript Entitled “Skillful Subseasonal Soil Moisture Drought Forecasts with Deep Learning-Dynamic Models” (Ref.: NCOMMS-24-67423A)\nWe thank you and the reviewers for your continued consideration of our manuscript entitled “Skillful Subseasonal Soil Moisture Drought Forecasts with Deep Learning-Dynamic Models\" (Manuscript ID: NCOMMS-24-67423A). We appreciate the opportunity to revise and resubmit, and we are grateful for the constructive feedback that has helped us improve the quality and clarity of our work. We have carefully addressed all comments and concerns raised by the reviewers.\nBelow, we provide a detailed point-by-point response. To aid clarity, original reviewer comments are in black, our responses are in blue, and changes in the manuscript are referenced where applicable.\nReviewer #1\nThe author has addressed all of my comments. I have no more questions.\nResponse: We thank Reviewer #1 for the positive evaluation.\nReviewer #2\nThanks for the point-by-point respond. Yet, the following questions remain to be clarified.\n1. The methodology for generating ensemble predictions lacks clarity. The uncertainty may come from a) meteorological forcings, which are in-deterministic beyond the weather range, b) errors in representing land processes. However, the authors only apply a perturbation to the neural nets' weights for generating ensembles, which is not justified or well clarified.\nResponse: We appreciate this important point and have significantly clarified the ensemble generation methodology in the revised manuscript (see Section Deep learning architecture and experiment setup, Lines 553–567).\n“Weather and climate ensemble forecast systems employ various perturbation methods to account for uncertainty in initial conditions and allow for the quantification of uncertainties due to errors in initial conditions, model formulation, and physical parameterizations (Wang et al. 2014, Zhu et al. 2005, Stenstrud et al. 2000). More specifically, when using a singular realization contrived from assimilated observations and processed through a physically-based land and atmospheric model, uncertainty can be imparted from instrument calibration error or model structural errors associated with parameterizations of complex processes. Small errors within initial conditions have been observed to rapidly amplify in weather forecasting and restrict the intrinsic predictability limit of forecast systems (Lorenz 1969);Click or tap here to enter text. therefore, perturbation methods allow for inclusion of a wider range of plausible initial conditions to evaluate uncertainty in ensemble forecasts through sensitivity experiments and probabilistic evaluations. The RISE-UNet DL architecture accepts an ensemble of stochastically generated initial conditions and passes these inputs through a series of operators (convolution kernels, squeeze-and-excitation, and max-pooling), non-linear activation functions (ReLU), and dropout layers to produce stochastic ensemble forecasts for each lead week. The CRP loss function is computed over each set of 11 Sexp ensemble members generated from the RISE-UNet which allows for gradient-based model weight adjustments based on the aggregated loss.”\nMeteorological Forcing Uncertainty: Our ensemble prediction framework already incorporates meteorological uncertainty by using ensemble forecasts from the dynamic model (GEFSv12 or ECMWF) as inputs. Each ensemble forecast include multiple members (11 for GEFSv12 and ECMWF), each with different initial conditions and perturbations. These variations reflect the inherent uncertainty in medium- and sub-seasonal atmospheric forecasts. We now explicitly describe this in the Methods section and cite relevant literature supporting this ensemble approach to meteorological input uncertainty.\nModel Structural (Epistemic) Uncertainty: To represent uncertainty arising from the deep learning model structure itself, we use Monte Carlo dropout during inference. This involves enabling dropout layers at test time and running multiple stochastic forward passes through the model. Each pass samples from a different subnetwork, thereby producing a distribution over outputs. This technique is a widely used and computationally efficient approximation of Bayesian inference (Gal & Ghahramani, 2016). It does not perturb the weights directly but introduces stochasticity in the activation pathways.\nIn summary, our ensemble generation reflects both input forcing variability (through meteorological ensemble members) and model uncertainty (through stochastic dropout at inference), and we have clarified these mechanisms and their contributions more rigorously in the revised manuscript.\nReferences:\nGal, Y. and Ghahramani, Z., 2016, June. Dropout as a Bayesian approximation: Representing model uncertainty in deep learning. In international conference on machine learning (pp. 1050-1059). PMLR.\nLorenz, E. N. The predictability of a flow which possesses many scales of motion. Tellus A: Dynamic Meteorology and Oceanography 21, 289 (1969).\nStensrud, D. J., Bao, J. W. & Warner, T. T. Using Initial Condition and Model Physics Perturbations in Short-Range Ensemble Simulations of Mesoscale Convective Systems. Mon Weather Rev 128, 2077–2107 (2000).\nWang, Y. et al. A New Method for Generating Initial Condition Perturbations in a Regional Ensemble Prediction System: Blending. Mon Weather Rev 142, 2043–2059 (2014).\nZhu, Y. Ensemble forecast: A new approach to uncertainty and predictability. Adv Atmos Sci 22, 781–788 (2005).\n2. The study still fails to provide detailed, process-based analysis of how early meteorological forecasts, current soil moisture status contribute to sustained soil moisture predictability. While these aspects are declared, we may need further analysis showing how deep learning help to enhance our understandings, and possibly lead to transparent, process-based understandings.\nResponse: We appreciate this important point. In the response, we have strengthened the manuscript by adding a new section (Section Results, Lines 147–160) focused on process-based interpretation of the drivers of predictability. Specifically, we now include:\nSensitivity Attribution: Besides performing permutation-based feature importance analysis, we included SHAP (SHapley Additive exPlanations) to quantify the relative influence of early dynamic model forecasts and antecedent soil moisture status on the DL-DM forecast performance. The analysis shows that initial soil moisture and early-week temperature anomalies (e.g., week-1) have a larger influence on predictions at week 1, aligning with physical expectations from prior hydrological memory studies (See Supplemental Fig. S2).\n“SHAP values were computed for week lead 1 inputs to further estimate each input feature’s contribution to the final model prediction (see Methods). Using GEFSv12 as the forecast, it was observed that the most influential features impacting model outputs included GLEAM RZSM at week lag 1 as the most influential feature, followed by ERA5 t at week lag 1, GEFSv12 max reforecast at week lead 1, GLEAM RZSM at week lag 2, and geopotential height (z ) from both 200 ERA5 and GEFSv12 (Supplemental Fig. S2). The higher density of SHAP values for GLEAM RZSM week lag 1 which corresponds with negative SHAP values indicates that wetter soil conditions will lead to a drying of soils in predicted RZSM. While this may seem counterintuitive, it is probable that the model has learned that saturated conditions lead to future drying. GLEAM RZSM at week lag 1 also identifies that wetter soils can increase the predicted RZSM output as seen with a slightly lower density of positive SHAP values. ERA5 t at week lag 1 identifies that max lower temperature values lead to an increase in saturation in predicted RZSM. Z results from 200 both ERA5 and GEFSv12 identify that below normal heights, associated with troughs, lead to an increase in RZSM predicted values (towards wetness) as observed with a slightly higher density of values with positive SHAP.”\nWe also included an additional figure (Supplemental Fig. S5) to show prediction skill as a function of initial soil moisture percentiles (wet/dry conditions). This helps to reveal which states of the land surface lead to higher predictability, and how long that influence persists. The results suggest that wet initial conditions can produce higher week-4 predictability than dry ones, consistent with physical expectations about moisture retention, although the skill is regionally dependent as identified in other studies (Koster et al. 2021, Koster et al. 2010).\nWe agree with the reviewer that transparent, process-based understanding is critical. While deep learning models are often criticized as black boxes, these analyses show that our model captures hydrologically meaningful mechanisms. We believe these additions help to link data-driven insights with well-known land-atmosphere dynamics, and we hope this addresses the reviewer’s concern regarding physical interpretability.\nReferences:\nKoster, R. D., A. M. DeAngelis, S. D. Schubert, and A. M. Molod, 2021: Asymmetry in Subseasonal Surface Air Temperature Forecast Error with Respect to Soil Moisture Initialization. J. Hydrometeor., 22, 2505–2519, https://doi.org/10.1175/JHM-D-21-0022.1\nKoster, R. D., et al. (2010), Contribution of land surface initialization to subseasonal forecast skill: First results from a multi-model experiment, Geophys. Res. Lett., 37, L02402, doi:10.1029/2009GL041677.\nOverall, I still consider the work here to be trivial application of deep learning for a prediction task. The novelty and depth could not meet the criteria of the journal.\nResponse: We respectfully disagree with this assessment and appreciate the opportunity to clarify the novelty and broader significance of our work. Our study goes beyond a routine application of deep learning by demonstrating that a hybrid framework combining dynamic model outputs and deep learning (RISE-UNet) can produce skillful root-zone soil moisture forecasts up to five weeks in advance, which is a known challenge for both dynamic and statistical models.\nKey contributions include:\n A novel hybrid forecast architecture that integrates dynamic model outputs (GEFS, ECMWF) with reanalysis-driven soil moisture memory via deep learning.\n Demonstration of consistent and substantial skill gains over both raw and postprocessed ECMWF/GEFS forecasts across major flash drought events in the U.S., China, and Australia.\n We introduce a DL-DM ensemble forecast scheme that captures forecast uncertainty using ensemble forecasts from dynamic models to represent meteorological forcing uncertainty and Monte Carlo Dropout to quantify DL uncertainty. This approach is effective on capturing forecast uncertainties (as evaluated by the probabilistic forecast metrics).\n Attribution of forecast skill at weeks 3 to 5 to both early dynamic model input and antecedent land surface conditions, offering a mechanistic understanding of S2S predictability.\n Use of SHAP analysis to provide interpretable insights into the roles of meteorological and land-surface predictors at different lead times, bridging statistical skill with physical reasoning.\nGiven the difficulty of achieving extended-range hydroclimate predictability, we believe this work offers a meaningful scientific contribution with clear novelty, especially in operational and process-informed drought forecasting.\nReviewer #3\nI have looked over the authors' replies to my comments, along with the revised version of the manuscript. I am not convinced that this apparent technical flaw has been addressed properly:\nMy original comment:\nOne apparent technical flaw needs to be addressed: The GEFSv12 and ECMWF dynamical model reforecasts (years 2000-2019) use different configurations which means that the same forecast start dates are not the same between both (L274). The main findings don’t depend on this, but it needs to be addressed because it (rightly or wrongly) indicates an incomplete understanding of subseasonal forecasting that otherwise casts doubt over the study. While agreeing with my comment, the revised MS (L412-438) still confuses real-time forecasts and reforecasts. Until the authors can demonstrate that they understand the distinction, and can provide a clear description of the data analysis procedure used to construct the combined GEFSv12 and ECMWF multimodal reforecast combination, I cannot recommend publication.\nResponse: We thank the reviewer for pointing this out and apologize for the confusion in the previous revision. We agree that the distinction between real-time forecasts and reforecasts, especially in the context of GEFSv12 and ECMWF, must be clearly articulated to ensure scientific rigor and clarity.\nClarification:\nGEFSv12 reforecasts are initialized once per day, providing consistent temporal coverage over 2000–2019. ECMWF reforecasts are initialized twice weekly (Monday and Thursday), each with a 20-year hindcast using ERA. In our original example, we mistakenly referenced a real-time ECMWF forecast instead of its corresponding reforecast (e.g., December 2, 2019). This has now been corrected.\nRevised Method Description (L458–L465):\n“To align ECMWF and GEFSv12 reforecasts, we matched lead times based on valid forecast dates, rather than attempting to align initialization dates. For example, an ECMWF reforecast initialized on Monday, Dec 2, 2019 (and its 20-year hindcast ensemble) would be compared against a GEFSv12 reforecast initialized on Dec 2, 2019, even if ECMWF did not produce a forecast on every intervening day. We verified that this approach preserves forecast lead equivalency (e.g., week-3 lead = days 15–21), and we tested sensitivity to potential day-of-week effects, which showed negligible differences in skill scores across the majority of grid cells across CONUS.”\nThis correction removes the confusion between real-time and reforecast data and ensures that the ECMWF and GEFS comparison is methodologically valid. We believe the revised explanation now accurately reflects subseasonal forecasting protocols.\nIn particular:\n1. L412-413: \"ECMWF extended range reforecasts are initialized from ERA5-reanalysis and produce a global 11-member ensemble twice weekly (Monday and Thursday) with a lead time of 45 days between 1995-present.\" This is incorrect. The *real-time forecasts* are initialized *2015-present* twice weekly (Monday and Thursday). Each forecast in that range has a set of reforecasts associated with it, initialized on the same *calendar date* as the forecasts, but for the previous 20 years, initialized from ERA. BUT these will not generally be on Mondays and Thursdays in those past years!\nResponse: We appreciate the reviewer’s comments. The confusion regarding the initialization of real-time forecasts versus reforecasts has been clarified in the revised manuscript. We have updated the manuscript to accurately reflect the distinction between real-time and reforecast initialization, as well as to correct the phrasing:\nRevised Statement (L439-445):\n“The operational ECMWF real-time extended range forecasts are a set of global 51-member ensemble produced twice weekly (Monday and Thursday) with a lead time of 45 days for years 2015-present. For our experiment, we utilized the ECMWF reforecasts which use the same model version and configuration of the ECMWF real-time forecast but are initialized with past known conditions from ERA5 and are available at 1.5° resolution with a smaller set of 11-ensemble members67. These 11-member ensembles are the produced for the previous 20 years for the same calendar date as the real-time forecast and are hereafter referred to as ECMWF reforecast.”\nWe have ensured that the revised explanation properly reflects the methodology behind ECMWF’s real-time forecasts and reforecasts. The revised manuscript now distinguishes between real-time forecasts and reforecasts to avoid any further confusion.\n2. L427-430: The example given is impossible: \"For example, ECMWF forecast initialized on Monday December 2, 2019 were subset to match the same 35-day forecast as GEFSv12 which was initialized on December 4, 2019.\" This would imply that the ECMWF real-time *forecast* was used, not the reforecasts. If that were the case, then the range of possible years would only be 2015-present (not 1995-present or 2000-present), and the full 51-member real-time ensemble could have been used, not the 11-member reforecast set which the paper says was used.\nResponse: We have now clarified our description of steps to ensure that there is no confusion about subsetting procedures. In our previous clarification, we addressed that we utilized ECMWF reforecasts which are not the ECMWF real-time forecasts. We have now changed the wording to explicitly indicate that we are using the reforecasts.\nRevised Description (L458–L465):\n“This process involved selecting ECMWF 45-day reforecasts which were not initialized on the same day as GEFSv12 and subsetting these reforecasts to be consistent with the 35-day forecast dates of GEFSv12. For example, ECMWF reforecast initialized on Monday December 2, 2019 were subset to match the same 35-day forecast as GEFSv12 which was initialized on December 4, 2019.”\n3. L433-438: \"While the impact on ECMWF skill would be negligible for short-lead forecasts due to the strong influence of initial conditions on forecast skill, longer lead forecasts (e.g., week-4 forecasts) would be more impacted due to two additional days of accumulated error growth. Research has identified that long-lead subseasonal forecasts are highly unskillful regardless of predictor or model; therefore, this methodology to compare models with different initialization dates can be viewed as the best approximator with negligible uncertainties.\" This is incorrect. The week-4 forecasts would *not* be impacted much by an additional 2 days of lead time; the accumulated error growth due to two additional days at such long lead (21 days ahead vs 23 days ahead) would be negligible.\nResponse: We appreciate this insightful correction. You are right that for subseasonal forecasts at week-4 lead times, the marginal impact of a two-day difference in initialization is minimal due to the saturation of forecast error growth at such extended ranges.\nRevised Text (L466-L470):\n\"While short-lead forecasts are sensitive to initialization timing due to the strong influence of initial conditions, at longer leads (e.g., week-4), the marginal degradation in skill from a two-day difference in initialization date is negligible. Therefore, aligning forecast periods across GEFSv12 and ECMWF, even with minor mismatches in start dates, remains a reasonable approximation for comparative skill analysis.\"\nFurther, the statement \"long-lead subseasonal forecasts are highly unskillful regardless of predictor or model\" seems to contradict the findings of the study.\nResponse: Thank you for noting this inconsistency. We agree that the prior statement was poorly worded and does not reflect the findings of our study. Our intention was to refer to the baseline skill of standalone dynamical models at subseasonal lead times, which is generally low, motivating the need for hybrid and deep learning-based approaches.\nWe have revised the statement to better align with our results (L470-L474):\n\"While baseline subseasonal forecast skill from dynamical models alone tends to decline significantly beyond two weeks, our results show that hybrid approaches, combining early-lead dynamical forecasts and antecedent soil moisture states within deep learning models, can substantially improve predictability at weeks 3–5.\"\nThis revision removes the overgeneralization and emphasizes the contribution of the hybrid method to extending forecast skill at long leads, which aligns with our findings."}, {"round": 3, "reviewer_comments": {"reviewer_2": "Most of my concerns have been addressed.", "reviewer_3": "I have looked over the authors' responses to my comments on the previous version. They have addressed them all to my satisfaction."}, "author_response": ""}]

Decision letter