Assimilating multi-satellite snow data in ungauged Eurasia improves the simulation accuracy of Asian monsoon seasonal anomalies

Properly initializing land snow conditions with multi-satellite data assimilation (DA) may help tackle the long-standing challenge of Asian monsoon seasonal forecasts. However, to what extent can snow DA help resolve the problem remains largely unexplored. Here we establish, for the first time, that improved springtime snow initializations assimilating the Moderate Spectral Imaging Satellite (MODIS) snow cover fraction and the Gravity Recovery and Climate Experiment (GRACE) terrestrial water storage data can improve the simulation accuracy of Asian monsoon seasonal anomalies. Focusing on the western Tibetan Plateau (TP) and mid- to high-latitude Eurasia (EA), two regions where multi-satellite snow DA is critical, we found that DA influences the monsoon circulation at different months depending on the regional snow–atmosphere coupling strengths. For the pre-monsoon season, accurate initialization of the TP snow is key, and assimilating MODIS data slightly outperforms jointly assimilating MODIS and GRACE data. For the peak-monsoon season, accurate initialization of the EA snow is more important due to its long memory, and assimilating GRACE data brings the most pronounced gains. Among all the Asian monsoon subregions, the most robust improvement is seen over central north India, a likely result of the region’s strong sensitivity to thermal forcing. While this study highlights complementary snow observations as promising new sources of the monsoon predictability, it also clarifies complexities in translating DA to useful monsoon forecast skill, which may help bridge the gap between land DA and dynamical climate forecasting studies.


Introduction
The Asian monsoon affects more than 60% of the world's population with profound socio-economic implications [1], yet its skillful seasonal forecast, with 6 Current address: Department of Civil and Environmental Engineering, Princeton University, Princeton, NJ, United States of America 7 Current address: University of California, Los Angeles, CA, United States of America 8 Author to whom any correspondence should be addressed. a lead time of one to several months, remains very challenging [2][3][4]. Current state-of-the-art dynamical climate forecasts put emphasis on the memory of tropical oceans and sea surface temperatures (SSTs) data assimilation (DA) for the monsoon predictability [5][6][7][8][9]. However, properly intializing land conditions to fully utilize the weeks' to months' of land memory [10] for the monsoon forecast is largely unexplored. Despite increased satellites monitoring of global land snow, soil moisture, and vegetation conditions [11][12][13], land is still the least-constrained earth system component by observations in dynamical climate models [14,15]. This has been limiting our understanding of the intrinsic limit of seasonal monsoon predictability [10], as well as the additional skill potentially attainable from land DA [16,17].
The link between Eurasian snow and the Asian monsoon has long been recognized [18]. Positive (negative) snow anomalies tend to be followed by negative (positive) monsoon anomalies [19], highlighting the need to properly initialize Eurasian snow for skillful monsoon forecasts. However, previous snow DA studies have mainly focused on assimilating satellite data to improve point-to basin-scale snow estimates for hydrological predictions [20][21][22], while global multi-satellite snow DA tailored to addressing dynamical climate forecasting problems have been relatively rare [23,24]. Station-based snow depth and satellite snow cover have been assimilated in some climate reanalyses [13,25], but they cannot be separated from these reanalyses to establish their contribution to climate forecasts. The lack of data has resulted in a lack of research efforts tailored to understanding monsoon forecast skill from the snow DA perspectives. Recently, Zhang et al [23] and Zhang and Yang [26] (hereafter Z1416) assimilated the Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover fraction (SCF) and the disaggregated Gravity Recovery and Climate Experiment (GRACE) terrestrial water storage (TWS) data into a state-of-the-art land surface model (LSM). The resulting two global snow DA products potentially contain improved snow cover information from optical observations and improved snow mass information from gravity measurements, providing valuable data sources to tackle the unresolved challenge of dynamical monsoon forecast. Lin et al [16] used these products to initialize the land components in a suite of ensemble-based climate modeling experiments, and found the temperature prediction can be improved for the Tibetan Plateau (TP) and the highlatitude Eurasia (EA) regions. However, to the best of our knowledge, no existing studies have utilized similar satellite-constrained snow data for monsoon forecast studies. This paper pioneers the use of multi-satellite snow DA data for seasonal forecast of precipitation in the Asian monsoon region. Our analyses put emphasis on snow DA over the TP and the EA, two regions featuring sparse rain gauges and pronounced snow DA updates, to demonstrate the mechanistic pathway for snow DA to influence the monsoon forecast. We also scrutinize the impacts of assimilating different satellite data, to understand how to maximize skill improvement from different sorts of snowrelated satellite observations that may be assimilated in the future [24,27]. Additionally, as the methodology development and refinement of global snow DA studies is an area of active research [22][23][24], we also offer a new perspective to assess the accuracy of DA products by independently analyzing the precipitation forecast skill, whereas conventional DA assessments commonly use highly uncertain snow observational data as references and thus have been proven problematical [24,26]. Our analyses focus on springtime snow initialization (March 1st), because the role of snow in seasonal climate prediction is much less understood during the cold-to-warm season transition [16,28]. This transitional time coincides with the so-called 'spring barrier' [29] when SSTs are limited in their predictability, thereby rendering it critical to seek for new sources of predictability. Although there is no guarantee that new sources of predictability exist, it is worth examining the issue since improved prediction would be quite important for people who live in the monsoon region.
This study shows, for the first time, that the simulation accuracy of Asian monsoon seasonal anomalies can be improved by assimilating complementary multi-satellite data on snow cover fraction and snow mass. The results suggest that multi-satellite snow DA may be an underutilized source of predictability, which can offer detectable improvements to the Asian monsoon seasonal forecast. Moreover, satellitedependent, region-dependent, lead time-dependent, and the DA uncertainty-related skill improvement are also carefully examined through our experimental designs. This study may help bridge the gap between land DA and seasonal climate prediction studies, with potentially wide socio-economic implications.

