Spatiotemporal distribution of seasonal snow water equivalent in High-Mountain Asia from an 18-year Landsat-MODIS era snow reanalysis dataset

Seasonal snowpack is an essential component in the hydrological cycle and plays a significant role in supplying water resources to downstream users. Yet the snow water equivalent (SWE) in seasonal snowpacks, and its space-time variation, remains highly uncertain, especially over mountainous areas with complex terrain and sparse observations, such as in High-Mountain Asia (HMA). In this work, we assessed the spatiotemporal distribution of seasonal SWE, obtained from a 10 new 18-year HMA Snow Reanalysis (HMASR) dataset, as part of the recent NASA High-Mountain Asia Team (HiMAT) effort. A Bayesian snow reanalysis scheme previously developed to assimilate satellite derived fractional snow-covered area (fSCA) products from Landsat and MODIS platforms has been applied to develop the HMASR dataset (at a spatial resolution of 16 arc-second (~500 m) and daily temporal resolution) over the joint Landsat-MODIS period covering Water Years (WYs) 2000-2017. 15 Based on the results, the HMA-wide total SWE volume is found to be around 163 km on average and ranges from 114 km (WY2001) to 227 km (WY2005) when assessed over 18 WYs. The most abundant snowpacks are found in the northwestern basins (e.g. Indus, Syr Darya and Amu Darya) that are mainly affected by the westerlies, accounting for around 66% of total seasonal SWE volume. Seasonal snowpack in HMA is depicted by snow accumulating through October to March and April, typically peaking around April and depleting in July-October, with variations across basins and WYs. When examining the 20 elevational distribution over the HMA domain, seasonal SWE volume peaks at mid-elevations (around 3500 m), with over 50% of the volume stored above 3500 m. Above-average amounts of precipitation causes significant overall increase in SWE volumes across all elevations; while an increase in air temperature (~ 1.5 K) from cooler to normal conditions leads to an redistribution in snow storage from lower elevations to mid elevations. This work brings new insight into understanding the climatology and variability of seasonal snowpack over HMA, with the 25 regional snow reanalysis constrained by remote sensing data, providing a new reference dataset for future studies of seasonal snow and how it contributes to the water cycle and climate over the HMA region.

