Validation of modeled snow properties in Afghanistan , Pakistan , and Tajikistan 1 2

2 Edward H. Bair1, Karl Rittger2, Jawairia A. Ahmad3, and Doug Chabot4 3 4 1Earth Research Institute, University of California, Santa Barbara, California, USA 5 6832 Ellison Hall, University of California, Santa Barbara, CA 93106-3060. correspondence 6 email: nbair@eri.ucsb.edu 7 8 2Institute for Arctic and Alpine Research, University of Colorado, Boulder, Colorado, USA 9 10 3Department of Civil & Environmental Engineering, University of Maryland, College Park, MD, 11 USA 12

ABSRACT: Ice and snowmelt feed the Indus and Amu Darya rivers, yet there are limited in situ measurements of these resources.Previous work in the region has shown promise using snow water equivalent (SWE) reconstruction, which requires no in situ measurements, but validation has been a problem until recently when we were provided with daily manual snow depth measurements from Afghanistan, Tajikistan, and Pakistan by the Aga Khan Agency for Habitat (AKAH).For each station, accumulated precipitation and SWE were derived from snow depth using the SNOWPACK model.High-resolution (500 m) reconstructed SWE estimates from the ParBal model were then compared to the modeled SWE at the stations.The Alpine3D model was then used to create spatial estimates at 25 km to compare with estimates from other snow models.Additionally, the coupled SNOWPACK and Alpine3D system has the advantage of simulating snow profiles, which provide stability information.Following previous work, the median number of critical layers and percentage of facets across all of the pixels containing the AKAH stations was computed.For SWE at the point scale, the reconstructed estimates showed a bias of -42 mm (-19%) at the peak.For the coarser spatial SWE estimates, the various models showed a wide range, with reconstruction being on the lower end.For stratigraphy, a heavily faceted snowpack is observed in both years, but 2018, a dry year, according to most of the models, showed more critical layers that persisted for a longer period.

INTRODUCTION
There are many parts of the world where little is known about the snowpack.This lack of knowledge presents a challenge for water managers and for avalanche forecasters.Afghanistan is particularly austere in this respect, as there have been no snow measurements available since the early 1980s.This lack of information about the snowpack potentially creates a humanitarian crisis, as snowmelt fed streams run dry in the fall without warning (USAID, 2008).Accurate historical estimates of basin-wide snow water equivalent (SWE) are crucial for creating a baseline of climatological conditions, which can then aid in predicting today's SWE.For example, climatological estimates of spatially-distributed SWE are the most important predictors in machine learning statistical models for this region (Bair et al., 2018b).
To improve our knowledge about the snowpack in these areas, we have developed an approach that requires no in situ measurements.Using satellite-based estimates of the fractional snowcovered area (fSCA) and downscaled forcings in an energy balance model, we build up the snowpack in reverse, from melt out to its peak, using a technique called SWE reconstruction (Martinec and Rango, 1981).This technique has been shown to accurately estimate SWE in mountain ranges across the world, including: the Sierra Nevada USA (Bair et al., 2016;Rittger et al., 2016); the Rocky Mountains USA (Jepsen et al., 2012;Molotch, 2009); and the Andes of South America (Cornwell et al., 2016)-all areas with relatively abundant independent ground validation measurements.For the so called Third Pole of High Mountain Asia, and especially the northwestern parts of this region, e.g.Afghanistan, Tajikistan, and Pakistan, ground-based validation is challenging.