Models and experimental design
We design five suites of ensemble-based climate modeling experiments (  [16]. Details of the DA implementations can be found in Z1416 [23,26].

Land initialization
Role of multi-satellite DA: OL (Exp1) 'Open-loop' experiment: land initialization without assimilating satellite data MOD (Exp2) "MODIS" DA experiment: land initialization with global assimilation of the MODIS SCF data [23] GRAMOD (Exp3) "GRAMOD and MODIS" DA experiment: land initialization with joint assimilation of the GRACE TWS and MODIS SCF data [26] Role of regional DA (figure S2 (stacks.iop.org/ERL/15/064033/mmedia)): Land initialization with GRAMOD applied to the TP only (figure 1); Other regions' land initialization is the same as OL EA_only (Exp5) Land initialization with GRAMOD applied to the EA only (figure 1); Other regions' land initialization is the same as OL initial condition. Z1416 used the Ensemble Adjustment Filter (EAKF) for DA; the prior distribution for the model states (e.g. snow) is estimated using 40member ensemble forcing to feed the model [23]. More DA details can be found in Z1416 [23,26]. We caution the readers that these experiments are different from actual dynamical seasonal forecasts that initialize all components for fully coupled runs; instead, the experiments are designed to provide an isolated investigation of the land impacts only. Koster et al [31] designed similar modeling experiments to assess land impacts, and our experiments are different from theirs in that the SSTs are also prescribed with reanalysis data while Koster et al [31] used climatology SSTs.
To consider atmospheric chaos in the coupled forecast, the atmospheric component is initialized with time-shifted ERA-Interim data at eight time points (0000 UTC within ±4 d of 1 March), where the data are pre-processed to CAM5-consistent grids. This ensemble generation technique is commonly used in operational seasonal climate forecasting [32], producing eight members for ensemble-based experiments. The same land initialization is consistently used across the eight members. More details of the model configurations can be found in Lin et al [16].

Regional DA experiments
In addition to the first three experiments to examine the benefits of assimilating different satellite observations (MODIS and GRACE), we also design two other suites of experiments to understand the relative importance of regional snow DA (Exp 4-5, table 1). To apply DA updates to a specific region, the columns, plant functional types (PFTs), and land units in the DA land restart files are first modified where variables outside the predefined geographic bounding box are replaced with those from the OL restart file. To ensure the replaced variables are physically consistent with the OL snow conditions without imposing model imbalance errors, we also replace snow-related variables for water balance calculation, and soil/vegetation/temperature variables for energy balance calculation with OL. The geographic boxes for regional DA are defined in figure 1; the initial snow depth difference between regional DA and OL is shown in figure S2 to illustrate the effectiveness of regional DA. In total, 280 six-month model runs (i.e. 5 suites × 8-member ensembles × 7 years). The model time step is daily, and the total simulation data size is~15 TB.

Statistical analyses and reference datasets
Due to the limited available snow DA products (7 years from 2003 to 2009) and possible noise from the chaotic atmospheric background, we jointly use the r 2 and RMSE metrics, bootstrapping statistical tests, and five reference precipitation datasets to ensure the robustness of the detected signal. Equations to calculate the skill metrics are shown below for one reference precipitation; the same equations are used for the other four additional reference datasets.
For each season (i.e. three month) into the forecast period, r 2 i for each ensemble member is first computed at each grid cell using 21 data samples (3 months × 7 years) obtained from modeled and observed precipitation, denoted as P m and P o , respectively (equation (1); δP m and δP o denote the standard deviation). r 2 is calculated by averaging r 2 i for the 8 ensemble members (equation (2); n refers to the ensemble size), and the skill difference between DA and OL is computed using equation (3). Negative correlation is set to zero before the square operation in equation (1) to reduce sampling noise [33], although we found the influence of this procedure is minimal to our results.
To test the null hypothesis that r 2 diff is not significantly different from zero, at each grid cell, a 1000member bootstrapping is performed to randomly resample the 21 values in the P m vector in equation (1). This leads to 1000 r 2 diff,shuffle values with equations (2)-(3). The null hypothesis is rejected if the percentage of r 2 diff,shuffle (α) is exceeded by r 2 diff with a twotailed distribution. Since we used 21 data samples (3 months × 7 years), the skill in this analysis incorporates both interannual and month-to-month variability; greater r 2 suggests better explained total variance. Thus, the interpretation for the skill should differ from those using anomaly correlation skill, which only deals with interannual variance (see discussions in Results and Text S3).
RMSE is additionally computed with equation (4) (N = 21 for sample size; n = 8 for ensemble size). Compared to r 2 , RMSE is more obscured with model and observational biases [31]. However, DelSole and Shukla [34] showed that model biases and skills are interrelated, which necessitates reduced RMSE as another target. Therefore, RMSE is jointly used in this study. To reduce the noise existing in the raw RMSE field, five reference datasets are used to compute RMSE with equations (4)- (5), and only the consensus percentage RMSE (PRMSE) is presented in figure S5-we define a consensus only when the signs of PRMSE diff are agreed among ≥4 references. Once a consensus is reached, PRMSE diff is averaged among the agreed observations as the final results. Regions with annual climatological precipitation ≤200 mm are masked (figure S5).
The five reference monthly precipitation datasets include data from the Global Precipitation Climatology Project (GPCP), the Global Precipitation Climatology Center (GPCC), the Tropical Rainfall Measuring Mission (TRMM), the University of Delaware gridded precipitation (UDEL), and the National Oceanic and Atmospheric Administration (NOAA) Precipitation Reconstruction over Land (PRECL). All datasets are interpolated to the model resolution before assessment.  figure 1(a)), and only confined areas of the Rocky Mountains and the Alaska show increased SWE (blue colors). This DA behavior is in line with alleviating the positive snow bias at the high-latitude [35] and the cold bias problem over the TP [36] as previously reported. Global area-weighted averages (bar plots, figure 1(a)) suggest DA is more impactful in certain years than others, such as 2004-2005 and 2008-2009. It is noteworthy that assimilating the two satellite observations leads to different spatial patterns-GRAMOD brings prominent snow depth (SD) updates to the high latitudes ( figure 1(b), right) while MOD does not. This is because GRACE uniquely measures the Earth's gravity change, and so it contains snow mass information that is absent in optical sensor observations, e.g. the MODIS snow cover fraction (SCF), especially when snow cover reaches~100%. Over the TP, however, MOD more effectively reduces SCF than GRAMOD ( figure 1(b), blue), which suggests inconsistent MODIS and GRACE DA updates as previously reported [16]. We note figure 1 and figure S1 only show the DA and OL difference, but it does not aim to address whether the snow estimates by DA are necessarily closer to the truth or not, or which satellite DA is better. This is partly to avoid repeating previous assessments that have reported MODIS DA issues over boreal forests [23] and GRACE DA issues over certain basins where TWS disaggregation is complex [26]. Moreover, due to the unknown 'truth' , conventional DA assessments have been ambiguous when different observational references were used [24,26]. Therefore, we will independently assess the DA uncertainty and effectiveness through the coupled modeling results of this study in Section 3.4.