Even though both seasonal snow and glaciers are crucial to hydrology and water availability, seasonal snow has arguably received less attention than glaciers in the HMA region. Many studies have addressed the status and changes in glaciers over 30 HMA (e.g. Bolch et al., 2012;Kääb et al., 2012;Sorg et al., 2012;Yao et al., 2012;Lutz et al., 2014;Bolch et al., 2019;Rounce et al., 2020;Shean et al, 2020). For seasonal snow, previous studies have examined the snow extent (e.g. Dahe et al., 2006;Pu et al., 2007;Immerzeel et al., 2009;Tahir et al., 2011;Basang et al., 2017;Wang et al., 2017;Notarnicola 2020), or snow mass and snow depth (e.g. Dahe et al., 2006;Che et al., 2008;Terzago et al., 2014;Dai et al., 2017;Stigter et al., 2017;Bookhagen, 2018. 2020;Ahmad et al., 2019;Kirkham et al., 2019;Xue et al., 2019). In the current literature 35 involving seasonal snow, most of the studies have focused on snow covered area (or extent, which is readily available from satellite-borne remote sensing) instead of snow mass, or have been applied at relatively localized scales or coarse scales. The seasonal snow water storage and its spatiotemporal distribution at finer scales are highly uncertain in HMA, primarily due to the lack of in situ observations and fine-scale snow water equivalent (SWE) datasets over this large domain (Takala et al., 2011;Kirkham et al., 2019). In fact, accurately estimating SWE (at scales of a few kilometers) remains a great challenge 40 worldwide, and it is even more difficult in mountainous regions due to the terrain complexity (Lettenmaier et al., 2015;Dozier et al., 2016;Bormann et al., 2018).
In situ measurements are usually expensive and difficult to install and maintain in HMA and are mostly located in low-lying valleys, thus resulting in a sparse and potentially nonrepresentative network (Winiger et al., 2005;Palazzi et al., 2013;Dozier et al., 2016;Kirkham et al., 2019). In recent decades, satellite observations can provide large-scale estimates of some 45 snowpack properties. However, most of these measured properties, such as snow-covered area (SCA) based on visible and near infrared bands (e.g. Dozier, 1989, Hall et al., 2002, Painter et al., 2009, are only indirectly related to snow mass. While SWE and snow depth can be directly estimated from passive microwave sensors (using retrieval algorithms based on the brightness temperature; e.g. Chang et al., 1987), these estimates are at coarse spatial resolution (e.g. 25 km), and are generally negatively biased in deep snowpacks (Takala et al., 2011;Dozier et al, 2016). Recent applications of C-band 50 synthetic aperture radar (SAR) techniques show promise for snow depth retrieval (Lievens et al., 2019) but are available only over recent years and do not directly provide SWE.
Global atmospheric reanalysis products provide another approach to large-scale SWE estimates as by-products of their land surface schemes. Examples include the Global Land Data Assimilation System (GLDAS, Rodell et al., 2004), Modern-Era Retrospective analysis for Research and Applications (MERRA, Rienecker et al., 2011;MERRA-2, Gelaro et al., 2017), 55 European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis products (ERA-Interim, Dee et al., 2011;ERA5, Hersbach et al., 2020), High Asia Refined analysis (HAR, Maussion et al., 2011;, Japanese 55-year Reanalysis (JRA-55; Kobayashi et al., 2015), and others. SWE estimates in these datasets are found to be generally consistent in their interannual and seasonal variations, but can differ significantly in their magnitudes when evaluated over different regions (Mudryk et al., 2015;Wrzesien et al., 2019), where the uncertainties come from different land surface models and 60 meteorological inputs (Mudryk et al., 2015;Mortimer et al., 2020;Kim et al., 2021). In addition, most reanalysis datasets are not specifically designed for SWE estimation, and usually do not assimilate snow observations. Bian et al. (2019) found https://doi.org/10.5194/tc-2021-139 Preprint. Discussion started: 17 May 2021 c Author(s) 2021. CC BY 4.0 License. many reanalysis datasets overestimate SWE compared to ground observations in the Tibetan Plateau, although part of the differences may come from inconsistent spatial resolution and elevations between in situ and gridded datasets. The performance of these large-scale reanalysis datasets over the full HMA domain has not been fully assessed due to the sparse 65 and uneven in situ station network.
Recent works have contributed to the development of SWE (or snow depth) estimates covering the HMA region based on passive microwave (e.g. Talaka et al., 2011;Dai et al., 2017;Pulliainen et al., 2020) or active microwave measurements (Lievens et al., 2019). Machine learning approaches have been explored to predict SWE using passive microwave data (e.g. Ahmad et al., 2019). Moreover, some studies have explicitly merged snow observations with 70 modeling via data assimilation (DA) to provide more realistic SWE estimates, which is an effective approach in reducing SWE uncertainties especially over the mountains (Largeron et al., 2020). For example, the JRA-55 product assimilates both ground snow depth and satellite retrieved snow cover observations (Kobayashi et al., 2015), and its snow depth estimate is found to have good performance when evaluated over the Tibetan Plateau (Orsolini et al., 2019). The ERA5 dataset also employs in situ snow depth and satellite observed snow cover (limited to elevation below 1500 m) in its snow assimilation 75 scheme. Xue et al. (2019) showed that assimilating SCA products can improve SWE estimates in the HMA region. Passive microwave retrieved SWE along with ground snow depth observations are also assimilated in the GlobSnow (Talaka et al., 2011;Pulliainen et al., 2020) products to provide SWE and snow extent estimates, with mountain areas of high terrain complexity masked out. These are promising approaches to improve the accuracy in SWE estimates over HMA, yet currently there is still a need for large scale SWE datasets at higher resolutions, over longer periods and covering 80 mountainous areas in this region.
To better understand the spatiotemporal pattern and variability in seasonal snowpack over HMA, the so-called High-Mountain Asia Snow Reanalysis (HMASR; Liu et al., 2021) dataset, available for the joint Landsat-MODIS era between Water Years (WYs) 2000 to 2017 (which will be extended to present in later versions), was developed as part of the NASA High-Mountain Asia Team (HiMAT) activities. HiMAT is a multi-investigator effort in developing new datasets to 85 understand cryosphere variability over HMA (Osmanoglu et al., 2017). The HMASR dataset provides daily estimates of SWE, fractional snow-covered area (fSCA) and other snow variables, at a 16 arc-second (~500 m) resolution. SWE estimates are derived by assimilating fSCA from Landsat and MODIS platforms using a previously developed snow reanalysis framework , where the method has been shown in previous applications to provide realistic 1) How much seasonal snow is stored across HMA?
2) How is this snow storage distributed spatially across the major watersheds of HMA?
3) What is the seasonal and interannual variability in amount of snow storage over HMA? 4) How is the amount of snow distributed across elevation? 100