AGA KHAN AGENCY FOR HABITAT (AKAH) STATIONS
In 2017, we received daily manual snow depth and other meteorological measurements from nearly 100 stations (Figure 1) in an operational avalanche network (Chabot and Kaba, 2016).These stations are funded by the Aga Khan Agency for Habitat (AKAH) and are the first snowpack measurements available, at least that we are aware of, in Afghanistan in nearly 40 years.Hence, we refer to the region as the AKAH study region and the weather stations as the AKAH stations.The AKAH stations contain manual daily snow depth (also called height of snow), height of new (24-hr) snow, daily high and low air temperature, instantaneous wind speed/direction, rainfall, and some text fields on weather and avalanche conditions.For mountainous areas, precipitation is the most uncertain term in the water balance (Adam et al., 2006;Milly and Dunne, 2002) because it exhibits high spatial variability and is difficult to measure with traditional gauges.Measuring snow on the ground has many advantages compared to using precipitation gauges, which suffer from undercatch, especially in the windy and treeless areas (Goodison et al., 1998;Kochendorfer et al., 2017;Lehning et al., 2002a) typical of this part of the world.Likewise, a strength of the SWE reconstruction technique is that it does not depend on precipitation measurements to build the snowpack.
Additionally, many of the AKAH stations are at high altitudes, with 64 stations above 2500 m and 17 stations above 3000 m.Unfortunately, most of these stations are located in deep valleys, where the villages are, rather than on exposed mountain sides or ridges and the daily resolution is too coarse to use in a snow model without temporal interpolation.Additionally, many of the stations are near glacierized areas which complicates spatially interpolated snow estimates, as some of the snow is on top of ice.The area covered by glaciers in Figure 1 is 7.8%.Asia, there are fewer studies examining snowfall in High Mountain Asia, which is odd since hydrologically in this region, snow on land melt provides the vast majority of runoff compared to snow on ice and melting glacier ice (Armstrong et al., 2018).Many of these studies are focused on the region to the east of the AKAH study area shown in Figure 1.To our knowledge, there have been no studies on snowpack stratigraphy in the AKAH study area and we were unable to obtain any snow pit measurements from this area.

LITERATURE REVIEW
A few studies have specifically examined snowfall in larger regions that include some of the AKAH stations, mostly for stations in the southern basins that flow into the Indus River; that is all of the stations in Pakistan.The rest of the stations in Afghanistan and Tajikistan flow into the Amu Darya River.The most comparable study (Shakoor and Ejaz, 2019) examines the Passu catchment in the Hunza River Basin, to the east of Figure 1.As in this study (Section 5.1), Shakoor and Ejaz (2019) also use the SNOWPACK and Alpine3D models.Model parameters were calibrated using a single weather station, Urdukas at 3926 m elevation near the Baltoro glacier (Ev-K2-CNR, 2014), with one year of precipitation measurements, using snow depth for validation.The authors report overestimation of the measured snow depth at the calibration station, even after questionable adjustments to the snow albedo and other model parameters.They attribute the overestimation to problems with the precipitation measurements, common for high altitude stations.One problem with the Urdukas station in particular is that the tipping bucket precipitation gauge is unheated, making it unusable for measuring solid precipitation.Temperatures at this station were well below freezing for the winter and most of the spring, which explains why no precipitation was recorded from January until sometime in March during 2012, the calibration year.Viste and Sorteberg (2015) study several gridded precipitation products throughout High Mountain Asia, including the Indus River Basin.They report that while total precipitation was similar across the products-including MERRA (Rienecker et al., 2011), APHRODITE, TRMM (Huffman et al., 2007), and CRU (Harris et al., 2014)-the total snowfall varied by 2 to 4 ×.Smith and Bookhagen (2018) used 24 years (1987to 2009) of satellite-based passive microwave SWE estimates to examine trends throughout High Mountain Asia, including the Amu Darya and Indus Basins.Their SWE estimates show most 25 km pixels in this region in the 50-100 mm range for December through February, with a few over 100 mm in the Amu Darya (i.e.all the AKAH stations in Afghanistan and Tajikistan) and none over 100 mm in the Indus (i.e.all the AKAH stations in Pakistan), likely too low by an order of magnitude for some pixels given our previous reconstructed SWE values and limited climate measurements in Afghanistan (Bair et al., 2018b).
For the AKAH stations in Tajikistan, the most comprehensive snow measurements come from Soviet snow surveys (mostly depth, but with some SWE and density measurements) that have been digitized (Bedford and Tsarev, 2001).Most of these measurements begin in the late 1950s and end around the fall of the Soviet Union, in either 1990 or 1992, making them useful for climatological studies, but not for validation of modern satellite-based estimates.
The sole source of snow measurements in Afghanistan that were accessible to us was a table of outdated WMO monthly climatological data from Kabul (el.1791 m) and North Salang (el.3366 m), showing the maximum monthly snow depth and the mean number of days with snow (Table 1 in Bair et al., 2018b).Again, these measurements are not useful to validate more modern snow estimates.
There have been many other studies that have attempted to estimate basin-wide precipitation (including snowfall) for larger areas that include the AKAH region, especially in the Indus.Several climate studies of the Indus have focused on using lower elevation precipitation gauges, which are then used to spatially interpolate basin-wide precipitation.Dahri et al. (2016) and Dahri et al. (2018) have assembled perhaps the largest collection of climatological measurements covering the AKAH region, mostly based on gauge measurements, as part of a study on the hydrometeorology of the Indus Basin.Using undercatch corrections based on wind, often from reanalysis, they increased precipitation estimates by 21% on average throughout the Indus Basin (Dahri et al., 2018).For example, in the Gilgit sub-basin, they find an unadjusted precipitation estimate of 582 mm/year, adjusted to 787 mm/year, a 35% increase.Although some of the measurements are taken from publicly available sources, as with most publications for this region, the comprehensive data used are not publicly accessible.A similar but less sophisticated approach was used by Lutz et al. (2014), who used a constant increase of 17% across the APHRODITE (Yatagai et al., 2012) precipitation dataset which covers all of High Mountain Asia.Immerzeel et al. (2015) used glacier mass balance estimates with streamflow measurements as validation to show that high-altitude precipitation in the upper Indus Basin is 2 to 10 × what is shown using gridded precipitation products like APHRODITE.Bookhagen and Burbank (2010) estimate that snowmelt contributes 66% of annual discharge to the Indus, and averages 424 mm across the basin.
In summary, quite a few studies have produced varying precipitation and snowfall estimates for the AKAH region, with no recent in situ snow measurements from Afghanistan or Tajikistan.