Multi-satellite snow
Globally, the most prominent DA updates occur in the western TP, snow-covered and snow-free transitional areas (50 • N-55 • N), and the higher latitudes (figure 1). These geographic regions overlap with regions where models have high uncertainty in simulating snow, due to their relatively high seasonal snow variations and sparse rain gauge observations to serve as accurate meteorological inputs. DA mathematically weighs the relative uncertainty between model simulations and satellite observations [23], and hence greater DA updates are introduced there. The next sections focus on the TP and the EA (green boxes in figure 1) for further analyses due to their relatively well-known mechanistic link with the Asian monsoon [37]. Moreover, the TP and the EA, without constraints from satellite DA, were previously shown to have limited land-derived climate predictability [31]. As DA provides critical constraints to the land snow states there, we hypothesize that additional gains may be expected to complement the best-attainable skill reported in previous studies.

Contrasting thermodynamic effects introduced by snow DA in Tibet and EA
To facilitate physical understanding, DA-introduced SWE change is viewed as external forcing brought into the system, which persists into the forecast months. DA reduces snow in TP and EA (negative figure 2(a)), leading to more absorbed shortwave radiation (reduced albedo) and less latent cooling (reduced warm-season snow melt). Hence, the net thermal effect of DA over these regions is equivalent to a positive sensible heating that warms up the local atmosphere ( figure 2(b)). Interestingly, this sensible heating starts to warm up the atmosphere in different months over the two regions. During March and April, the dominant DA forcing is seen over the TP ( figure 2(b), blue lines in MOD and GRAMOD), whereas in May and June, the dominant DA forcing is shifted to the EA ( figure 2(b), red lines, only prominent in GRAMOD). This timing difference is explained by the regional snow-atmosphere coupling differences [16,38]-the TP receives strong incident solar radiation to immediately translate changes in snow to atmospheric responses. By contrast, strong snow-atmosphere coupling in the EA only becomes established when solar radiation increases in warmer months to drive snow melts (although this EA forcing is only prominent and in GRAMOD still summer; see figure S3 for the forced air temperature responses to July). Evidently, although the springtime snow DA is applied globally, due to the snow-atmosphere coupling that differs from region to region, the period over which DA has a significant thermal effect to the atmosphere also differ, so as to the assimilation of optical versus gravity satellite data.
We next analyze the monsoon circulation changes because the thermal forcing by snow DA is expected to influence the land-sea thermal contrast (LSTC), which is the fundamental driver of the monsoon circulation [41] (additional assessments on the monsoon onset timing are provided in table S1 and text S1). To simplify the characterization of the highly complex Asian monsoon dynamics, two widely used monsoon circulation indices are computed (more in Text S2). One is the WY index [39] (defined as the zonal wind shear between U850 and U200) and the other is the monsoon Hadley circulation index, MH index [40] (defined as the meridional wind shear between V850 and V200). Among the three experiments, GRAMOD generally shows the least biased circulation indices (figures 2(c)-(d)) with the least uncertainty in May and June (figures 2(c)-(d)); this overlaps with the time when the EA forcing dominates ( figure 2(b)) and the circulation shifts direction for monsoon onset. Compared to GRAMOD, the contribution of MOD is more consistently seen in the meridional (MH index) than in the zonal direction (WY index). This may be due to that the MOD signal is geographically confined to the TP, different from the widespread, zonally aligned EA forcing. Yet due to the high elevation of TP, the additional 4.8 W m −2 equivalent sensible heating (5%-6% of the total TP springtime sensible heating) can still exert strong controls on the monsoon meridional circulation [1,42], although this signal quickly decreases after May ( figure 2(b)). In June and July, GRAMOD becomes more effective in reducing the monsoon circulation bias than May and June, which may be ascribed to the latent heat release aloft the Indian subcontinent ( figure S4). This condensationrelated mid-troposphere warming was known to play important roles in monsoon strengthening in a positive feedback [43], which potentially explains the pronounced impact of GRAMOD on the monsoon circulation.