Data and method
This section describes the data and methods used in this study. Section 2.1 introduces the study domain, including the major river basins and mountain ranges in the region. Section 2.2 and 2.3 provide a brief description of the reanalysis method, input data and models used in the development of the HMASR. Finally, a non-seasonal snow and ice mask applied to mask out semi-permanent snow and ice for the assessment of seasonal snow is explained in Sect. 2.4. 105

HMA domain
The HMA domain used in this work is bounded from 27° N to 45° N, and from 60° E to 105° E ( Fig. 1), covering the highest mountain ranges and plateaus (Tien Shan, Pamir, Hindu Kush, Karakoram, Himalayas, and Tibetan Plateau), as well as the headwaters of the main river basins (Syr Darya, Amu Darya, Indus, Ganges-Brahmaputra, Yangtze, and Yellow). Winter westerlies and the summer monsoon are the major moisture sources in this region, significantly influencing the 110 spatiotemporal patterns in snowfall and glacier mass balance. More specifically, the northern and western HMA is dominated by westerlies and receives abundant winter snowfall, while the southern and eastern HMA is dominated by the Indian monsoon from June to September and receives a considerable amount of summer snowfall; the eastern edges of HMA are affected by the East Asia monsoon but with limited impact (Bookhagen and Burbank 2010;Yao et al., 2012;Bolch et al., 2019). Note that in HMASR, outputs are provided for each regular 1° by 1° latitude-longitude tile (within which a regular 115 computational grid of 16 arc-seconds is used), and only tiles likely to contain snow (with a tile-averaged elevation above 1500 m) were selected and processed in the dataset (Fig. 1).

Snow reanalysis scheme 130
A previously developed snow reanalysis methodology  is employed in deriving the HMASR. For brevity, only the key details are repeated here. Prior model and measurement estimates are obtained via the coupled Simplified Simple Biosphere model, version 3 snow model (SSiB3; Sun and Xue, 2001;Xue et al., 2003) and the Liston (2004) snow depletion curve (SDC). Specifically, an ensemble approach is used whereby the model generates prior estimates of snow states (i.e. SWE, snow depth, etc.) and predicted measurements (fSCA) and their uncertainties. Meteorological forcing inputs are bias-corrected, downscaled to the modeling grid (16 arc-second) and perturbed with uncertainty in the ensemble approach, using the methods described in Durand et al. (2008) and Girotto et al. (2014).
To constrain the prior snow estimates on the remotely sensed fSCA observations, a Bayesian update is performed using the Particle Batch Smoother (PBS; Margulis et al., 2015) approach, resulting in posterior snow estimates that are more consistent with the batch of observed fSCA in a given water year. More details on the method are described in Margulis et al. (2015) 140 and Margulis et al. (2019). The lack of in situ SWE data over HMA prevents a thorough verification of the HMASR.
However, previous applications of the snow reanalysis method in similarly complex terrain in the Sierra Nevada of the Western U.S. and the South American central Andes thoroughly compared reanalysis estimates vs. in situ and airbornederived SWE data. Performance in both domains were positive relative to in situ data with values of mean error, root-meansquared error and correlation coefficients of: ~ 3 cm, 13 cm and 0.95 for the Sierra Nevada (Margulis et al., 2016) and ~ 1 145 cm, 29 cm and 0.73 for the Andes (Cortés and Margulis, 2017), respectively. In Margulis et al. (2019), comparison with the Airborne Snow Observatory (ASO) data yielded similar results (mean error, root-mean-squared error, and correlation coefficients of ~ 5 cm, 23 cm, and 0.84).