PREVIOUS WORK WITH AKAH SNOW MEASUREMENTS
Our previous work (Bair et al., 2018b) used a simple density model (Sturm et al., 2010) based on snow climatology (Sturm et al., 1995) and day of year to model SWE from the manual snow depth measurements.The density model itself has -12 to 26% bias in predicting SWE.When taking into account geolocational uncertainty of the reconstructed SWE estimates and uncertainty in the density model, errors are on the order of 11-13% Mean Absolute Error (MAE) and -2 to 4% bias, depending on the date.However, we only examined one year of the AKAH station data (2017) and the high uncertainty in the density model itself begs a more sophisticated approach.
From recent work (Bair et al., 2018a), we have shown that the SNOWPACK (Bartelt and Lehning, 2002;Lehning et al., 2002a;Lehning et al., 2002b) model is capable of accurate SWE prediction when supplied only with snow depth for precipitation, as well as the other requisite forcings (i.e.radiation, snow albedo, temperatures, and wind speed).Over a 5-year period using hourly in situ measured energy balance forcings and a snow pillow for validation at a high elevation site in the western US, the SNOWPACK modeled SWE showed a bias of -17 mm or 1% (Bair et al., 2018a).

METHODS
Our modeling approach consisted of: a) downscaling forcings in ParBal and reconstructing SWE; b) combining the downscaled forcings for each AKAH station with temporally interpolated manual snow measurements; c) running SNOWPACK for each of the AKAH stations with the downscaled and interpolated measurements from a) and b); and d) running Alpine3D using the output from SNOWPACK, notably the hourly precipitation.In addition to predicting SWE, the SNOWPACK/Alpine3D coupled model also predicts stratigraphic parameters useful for avalanche prediction, thereby giving us an idea of the layering and stability in this region.For comparison, we also ran the NOAH-MP land surface model over the region with widely-used forcings.We also compared spatial estimates of SWE from GLDAS-2.Methods are summarized in Table 1 and explained below, with more detail provided in Appendix A.