Robust skill improvement for monsoon forecast for central north India
We then quantify the monsoon forecast skill obtained from snow DA. To ensure the results robustness, we jointly use r 2 and root-mean-square-error (RMSE) metrics, five precipitation reference datasets, and bootstrapping techniques (Methods) for the analysis. Since the results are in general consensus in terms of spatial pattern (see figure S5), we only present r 2 with one reference data for simplicity in figure 3. Across four forecast seasons, significantly improved forecasts (dotted, p < 0.05) can be seen over the Southeast Asia (Myanmar and Thailand), southern China (near the Mei-Yu front), and the Indian subcontinent. The onset of the Asian monsoon typically occurs over the Indochina peninsula in late April, over southern China in mid-May, and over the Indian summer monsoon (ISM) in May to June [44]; the improvement for each of these regions covers a time span after the regional monsoon onset. The only exception is for the JJA ISM forecast in MOD, which is a satellite-dependent behavior that will be discussed in Section 3.4.
Among all the Asian monsoon sub-regions, central north India (CNI) features the most robust improvement independent of the metrics and precipitation references used (green box in figure S5). CNI is the region that receives the most rainfall in the ISM [46]. Hence, an improved area-averaged precipitation forecast by~8% ( figure S6(a)) could have profound agricultural and socio-economic implications. CNI is part of the South Asian Summer Monsoon, which was shown to exhibit the strongest sensitivity to LSTC changes among all global monsoons [45]. Therefore, it may be subject to less complex interplays among boundary forcings [43,47] compared to the East Asian Summer Monsoon, potentially explaining its most robust improvements. For CNI, DA can reduce the RMSE by up to 30% and increase r 2 by up to 15% (Taylor diagram in figure S6), however, it is important to note that the increased variance consists of both interannual and month-to-month variability components (Methods). Further decomposition (figure S7) shows that DA increases the explained interannual variance only for the short-lead but not for the long-lead forecasts, while bias is consistently reduced at all lead times (figure S6). Although interannual variance is arguably the most critical target for skill improvement and our limited data samples prevent a clear separation of this, DelSole and Shukla [34] showed that model bias is intrinsically related to its forecast skill. To this end, our results are still promising for the CNI forecast as a pilot demonstration with snow DA, but future studies should use more DA samples to improve this skill quantification (see text S3 for a discussion). Some locations with degraded skill are also noted, such as the edge of the Asian monsoon boundaries (blue colors, figure 3(a)). This may be partly due to the known issues with Z1416 [23,26], and partly be related to the possible atmospheric/sampling noise. Since a perfect improvement is by no means expected, we define δDA as percentage area with significant improvement minus that with degradation [28], to more objectively summarize the DA contribution. A positive (negative) δDA suggests DA improves (degrades) a forecast, while zero suggests no net contribution from DA. For the Asian monsoon with an area extent of 11.29 × 10 6 km 2 ( figure 3(a)), δDA is +20% to +30%, with slightly higher values in GRAMOD ( figure 3(b)). For the ISM (areal extent: 4.3 × 10 6 km 2 ), the strongest signal (δDA exceeds +50% in GRAMOD) is seen during the monsoon transitional period, overlapping with when the EA forcing signal dominates (figure 2). δDA is generally small for non-monsoonal areas; highlights for some other global monsoon regions are summarized in figure S8.