Meteorological, topographic and land cover data 150
In HMASR, prior surface meteorological inputs were obtained from MERRA-2 at its raw resolution (0.5° by 0.625° latitude-longitude), including precipitation, air temperature, solar radiation, specific humidity, surface pressure and wind speed. The uncertainty models and their parameters used to perform bias-correction and uncertainty perturbation are specified in Margulis et al. (2019) for the HMA region, except that prior ensemble precipitation is perturbed by a lognormal distribution with mean of 1.54 and coefficient of variation (CV) of 0.83 based on the results from Liu and Margulis (2019). 155 Digital elevation model (DEM) data were obtained from the Shuttle Radar Topography Mission (SRTM, http://www2.jpl.nasa.gov/srtm/) 1 arc-second product, and aggregated to 16 arc-second (~500 m) resolution. Gaps in DEM data were filled by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM, version 2) 1 arc-second product (https://asterweb.jpl.nasa.gov/). Land cover data were obtained from the AVHRR global land cover classification dataset (Hansen et al., 2000). Forest cover information was obtained from 160 the Tree Canopy Cover (TCC) product containing the Landsat Vegetation Continuous Fields (https://lcluc.umd.edu/metadata/global-30m-landsat-tree-canopy-version-4; Sexton et al., 2013).
The HMASR dataset used fSCA from both satellites to increase the number of observations with good quality, because of 170 the frequent cloud contamination in some areas of HMA . Cloud screening was performed for both Landsat and MODIS images (as illustrated in Margulis et al., 2019), which were then aggregated to the modeling resolution (16 arc-second). For MODSCAG, erroneous fSCA observations are removed when the viewing angles are significantly offnadir, and a CDF-matching step is applied to reconcile the inconsistency with Landsat images. More details about the fSCA data used in HMASR are described in Margulis et al. (2019). 175 It is worthwhile to note that the fSCA data availability is significantly affected by cloud contamination in some areas of HMA region, especially during the monsoon seasons (June-September) where fSCA measurements are limited over regions such as the Himalayas (Fig. 2). The lack of abundant fSCA data can be a potential limitation in assimilating fSCA observations for these monsoon-affected regions, and therefore leads to higher uncertainty and potential errors in posterior SWE estimates (i.e. where in the limit of no available observations, the posterior will, by construct, equal the prior estimate). 180

Non-seasonal snow and ice mask
A significant fraction of HMA is covered by glacier and semi-permanent snow owing to its extremely high elevation. Thus, 185 it is important to distinguish seasonal vs. non-seasonal snow over land or glacier surfaces. In particular, the reanalysis method used in the development of the HMASR is best suited for seasonal snow characterization, because it relies on the Herein a combination method was used to exclude the non-seasonal snow and ice pixels in HMASR, based on: 1) a glacier mask derived from GLIMS to identify glacierized pixels and 2) a persistent snow mask derived from the HMASR dataset itself. To be more specific on the second mask, pixels with a significant amount of persistent snow were identified, by comparing the annual minimum SWE at a particular pixel to its annual maximum SWE in each year. If the minimum SWE 200 exceeds 10% of the maximum SWE for more than once out of the 18 years, this pixel is considered to be a persistent snow pixel to be masked out in the computation of seasonal snow estimates. The derived glacier and persistent snow masks are combined into a non-seasonal snow and ice mask, which is applied when presenting the spatiotemporal patterns of seasonal SWE and SWE volumes in the following section.

Results and discussion 205
The HMASR dataset is designed to provide a reliable and consistent SWE product that can be used for assessing the spatiotemporal distribution of seasonal SWE over the recent remote sensing record. To present an overall assessment of seasonal snowpack variability in the HMA region using the HMASR dataset, the results are organized as follows: 1) the spatial distribution of seasonal snowpack climatology, at annual and seasonal scales; 2) the temporal distribution of seasonal snowpack volume at basin and domain-wide scales; 3) the elevational distribution of seasonal snowpack storage at the basin 210 scale.