SNOWPACK and Alpine3D
SNOWPACK and Alpine3D are freely available (https://models.slf.ch)point and spatially distributed snow models, courtesy of the Swiss Federal Snow Institute.SNOWPACK is the older of the two and uses finite elements to model all of the layers in a snowpack at a point.SNOWPACK has shown promising results in both operational (e.g.Lehning et al., 1999;Nishimura et al., 2005) and research applications (e.g.Bellaire et al., 2011;Hirashima et al., 2010).Alpine3D (Lehning et al., 2006) is essentially a spatially-distributed version of SNOWPACK with a number of additional modules including: terrain-based radiation modeling, blowing snow, and hydrologic modeling.Integral to Alpine3D is SNOWPACK, which is run for each pixel, as well as the MeteoIO library (Bavay and Egger, 2014), which provides a large number of temporal and spatial interpolation functions that can be used on forcings for Alpine3D and SNOWPACK.

The Parallel Energy Balance Model
The Parallel Energy Balance Model (ParBal) was created at UC-Santa Barbara and designed for reconstruction of SWE.It is also publicly available (https://github.com/edwardbair/ParBal/releases/tag/v1.0).Currently, ParBal is designed to use: downscaled temperature, pressure, and humidity from version 2 of the Global or National Land Data Assimilation System (GLDAS-2/NLDAS-2, Rodell et al., 2004;Xia et al., 2012); shortwave and longwave radiation from edition 4a of the Clouds and the Earth's Radiant Energy System (CERES, Rutan et al., 2015) SYN product; and time-spaced smoothed (Dozier et al., 2008;Rittger et al., in press) snow surface properties from MODIS Snow Covered Area and Grain Size (MODSCAG, Painter et al., 2009) and MODIS Dust and Radiative Forcing in Snow (MODDRFS, Painter et al., 2012).ParBal is run hourly at 500 m spatial resolution and forcings are adjusted for terrain and elevation.The main output is the residual energy balance term, which is assumed to go into melt when positive during the ablation phase after cold content is overcome (Jepsen et al., 2012).This residual melt term is then summed in reverse during periods of contiguous snow cover and multiplied by the fSCA to spread the snow spatially.Details on ParBal are covered extensively in Bair et al. (2016) and Rittger et al. (2016).

NOAH Multi-Parameterization (MP)
The NOAH-MP v3.6 (Ek et al., 2003;Niu et al., 2011) land surface model, forced using MERRA-2 (Gelaro et al., 2017), was used to simulate the hydrologic cycle over the study area and provide SWE estimates for comparison with ParBal and the Alpine3D output.NOAH-MP was selected due to its detailed representation of the snowpack relative to other land surface models.The model subdivides the snowpack into up to three layers with associated liquid water storage and melt/refreeze capability (Niu and Yang, 2004;Yang and Niu, 2003).It incorporates the exchange of heat and moisture through the snowpack between the land surface and the atmosphere.In a model intercomparison study using a 2 km spatial resolution regional climate model for forcings, Chen et al. (2014) show that NOAH-MP modeled peak SWE at SNOTEL sites in Colorado, USA with a -7% bias.

Use of AKAH station measurements
We modeled daily SWE at the AKAH stations during the 2017 and 2018 water years primarily using the manually measured height of snow (HS), also called snow depth, combined with our downscaled energy balance parameters (for downscaling methodology see Bair et al., 2018b;Bair et al., 2016;Rittger et al., 2016).We choose the manual HS and new (24-hr) snow (HN) as the only variables to use from the AKAH stations.The HS appeared to be the most reliably measured, as that only requires reading a value from a master snow depth stake.The HN was used to correct a data entry problem in 2017 that we discuss below.The reliability of the other measurements (instantaneous wind speed/direction, maximum/minimum temperature, and rainfall) was questionable.For example, we were not provided with sensor or measurement metadata, e.g.sensor make/model, measurement height, and whether or not the temperature sensor was shielded from shortwave radiation.These other measurements taken daily were also of limited value for interpolation to hourly values (see item 3 below).Thus, these other measurements were not used.
The AKAH dataset had a number of shortcomings that we list here along with how we addressed them.
1) Some of the stations recorded no snow at all, especially in the dry 2018 year, or had obvious problems, so they were excluded.For 2017, 52 (54%) of stations were used.For 2018, 41 (46%) stations were used.
2) There were spurious drops in the HS measurements.The drops were clearly cases of missing values being filled with zeros.These measurements were manually flagged and converted to null values for interpolation, see below.
3) The daily measurements had to be interpolated to hourly values.For the most part we used linear interpolation, although this is not ideal during snow accumulation since it's almost never the case that snowfall is uniform over a 24-hr period.This is a problem that affects the accuracy of snow settlement estimated by SNOWPACK.There were two cases where other interpolation methods were used.If there were several days of missing values, we used a nearest neighbor interpolation to fill in the missing daily values, followed by a linear interpolation from daily to hourly measurements such that we assumed all the new snow fell in a 24-hr period.The other case was for days where the linear interpolation would yield a value below the minimum threshold hard coded into SNOWPACK (0.5cm/hr) for the first accumulating snowfall on bare ground.In this case, a previous neighbor interpolation was used in such a way that the entire snowfall occurred in the last hr prior to the next day's measurement.
4) We found the AKAH stations suitable for snow on the ground measurements, but not for rainfall or total (solid+liquid) precipitation.This was only an issue for the Alpine3D snow modeling, as snow measurements were being extrapolated to higher elevations than the AKAH stations (Section 6.2), thus at these higher elevations, snow accumulated earlier and melted later than at the lower AKAH stations.
Given the near total lack of canopy cover in the region, we suspected substantial undercatch from rain gauges.Using the wind speed, an undercatch correction would have been possible given more information on the gauges (e.g.orifice opening diameter and whether or not a shield was present), however this instrument metadata was not available to us.Likewise, we did not know if the gauges were heated or not.
Further, the time period for recording measurements from the stations was not consistent.In WY 2017, measurements began being reported on 10 November 2016 and were reported until 24 November 2017.However, in WY 2018, measurements weren't reported until 1 December 2017 and no station measurements were reported past 1 April 2017.The reporting period likely covered all the snowfall events, but not all the precipitation events.
To address the rainfall measurement and reporting issues, we used GLDAS NOAH v2.1 (Rodell et al., 2004) rainfall + snowfall from the nearest grid cell (1/4º spatial / 3 hr temporal resolution) to fill in precipitation prior to the first measurements in each water year, and after 4-1 for both water years.We did not account for rain from 10 November 2016 to 1 April 2017 and from 1 December 2017 to 1 April 2018; instead we relied on the modeled precipitation from SNOWPACK using snow depth.The AKAH station observations show that rain during this time period was rare.
5) A database problem prevented snow heights > 100 cm from being entered into the database for a few days in 2017.This problem became apparent during February 2017, when the Nuristan avalanches took place, as that is the first time that most stations recorded values > 100 cm.Values were shown as 100 cm on multiple days followed by values > 100 cm.To address this issue, we flagged all the values equal to 100 cm prior to peak depth in 2017, then marked those as null values.We then filled those null values using the cumulative sum of new snow during that time.