Relative importance of satellite data (MODIS and GRACE) in the ISM forecast
To understand the DA uncertainty and its impact on monsoon forecast, here we focus on the precipitation skill over the Indian subcontinent, where the most robust signal is located. In addition to MOD and GRAMOD (figure 4(a)), we isolate the exclusive contribution of assimilating GRACE data in figure 4(b), which shows GRACE DA brings significant improvement to the two-to three-month lead forecast (MJJ and JJA). However, GRACE DA barely offers improvement, or even slightly offsets the contribution from MODIS, at short-lead time in spring (blue in figure 4(b), MAM and AMJ). Since the springtime improvement is mostly related to the forcing over the TP ( figure 2(b)), the slightly degraded skill indicates that GRAMOD may have provided inferior snow estimates over the TP, compared to MOD. In fact, when conducting multi-satellite snow DA, Zhang and Yang [26] applied no rules of preference if MODIS and GRACE data show opposite signs in updating snow (see figure 1(b) for the inconsistency problem over the TP). However, due to the unknowability of 'truth' [26], it was inconclusive as to which DA method produced better snow estimates. Here, by examining the precipitation forecast, and similarly the temperature forecast by Lin et al [16], we are able to infer that over the TP, only assimilating MODIS data produced better snow estimates than jointly assimilating GRACE and MODIS data. We further hypothesize that the inferiority of the latter approach may be related to the complex GRACE TWS composition over the TP, including glacier melt, lake changes, and crustal movement [48][49][50], while a simplified TWS disaggregation only considering snow, soil moisture, and groundwater variations was adopted by Z1416. Recent studies suggested that additional considerations for surface water variations in TWS disaggregation can yield better results for GRACE DA in some river basins [51,52]. Our results suggest that similar region-specific modifications, or rules to reject unreliable observations [27], may be needed to improve the GRACE DA method over the TP, in order to achieve more skill for the pre-monsoon forecast.
Similarly, we can infer which DA method produced better snow estimates over the EA by focusing on the long-lead forecast in MJJ and JJA. As GRAMOD significantly outperforms MOD, the critical role of assimilating GRACE data in producing better snow estimates over the EA, and in improving the peak-monsoon ISM forecast is highlighted. This contribution of GRACE DA to the ISM forecast is outstanding, suggesting an overlooked source of the monsoon predictability, upon which better disaggregation of the GRACE TWS signal [51,52] should further improve the forecast. Differing from conventional DA assessments, we use independent precipitation skill analysis to infer which satellite data/methods are more effective. This should have implications in guiding the multi-satellite DA methodology refinements in the absence of 'truth' snow conditions in the future.