Spatial distribution of peak seasonal SWE climatology
The spatial distribution of SWE is valuable in assessing the regional water storage and how it varies in time. Given the strong seasonal signature of snowpack processes over much of the domain, the pixel-wise peak SWE is a useful metric to quantify the distribution of the maximum amount of snow water mass held in the seasonal snowpack in a given water year. 215 The climatology (18-year average) of pixel-wise peak SWE over the HMA region is depicted in Fig. 3, where Fig. 3a presents only the results for seasonal snow pixels (where non-seasonal snow and ice pixels have been masked out), and snow and ice mask pixels, corresponding to glaciers or permanent snow). The non-seasonal SWE values (Fig. 3b) are expected to be unreliable because the initial conditions for SWE at those locations at the beginning of the dataset are 220 unknown and the lack of full melt-out makes the relationship between fSCA depletion and peak SWE much less direct. In general, seasonal snow is most abundant in the NW region that is directly exposed to westerlies (Fig. 3a). Among the 225 northwestern mountain ranges, the highest climatological peak SWE values are found in Pamir, Karakoram and the western Himalayas, with more than 1 m of peak SWE estimated. A significant amount of peak SWE is also estimated in Tien Shan 3a). The non-seasonal snow and ice is most notable in Karakoram but also evident in a few locations over the Pamir, Tien Shan and western Himalayas (Fig. 3b). 230 In contrast, seasonal snowpack is less abundant in the SE HMA (Fig. 3a), in part because it only receives much of its precipitation in summer from the Indian and East-Asia monsoons, while the winter westerlies have minimum impact.
Shallow snowpack exists over the Hengduan Shan and Tanggula Shan, with low values of SWE estimated (less than 0.2 m).
For the Himalayas and Nyainqentanglha mountain ranges (Fig. 1), which exhibit extremely high elevation and receive significant summer precipitation from the monsoons, high values of SWE are estimated in some locations (Fig. 3b). 235 However, those locations are masked out in the reanalysis through the non-seasonal snow and ice mask (Fig. 3a), because the fSCA observations are persistently high throughout the year (no observed melt-out), show irregular temporal patterns without a clear accumulation-depletion cycle (non-seasonal), or are obscured by clouds between June-Sept. (insufficient measurements), any of which can contribute to unreliable estimates of SWE.
The least abundant seasonal snowpack is estimated in the NE (Fig. 3a), where SWE is only notable over a few mountain 240 ranges such as the Qilian Shan, Kunlun Shan and Eastern Tibetan mountains. Despite their high elevations, most of the NE areas are snow-free or only have shallow and intermittent snow, as a result of being further away from the primary atmospheric moisture sources.
Previous studies have also examined the spatiotemporal distribution in seasonal snowpack, regarding SCA (e.g. Pu et al., 2007;Basang et al., 2017), snow depth and SWE (e.g. Terzago et al., 2014;Bian et al., 2019;Orsolini et al., 2019), and the 245 overall finding is that most existing datasets present consistent spatial pattern at large scales (e.g. regional) but differ greatly in the magnitudes of SWE and snow depth, which implies large uncertainties in snow mass estimates over this data scarce region. Similarly, our dataset exhibits coherent spatial patterns compared to these previous efforts, yet the magnitudes of SWE still show significant uncertainty. A more comprehensive analysis of HMA SWE between multiple products will be addressed in an upcoming intercomparison paper using HMASR. 250