Analysis of modeled snow profiles
For holistic measures of the snow profiles modeled in Alpine3D, we used two metrics from Seasonal Snow on the Ground primary codes FC, DH, and SH (Fierz et al., 2009), divided by the height of the snowpack.The number of critical layers was computed using a threshold sum approach (Schweizer and Jamieson, 2007) with modifications for simulated profiles (Monti et al., 2014 Table 1).In each profile, 6 different variables (grain size, difference in grain size, hardness, difference in hardness, grain type, and depth) in the top meter of height (from the surface) were checked against threshold values.Layers exceeding 5 or more thresholds were classified as critical.
The fraction of facets metric does not have a validation study, but faceted layers are a weak crystal form and are responsible for 43% (Bair et al., 2012) to 67% (Schweizer and Jamieson, 2001) of investigated avalanches.Layers classified as critical using the threshold approach above corresponded to failure layers about ½ of the time (Monti et al., 2014) in Compression Tests (Jamieson, 1999;van Herwijnen and Jamieson, 2007), an in situ snowpack stability test.

Spatial scale for comparisons
Because ParBal is the only model run at 500 m spatial resolution and all the other models were run at ~ 25 km, it is the only model appropriate for point comparisons, although point to area problems are still an issue.To address the geolocational uncertainty for the gridded MODIS products, which can be up to one ~500 m pixel (Tan et al., 2006;Xiaoxiong et al., 2005) and spatial variability of the snow, we used a 9-pixel neighborhood centered on each AKAH station and chose the best fit to the SNOWPACK modeled SWE.This approach has been used in previous work (Bair et al., 2018b;Rittger et al., 2016).
For all of the other model comparisons, we resampled all of the model output to a UTM (Zone 43S) grid with 25 km pixels, close to the native resolution of the NOAH-MP and GLDAS2 grid used (0.25º).The ParBal output had to be significantly upscaled from 500 m to 25 km using Gaussian Pyramid reduction (Burt and Adelson, 1983) in steps with bilinear interpolation for the final step.

RESULTS AND DISCUSSION
The relationships between the components are summarized in Figure 2. The results discussed below are comparisons of: 1) SWE and 2) snow stratigraphy across a) all of the AKAH stations (points) and b) the entire study region.