Relative importance of regional snow DA in the ISM forecast
We also further assess the relative importance of regional snow DA in the ISM forecast by conducting two more modeling experiments (see Methods & figure S2 for model configuration details). These regional DA experiments (figure 5) show that the TP snow DA is more important for the CNI rainfall forecast in spring (pre-monsoon season), while the EA snow DA is more important for the longlead ISM forecast (the peak monsoon season). This is consistent with figures 2(b) and 4, which is ascribed to the different regional snow-atmosphere coupling strengths, as well as the intrinsic memory difference of the TP and the EA snow. The elevated TP is a region of strong sensible heating [1,42] receiving the primary research focus for the Asian monsoon [37,53,54]. Our results suggest that due to the relatively shorter memory of the TP snow (figure 2(b): forced atmospheric response disappears by May; see similar discussions on the TP snow by Orsolini et al [55]), the contribution of snow DA over the TP may be limited to pre-monsoon season. In comparison, properly initializing the EA snow conditions, especially with the snow mass estimates from the GRACE DA, can offer the most gains to the peak-monsoon ISM forecasts. Despite the relatively simplified GRACE TWS disaggregation method used [26], assimilating GRACE data brings clear skill improvements for the long-lead ISM forecast. Therefore, our results here may provide a lower limit for the best attainable monsoon forecast skill, while future studies may improve the EA snow estimates through using improved GRACE DA methodology [51,52] and/or incorporating other sensor observations for snow DA [24,27].
Interestingly, beyond the Asian monsoon region, we also find some significant improvements in the North African Monsoon (NAF) region (see global monsoon delineations in figure S8), for which δDA can reach +35% in MJJ and JJA. Some recent studies reported the TP surface heating can influence its upstream NAF climate [56], while some viewed the Afro-Asian monsoon as a planetary-scale system that can respond to the same boundary forcings [57]. These suggestions may help explain the detectable improvements in the NAF forecast. However, given our snow DA was applied to a large Eurasian landmass beyond the elevated TP ( figure 1(b)), future idealized modeling experiments tailored to the TP and the NAF are still warranted to better understand this complex interplay.

Conclusion and discussion
Overall, this study offers promising evidence as to the contribution of snow DA to the simulation accuracy of Asian monsoon seasonal anomalies, highlighting an underutilized source of the monsoon predictability. More importantly, we clarify the complexities associated with translating snow DA to useful monsoon forecast skill, including but not limited to the regional snow-atmosphere coupling, the utility/limitation of different satellite observations, and the DA uncertainty, which may be helpful to bridge the gap between the land DA and the dynamical seasonal forecasting studies. Despite such complexities, the snow DA contribution is still very clear, and thus the results presented may be arguably viewed as a lower limit for the attainable skill, upon which future studies could improve.
One immediate future step would be to increase the number of data samples (i.e. using longer time period's snow DA products to serve as the land initializations), to better quantify snow DA's contribution to the interannual variability of the monsoon precipitation. To the best of our knowledge, our study is the pilot one and very few studies have been conducted along this line. To this end, future work is needed to better address this question which may have great implications for the huge population in the monsoon region.