Spatial distribution of peak SWE occurrence dates
The timing of peak seasonal SWE occurrence is associated with climatological (e.g. precipitation) and topographic factors (e.g. elevation), and therefore shows significant heterogeneity over HMA. Figure 4 depicts the pixel-wise peak SWE day of water year (DOWY) climatology map. Highly intermittent snow pixels were excluded, as well as permanent snow and ice pixels via the non-seasonal snow and ice mask. Peak SWE generally occurs between DOWY 100 and DOWY 250 for 255 seasonal snow. The median date of peak SWE is found on March 18th (DOWY 169), where the peak SWE occurs later than February 10th (DOWY 133) for over 90 percent of the domain, with peak SWE occurring later than May 5th (DOWY 217) for only 10 percent of the domain (Fig. 4). The peak SWE DOWY shows a bimodal distribution (Fig. 4, inset) with the earlier peak centered on DOWY 145 and the later peak centered on DOWY 192.
For those mountain ranges in the NW, the northern and western mountain slopes of Tien Shan, the western foothills of occurrences between February 10th and March 18th (Fig. 4). In contrast, the southern mountain slopes of Tien Shan, the majority of the Pamir, Karakoram, and western Himalaya, show a relatively late peak SWE occurrence between March 18th and May 5th (Fig. 4).
For those mountain ranges in the SE and NE, the peak SWE occurrence dates are more diverse (Fig. 4). In the SE, the central 265 and eastern Himalayas, Nyainqentanglha, and Hengduan mountains generally have later peak SWE occurrences (between March 18th and May 5th), except in the southern foothills, where peak SWE tends to occur earlier (between February 10th and March 18th). In the NE, the Eastern Tibetan mountains show the earliest peak SWE occurrence dates (before February 10th), while the Qilian Shan and Kunlun Shan show the latest peak SWE occurrences (after May 5th). Throughout the year, mountains in NW hold the maximum amount of SWE compared to other regions. In SON, the entire HMA region exhibits minimal SWE magnitudes (0.1 m or below) and most regions are snow free (Fig. 5). During this period, SWE starts accumulating in the Tien Shan, Pamir and western Himalayas that are directly facing the westerlies. SWE is also evident in Nyainqentanglha and Hengduan Shan that are associated with the summer monsoons. In DJF, both the overall magnitude and extent of SWE grow significantly, with mean SWE values up to 0.5 m found in Tien Shan, Pamir and western Himalayas (Fig. 5)

Temporal distribution of seasonal SWE volume
The temporal variations in integrated seasonal SWE volume across the major river basins are quantified and examined from 295 the 18-year average SWE seasonal cycle (i.e., daily time series shown in Fig. 6 and Fig. 7). Note again that the non-seasonal snow and ice mask has been applied when calculating the aggregated SWE volumes. The key statistics of annual peak SWE volumes are summarized for the entire HMA region and each basin in Table 1. Despite the significant literature on seasonal snowpack in this region, quantification of the regional scale SWE volume are more difficult to obtain, partly due to the large uncertainties in SWE estimation over this region. 300

Northwestern (NW), (c) Southeastern (SE), and (d) Northeastern (NE) basin totals.
Over the record examined, the HMA-wide annual peak SWE volume ( Fig. 6a and Fig. 7a) is found to be largest in WY2005 with a value of 227.12 km 3 , smallest in WY2001 with a value of 114.10 km 3 , and an 18-year climatological mean value of 162.58 km 3 . The climatological peak SWE volume was further assessed in each subregion, and compared against that over the entire HMA (Table 1; Fig. 7). The results show the highest peak SWE volume occurs in NW basins (107.42 km 3 , ~ 66% 310 of domain-wide total), followed by SE basins (29.49 km 3 , ~18%), and NE basins (14.78 km 3 , ~ 9%), which is coherent with the spatial pattern shown in Fig. 3a. Note that around ~7% of HMA-wide SWE volume falls in the regions outside of the watersheds examined (mainly in the northmost regions shown in Fig. 1), which is why these basin-wide quantities do not sum up to 100% of the HMA-wide totals.
For the NW basins, the maximum amount of SWE volume is found in the Indus basin, followed by Amu Darya and Syr 315 Darya ( Fig. 6b and Fig. 7b). Meanwhile, the peak SWE volume is found to occur earlier and disappear faster in the Syr Darya basin, followed by Amu Darya and Indus, with later melt-out corresponding to larger peak SWE volume magnitudes.
When examining the peak SWE volumes in each basin ( Fig. 6b and Fig. 7b), the Indus was found to have the highest value of 63.97 km 3 in WY2009, and an average of 48.95 km 3 across the 18 years; Amu Darya shows the highest peak SWE volume of 48.31 km 3 in WY2017, and an average peak SWE volume of 37.31 km 3 across all 18 years; Syr Darya shows a peak SWE 320 volume of 21.16 km 3 on average (Table 1).
For the SE basins, higher SWE volumes are found in Ganges-Brahmaputra, followed by Yangtze, Salween and Mekong ( Fig. 6c and Fig. 7c). As shown in Table 1, Ganges-Brahmaputra has an average peak SWE volume of 15.59 km 3 , with a https://doi.org/10.5194/tc-2021-139 Preprint. Discussion started: 17 May 2021 c Author(s) 2021. CC BY 4.0 License. carry-over SWE volume of around 2 km 3 at the end of the WY (Fig. 7c); the Yangtze basin has an average peak SWE volume of 8.02 km 3 , with a maximum of 14.79 km 3 in WY2005; Salween and Mekong have much lower SWE volumes that 325 are generally less than 5 km 3 , partly due to their small basin sizes (Fig. 1).
For the NE basins over HMA, the overall magnitude of SWE volumes is smallest (Fig. 6d and Fig. 7d), with an average peak SWE volume of 8.78 km 3 in Tarim, 3.65 km 3 in Yellow, and 2.35 km 3 in Inner Tibetan Plateau respectively (Table 1). These basins all have relatively large areas, but are mostly snow-free or covered by shallow snow as depicted in Fig. 3a.