Point comparisons between SNOWPACK and reconstructed SWE 308
A first step for any SWE reconstruction comparison is to determine when the ablation season 309 starts.This varies for different years and at different sites (e.g.Margulis et al., 2016).Using the 310 SNOWPACK modeled SWE, we can examine the peak SWE dates for both years for all of the 311 AKAH stations (Figure 3ab).Peak SWE dates vary across the stations and years, but the median 312 values are between years are a week apart, 19 February 2017 and 26 February 2018.Thus, we 313 use those dates for our comparisons.314  were computed and plotted for each day during the ablation season (Figure 4).For the reconstructed SWE on 19 February 2017, the bias is -77 mm (-28%).For the reconstructed SWE on 26 February 2018, the bias is -6 mm (-9%).Thus, together these biases average to -42 mm (-19%).The increases in reconstructed SWE during the ablation season are caused due to differences in how melt is summed for any given pixel.In ParBal, melt is only summed during periods of contiguous snow cover.This means that if a pixel containing an AKAH station has no snow on it at some point during the ablation season, but then snow is detected, it causes an increase in the mean SWE.This is called an ephemeral snow event, i.e. snow that disappears and reappears.For a more in depth examination of the error at individual stations, a box plot is shown for the median peak SWE dates for both years (Figure 5).The median bias of the reconstructed SWE is -11 mm (-14%).

Four model spatial comparisons
The AKAH stations are lower than the average elevation for the region.The average elevation of the AKAH stations is 2619 m (1735 to 3410 m).But when the 500 m DEM is upscaled to 25 km, the average elevation of the pixels containing the AKAH station is 3858 m with a range of 2517 to 4764 m.This has two important implications: 1) much of the higher elevation snowfall is being extrapolated and 2) the higher elevation causes the peak SWE date to move forward in time.The median peak SWE dates for the (N=169) 25 km pixels encompassing the study area are 5 May 2017 and 3 May 2017.Thus, we use the median of the two to compare our reconstructed SWE values (Figure 6ab and Figure 7a-d).
Striking is the range between models.NOAH-MP has the highest peaks ( 562  the melt if a pixel were 100% snow covered), and the melt from glacierized areas (light blue in 357 Figure 1).We did not find any errors in the model, its parameters, or its forcings.Thus, it is 358 possible that the ParBal SWE is low-biased in 2017 for reasons that we could not discern, or that 359 the other models are high biased.Of note is that the 2017 & 2018 SCA (Figure 8 purple and 360 orange) is very similar for both years during the ablation period, especially at the end of the ablation season.Since pixels do not contribute uniformly to melt, SCA alone cannot be used to predict SWE, but in general years with less snow have lower SCA values towards the end of the ablation season.Figure 8 shows that 2017 and 2018 were similar in terms of SCA from April through melt out.Thus, the large difference between 2017 and 2018 for the AKAH station SWE, but small differences in SCA and spatially-averaged reconstructed SWE, suggest that 2017 may have been a larger snow year at the lower elevations where the AKAH stations are, but similar to 2018 at the higher elevations.

Stratigraphy and stability
The simulated snow profiles from the AKAH stations (Figure 9ab) and the 25 km pixels containing the AKAH stations (Figure 10ab) show very different snowpacks.Because of the induced increase in elevation from scaling (e.g. from an average of 2619 m to 3858 m, Section 6.2), the 25 km pixels show a deeper, but more faceted snowpack with critical layers that persist for a month or longer.In 2017, for the median AKAH station values, the snowpack reaches a maximum of 76% facets on January 21 (Figure 9a).In 2018, the snowpack reaches a maximum of 71% facets (Figure 9b).There were no critical layers simulated.In contrast, for the median values in the 25 km pixels for both years, the height of snow (HS) is approximately 2 × that for the stations (Figure 10ab).The snowpack reaches a maximum of 94% facets in 2017, with one critical layer persisting for 35 days (Figure 10a).The snowpack in 2018 reaches 95% facets with 1 or 2 critical layers persisting for 80 days (Figure 10b).During the Nuristan avalanches on 4 February to 7 February 2017 that killed over 100 people (United Nations, 2017), the AKAH stations show the largest 3-day snowfall of the study period (Figure 9a) and the results for the 25