Elevational distribution of basin-wide peak SWE volume
Finally, the domain-wide and basin-wide peak seasonal SWE and SWE volume climatology distribution vs. elevation was assessed (Fig. 8), obtained at the DOWYs when domain-wide or basin-wide SWE volumes reach their annual maximum 335 (averaged over 18 years). The HMA-wide domain (Fig. 8a) and each basin (Fig. 8b, Fig. 8c and Fig. 8d)  in Smith and Bookhagen (2018). The non-seasonal snow and ice pixels were removed when calculating peak SWE volume, and its fractional areal coverage within a given elevation band, is computed to assess the relative elevational contributions to total seasonal SWE volume (Fig. 8). Similar work has been done in Smith and Bookhagen (2018) that assessed the SWE 340 trends in different elevation bands over HMA-wide and basin-wide scales, while here we focus more on a quantitative assessment in the SWE volume distribution itself against elevation. When examining the full HMA domain, the seasonal SWE volume was found to be most abundant at mid elevations (3000 -4000 m), with peak SWE values occurring at elevations around 3500 m (Fig. 8a). Meanwhile, the presence of non-seasonal snow and ice becomes evident at elevations above 3500 m, and it increases dramatically above 5000 m with a value up to 350 35% (Fig. 8a). When assessing the cumulative fraction of SWE volume as a function of elevation, it was found that over https://doi.org/10.5194/tc-2021-139 Preprint. Discussion started: 17 May 2021 c Author(s) 2021. CC BY 4.0 License. 50% of seasonal SWE volume is stored at elevations above 3500 m, and less than 10% of seasonal SWE volume is stored at elevations below 2000 m (Fig. 8a).
Within all NW basins (Syr Darya, Amu Darya and Indus; Fig. 8b) as well as Tarim in NE (Fig. 8d), and Ganges-Brahmaputra and Yangtze in SE (Fig. 8c), both the SWE depth and SWE volumes increase with elevation (below ~4000 m) 355 and then decline with elevation (above ~4000 m). This corresponds to an increase in non-seasonal snow and ice coverage (up to 60%) above 4000 m, which may be primarily responsible for the reduction in seasonal SWE volumes seen at high elevations. Other basins in SE (Salween, Mekong; Fig. 8c) and NE (Yellow, Inner Tibetan Plateau; Fig. 8d) generally show monotonically increasing SWE and SWE volumes against elevation. This is partly due to a much lower non-seasonal snow and ice coverage (mostly under 25%) at high elevations in those basins. For the cumulative fraction of SWE volumes above 360 a specified elevation, the median values are found within 3000 to 4000 m in NW basins (Fig. 8b). Meanwhile, the median values are found at higher elevations for SE and NE (between 4000 -5200 m) basins ( Fig. 8c and Fig. 8d).