CONCLUSION
Thanks to a novel operational avalanche observation network, there are now daily snow measurements at a number of operational weather stations in an austere region of High Mountain Asia.In this study, two years of daily snow depth measurements from these stations were combined with downscaled reanalysis and remotely-sensed measurements to force a point and spatially distributed snow model.Compared to a previous effort (Bair et al., 2018b), this study represents a substantial improvement in SWE modeling for the region, and a first attempt to characterize region-wide snow stratigraphy.At the point scale, SWE estimates from a reconstruction technique that does not use precipitation or in situ measurements compared favorably.At the regional scale, four models showed a wide spread in both peak SWE and melt timing.The simulated profiles showed very different snowpacks.At the point scale at lower elevations in the valleys, profiles showed fewer facets and almost no critical layers, while at the regional scale for higher elevations, the profiles showed heavily faceted snowpacks with critical layers that persisted throughout the winter and spring.

PARBAL
ParBal was configured and forced as described in Bair et al. (2018b); Bair et al. (2016).The model time step was 1 hr.The DEM used was the ASTER GDEM version 2 at 1 arc sec (NASA JPL, 2011), while the canopy type and fraction were taken from the Global Land Survey at 30 m (USGS, 2009).The shortwave and longwave forcings were downscaled from the CERES SYN edition 4a 1º/1 hr product (Rutan et al., 2015), while the air temperature, specific humidity, air pressure, and wind speeds were downscaled from the GLDAS NOAH version 2.1 0.25º/3 hr product (Cosgrove et al., 2003).Time-space smoothed (Dozier et al., 2008;Rittger et al., in press) fSCA and grain size from MODSCAG (Painter et al., 2009) was combined with the visible albedo degradation from dust in MODDRFS (Painter et al., 2012) to produce snow hourly snow albedo.

NOAH-MP
NOAH-MP v3.6 was run in retrospective mode within the NASA Land Information System (LIS) framework.A state vector ensemble (total 30 replicates) was generated by perturbing the forcings to account for the state uncertainty during forward propagation of the model.MERRA-2 (Gelaro et al., 2017) forcings were utilized with bilinear spatial and linear temporal interpolation.The model was run on an equidistant cylindrical grid with 0.25º spatial resolution and a 15 min model timestep.The spin-up time extended from May 2002 to May 2016 while the study period was from June 2016 to October 2018.The number of maximum layers in the snowpack was 3.  Plots were assumed to be level, so forcings without terrain correction were applied except for shading when the sun was below the local horizon, e.g. a mountain blocking the sun (Dozier and Frew, 1990).The wind direction, which is not available in GLDAS-2, was fixed at the mean value from the daily AKAH instantaneous values.The ground temperature was set as the minimum of the air temperature or -1.5ºC when snow cover was present.
Aside from setting required parameters and values for inputs and outputs, changes to default parameters that affected model output are provided in

ALPINE3D
Alpine3D version 3.10 was run using with the outputs produced by SNOWPACK as forcings for each of the AKAH stations at 25 km resolution.The DEM and land cover (incorrectly labeled land use in the Alpine3D documentation) data were upscaled from the ParBal data.Alpine3D was run at an hourly timestep using hourly forcings, with daily outputs using the "enable-eb" switch.Other switches were set to off, the defaults.The "enable-eb" switch computes the terrain radiation with shading and terrain reflections (see Alpine 3D documentation at https://models.slf.ch for a description).
To extend the length of the model runs, for each AKAH stations, GLDAS-2 precipitation was appended to periods prior to the first AKAH observation for the year and after the last, as described in Section 5.4.
The forcings were hourly: incoming shortwave, incoming longwave, air temperature, relative humidity, wind speed, wind direction, reflected shortwave, accumulated precipitation, and ground temperature.
Critical to Alpine3D are the interpolation methods from MeteoIO to spatially distribute precipitation and other forcings.We found the modeled SWE to be highly dependent on the spatial interpolation of precipitation.Our initial approach was to explore local (i.e. with a given radius from a station) and regional (i.e.all AKAH stations) lapse rates in the measured snow depth and modeled precipitation from SNOWPACK.We found almost no correlation in many of the measurements, not surprising given the complexity of the terrain and likely existence of microclimates with substantial influence on precipitation.Without having a good validation source for spatial precipitation (as is the case for all of High Mountain Asia), we selected an interpolation method that yielded relatively smooth results, but showed increases in precipitation with elevation.
Ultimately, we decided to use an inverse distance weighting scheme with elevation detrending (IDW_LAPSE) and a multilinear option.For this method, the input data are detrended, then the residuals are spatially interpolated according to an inverse distance weighting scheme.The detrending uses a multiple linear regression with northing, easting, and altitude.The linear regression has an iterative method for removing outliers.Finally, values at each cell are retrended using the multiple linear regression and added to the interpolated residuals.
A summary of the interpolation methods, all of which are defined in the MeteoIO documentation (Bavay and Egger, 2014), is given in Table A3.

Spatial interpolation method Description and notes
Air temperature IDW_LAPSE Inverse distance weighting with elevation detrending.
https://doi.org/10.5194/tc-2019-150Preprint.Discussion started: 16 July 2019 c Author(s) 2019.CC BY 4.0 License.Although there have been a large number of studies examining the glaciers of High Mountain

Figure 1
Figure 1 Study region with AKAH stations (green dots) overlaid on a MODIS true color image from 13 April 2018.Also shown are the country boundaries (red) and glacierized areas (light blue) from the Global Land Ice Measurement from Space dataset(Raup et al., 2007).All of the stations in Afghanistan and Tajikistan eventually flow into the Amu Darya River.All of the stations in Pakistan eventually flow into the Indus River.

Figure 2
Figure 2 Summary of relationships between the various components.Forcings are shown in red, models and selected outputs are shown in blue, and the comparisons discussed below are shown in green.The black arrows show the direction of inputs.
/doi.org/10.5194/tc-2019-150Preprint.Discussion started: 16 July 2019 c Author(s) 2019.CC BY 4.0 License.To create a holistic comparison for all the stations across the ablation period, mean SWE values

Figure 3 Figure 4
Figure 3 Peak SWE dates, modeled by SNOWPACK for 2017 (a) and 2018 (b) for each of the AKAH stations.The median peak SWE dates are 19 February 2017 and 26 February 2018.N=52 and 41 AKAH stations used for 2017 and 2018.

Figure 5
Figure 5 Bias (a) and relative bias (b) error for ParBal reconstructed SWE vs Alpine 3D modeled SWE at AKAH stations the median peak SWE date for both years, where bias here is ParBal SWE -Alpine 3D SWE.
mm in 2017 and 331 mm in 2018), but is among the first to melt out.The reconstructed SWE from ParBal only shows minor variation between the 2017 peak (240 mm) and the 2018 peak (206 mm).ParBal and GLDAS-2 melt snow out latest in both years.The Alpine 3D model shows the second highest peak SWE in 2017 (469 mm), but the lowest peak (165 mm) in 2018.The comparatively higher values from NOAH-MP could result from relatively high precipitation estimates from its MERRA2 precipitation forcings.Similarly,Viste and Sorteberg (2015) report that MERRA (version 1) showed higher snowfall in the Indus Basin than any other reanalysis or observationbased forcings dataset.

Figure 6
Figure 6 Time series of mean SWE for four snow models across the study area (13x13x25 km pixels) show in Figure 1 for 2017 (a) and 2018 (b).The reconstructed SWE from ParBal (yellow) goes back to 4 May, the median peak SWE date for both years, since reconstruction is only valid during the ablation season.

Figure 7
Figure 7 Four model (a-d) spatial comparison for the region on 4 May 2018.The white letters are the 3 letter ISO country codes (AFG-Afghanistan; TJK-Tajikistan; PAK-Pakistan).Also shown in (a) are the locations of the AKAH stations (orange points).This is a frame from a video sequence available as supplementary material.

Figure 8 Figure 9
Figure 8 Time series of snow covered area from spatially and temporally interpolated MODSCAG (Rittger et al., in press), an input for ParBal, for four selected years across the region.Years 2008 and 2009 had the lowest and highest values on July 1 over the period of record from 2001 to 2018, while 2017 and 2018 comprise the AKAH station study period.

Figure 10
Figure 10 Stratigraphy summary of the (13x13) 25 km pixels containing AKAH stations for 2017 (a) and 2018 (b).Plotted are the median: height of snow (HS); fraction of the snowpack containing facets; and number of critical layers.
Table A1 Noah-MP v3.6 physical parametrization scheme options utilized in this study.SNOWPACKSNOWPACK v3.50 was run in research mode at a 15 min timestep with hourly outputs for each of the AKAH stations.Hourly forcings were computed by combining temporally interpolated snow depth from the AKAH manual measurements with: air temperature, incoming shortwave, reflected shortwave, incoming longwave, wind speed, and relative humidity from the downscaled ParBal outputs, as described in Section 5.2.SNOWPACK was only run for periods when measurements from the AKAH stations were available, Nov/Dec to April/May, depending on the year.