Conclusions
A first-order spatiotemporal analysis of seasonal SWE over the HMA region is presented in this paper, using a new 18-year snow reanalysis dataset (HMASR; Liu et al., 2021). This HMASR dataset is derived based on a previously developed snow 365 reanalysis scheme ) that jointly assimilates fSCA observations from both Landsat and MODSCAG products, which has daily outputs of SWE and other snow variables, with a spatial resolution of 16 arc-second (~500 m), over the joint period of Landsat and MODIS from WYs 2000 to 2017.
Several scientific questions were addressed in this work using the new HMASR dataset, such as how seasonal SWE and snow storage distribute with space, time and elevation. In terms of the spatial distribution, seasonal snow is most abundant in 370 the NW, with over 1 m of peak SWE observed over the mountain ranges. Seasonal snow is also significant in the SE, where both relatively deep snowpacks (with peak SWE values up to 1 m or above) and shallow snowpacks (with peak SWE up to 0.2 m) are found. Seasonal snow is less abundant in the NE where most areas are snow free, or only covered by shallow snowpacks (with peak SWE values below 0.2 m). The domain-wide median date of peak SWE is estimated to be March 18th with significant heterogeneity across this region, linked with climatological drivers and topography. 375 When aggregating the total SWE volumes across the full HMA domain and its basins, the climatological peak seasonal SWE volume was found to be 163 km 3 , with NW basins accounting for around 66% of that volume, followed by SE (~18%) and NE (~9%) basins. This peak seasonal SWE volume showed a maximum value in WY2005 (227 km 3 ) and a minimum value in WY2001 (114 km 3 ) when assessed over the entire HMA region in the study period. The elevational distribution of seasonal peak SWE volume shows that SWE volume is most abundant at mid-elevations (3000 -4000 m), with over 50% of 380 the seasonal SWE volume stored at elevations above 3500 m, when assessed over the entire HMA region.
This HMASR dataset is presented to augment the spatiotemporal gaps in previous SWE datasets and provide better characterization of spatiotemporal patterns in seasonal snowpack over the HMA region, especially over the mountainous https://doi.org/10.5194/tc-2021-139 Preprint. Discussion started: 17 May 2021 c Author(s) 2021. CC BY 4.0 License. areas with complex terrain where existing products tend to underestimate SWE and present large uncertainties (Wrzesien et al., 2019;Kim et al., 2021). It should prove useful in providing more insight into the role of seasonal snowpack in the 385 regional hydrological cycle, as a verification dataset for atmospheric and other models, and in other applications where a space-time continuous snow dataset constrained by remote sensing data is needed.
A key limitation in the dataset is that, it is most effective in estimating seasonal snowpack, as the method requires a seasonal cycle of fSCA with a clear depletion signal. It is therefore not designed to provide estimates of SWE associated with glaciers or semi-permanent snow, and is generally expected to be less accurate in areas of intermittent and very shallow snowpack. In 390 addition, limited fSCA observations due to cloud contamination during monsoon seasons may result in higher uncertainty in SWE estimation over affected sub-regions like the Himalayas. Other remote-sensing approaches (e.g., active microwave measurements) that could penetrate clouds may potentially aid in reducing the uncertainties for SWE estimation over those areas. More research can be done to address such issues and improve the accuracy of SWE estimates for those regions in the future. 395 Data availability. The HMASR dataset used in this paper, is publicly available on National Snow and Ice Data Center (NSIDC) HiMAT data repository, entitled: 'High-Mountain Asia UCLA Daily Snow Reanalysis, Version 1'. It can be accessed through https://nsidc.org/data/HMA_SR_D/ or https://doi.org/10.5067/HNAUGJQXSCVU (Liu et al., 2021). The dataset is provided as NetCDF files for each 1° x 1° tile shown in Fig. 1, available at 16 arc-second (~ 500 m) and daily resolution from WYs 2000 to 2017. Posterior estimates of other key snowpack properties (i.e. in addition to SWE) not 400 focused on herein (e.g. snow depth, fSCA, snowmelt, sublimation, snow albedo, etc.) along with posterior forcing variables are included in this dataset. Data quality information, containing a classification mask and the non-seasonal snow/ice mask, can be found in the dataset as well.
Author contribution. YL led the production of the dataset and led the analysis in the manuscript. YF contributed to the production of the dataset and contributed to the analysis in the manuscript. SM supervised the project and provided guidance. 405 All authors contributed to writing the manuscript.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. We would like to thank those responsible for the creation of the datasets used in this study as well as members of the NASA High-Mountain Asia Team (HiMAT) who helped shape the overall direction of the work. G. Cortes is acknowledged for pre-processing some of the input data used to develop the HMASR dataset. This research was funded by 410 the NASA High-Mountain Asia Program Grant #NNX16AQ63G with additional support provided by NSF Grant #1641960.