Monitoring Spatial and Temporal Differences in Andean Snow Depth Derived From Satellite Tri-Stereo Photogrammetry

Quantifying the high elevation winter snowpack in mountain environments is crucial for lowland water supply, though it is notoriously difficult to accurately estimate due to a lack of observations and/or uncertainty in the distribution of meteorological variables in space and time. We compare high spatial resolution (3 m), satellite-derived snow depth maps for two drought years (2017 and 2019) in a high mountain catchment of the central Chilean Andes, applying a recently updated methodology for spaceborne photogrammetry. Regional weather station observations revealed an 80% reduction in precipitation for 2019 (the second driest winter since 1950) relative to 2017, though only a 10% reduction in total snow-covered area is seen in the satellite imagery. We threshold surface height changes based upon uncertainty of stable (snow-free) terrain differences for topographic characteristics of the catchment (slope, aspect, roughness etc). For a conservative analysis of change, outside of the topographically-derived confidence intervals, we calculate a mean 0.48 ± 0.28 m reduction of snow depth and a 39 ± 15% reduction in snow volume for 2019, relative to 2017 (for 23% of the total catchment area). Our findings therefore quantify, for the first time in the Andes, the relationship of high-resolution mountain snow depth observations with low elevation precipitation records and characterise its inter-annual variability over high elevation, complex terrain. The practical use of such detailed snow depth information at high elevations is of great value to lowland communities and our findings highlight the clear need to relate the high spatial (Pléiades) and temporal (in-situ) scales within the available datasets in order to improve estimates of region-wide snow volumes.


INTRODUCTION
Seasonal snow is a crucial fresh water resource in mountain regions and plays a highly important role in the socio-economic well-being of lowland communities for drinking water, agriculture, hydropower and mining (Mankin et al., 2015;Meza et al., 2015;Arheimer et al., 2017;Sturm et al., 2017;Biemans et al., 2019). Despite this, quantifying the spatio-temporal variability of the winter snowpack remains difficult and/or costly (e.g., Bühler et al., 2016;Painter et al., 2016;Möller et al., 2017), often relying on intensive ground and/or airborne surveys (Harpold et al., 2014;Deems et al., 2015;Nolan et al., 2015;Currier et al., 2019). Our understanding of large scale snow processes has been greatly aided by decades of global satellite observations of snow cover extent (e.g., Hall et al., 2010) and more recently of snow volume by leveraging highly detailed airborne LiDAR surveys (Painter et al., 2016;Hedrick et al., 2018;Margulis and Fang, 2019). Studies have found that improving representation of snow water equivalent (SWE) from such detailed snow depth observations can increase the predictive capability of seasonal streamflow forecast models (Li et al., 2019;Margulis and Fang, 2019) and help to quantify the inter-annual variability and persistence of the snow depth in space (Hedrick et al., 2018). Nevertheless, the spatial extent and cost of such intensive campaigns, despite their invaluable contribution to process understanding, limits the wider applicability to catchments in many parts of the world where snow water resources are highly important, yet more poorly understood (e.g., Cortés et al., 2016;Fayad et al., 2017;Smith and Bookhagen, 2018;Baba et al., 2018).
The advent of high resolution, optical satellite constellations has lead to many valuable contributions in the discipline of Earth sciences (Zhou et al., 2015;Bagnardi et al. 2016), including quantification of glacier mass balance (Berthier et al., 2014;Melkonian et al., 2016;Belart et al., 2017;Brun et al., 2018;Bƚaszczyk et al., 2019;Ren et al., 2020) and spatial snow volumes (Marti et al., 2016;McGrath et al., 2019;Deschamps-Berger et al., 2020;Shaw et al., 2020a). The application of sub-metre resolution stereo imagery for elevation model generation, namely that of Pléiades and WorldView products, has demonstrated decimetre accuracy when comparing to high resolution, albeit typically small scale, reference data (Marti et al., 2016;Shaw et al., 2020a). Recent work by Deschamps-Berger et al. (2020) found a sub-metre random error and low bias of Pléiades-derived snow depths when comparing to the Airborne Snow Observatory ("ASO") over an extensive (137 km 2 ) area of the Tuolumne Basin, United States. The direction of this research suggests that the use of such high spatial and temporal resolution (tri-) stereo satellite imagery will continue to be a valuable tool for understanding the mountain snowpack, regardless of the limitations that it presents, notably that of cloud cover occlusion and a reduced accuracy for steep terrain (Deschamps-Berger et al., 2020;Shaw et al., 2020a).
This study focuses upon the application of Pléiades-derived snow depths for two distinct years in a high mountain catchment of the central Chilean Andes. The winter snowpack in the central Andes is a vital resource for lowland communities as it offers a storage of water to be released during the spring/summer when typically <15% of the annual precipitation occurs (Garreaud et al., 2009). The very high (up to 5,000-6,000 m) elevation differences between the mountain peaks and populous lowland regions further emphasise the strong dipole of the source and sink of snow which is exacerbated by prolonged severe drought conditions for the region (Garreaud et al., 2017;Garreaud et al., 2019). Quantifying the snow volume of these high mountain catchments is therefore of great importance (Masiokas et al., 2006;Gascoin et al., 2013;Ragettli et al., 2014;Cortés et al., 2016;Masiokas et al., 2020;Shaw et al., 2020a;Shaw et al., 2020b), yet the lack of high elevation observations is a limiting factor in understanding the spatiotemporal behaviour of the snowpack (Alvarez-Garreton et al., 2018;Scaff et al., 2018;Masiokas et al., 2020). This is particularly evident when there are insufficient observations at very high elevations (>∼3,000 m a.s.l.) (Scaff et al., 2018;Shaw et al., 2020a) or for observing precipitation events that originate from the eastern side of the Andes .
We utilised 2 years of Pléiades imagery to explore the interannual variability of snow depth in a high elevation Andean catchment under prolonged regional drought conditions. Optical high-resolution satellite images were first evaluated in the Andes by Shaw et al. (2020a), though this work aims to extend that initial investigation with recently updated workflows and a previously unexplored inter-annual comparison of Pléiades-derived snow depths and their topographically-varying uncertainties. Specifically, we aimed to i) quantify the mountain snowpack using high spatial resolution (3 m) snow depth maps generated from Pléiades imagery in separate years, ii) explore the spatial uncertainties in Pléiadesderived snow depth and its dependence on the topography of the catchment, and iii) evaluate the relationship between spatial snow depth differences and persistence of snow to topographic parameters of the catchment and the meteorological conditions of the preceding winter months. Considering that this is the first time that such an inter-annual comparison of Pléiades snow depth maps has been performed, we place some emphasis on the steps taken to evaluate multi-year digital elevation model (DEM) uncertainties in order to provide generalised guidelines for future work.

STUDY SITE
The Río Yeso catchment is located in the semi-arid Andes of central Chile (33.44°S, 69.93°W) and is at the headwall of one of the many important tributaries to the Maipo River, which serves as the main freshwater resource for the country's capital, Santiago ( Figures 1A,B). The glacierised catchment has an area of 102 km 2 with an elevation range of ∼2,900-5,400 m a.s.l. ( Figure 1C). There are three main glaciers in the basin: Bello (4.6 km 2 ), Yeso (2.9 km 2 ) and Pirámide (4.7 km 2 , debriscovered) which have been estimated to contribute between 3-32% of the basin's spring-summer streamflow through icemelt , though the hydrological regime is typically dominated by contributions from snowmelt (Ayala et al., 2016;Burger et al., 2019;Shaw et al., 2020b). The catchment is characterised by steep terrain (average slope of 27°) that supports small sclerophyllous vegetation. Snow is typically present between early May and until late December at the highest elevations, with an estimated peak SWE in early to mid-September (Cornwell et al., 2016). Drought conditions since 2010 have resulted in a reduction in the number of snow cover days and increase in the snow line altitude (Mernild et al., 2016;Garreaud et al., 2017;Saavedra et al., 2018), with 2019 being the second driest year on record (as measured at Quinta Normal station, Santiago: 1950-present) up to the date of the 2019 Pléiades acquisition (Digital Elevation Model Generation).

Data and Methods
The following methodology for derivation of snow depth from Pléiades optical satellite imagery builds upon the initial work of Shaw et al. (2020a) who evaluated its use in the same catchment for a single year. This paper, however, aims at observing and understanding the key detectable snow depth differences between two end-of-winter scenes for the first time and characterising the expected uncertainties in an updated workflow that was not FIGURE 1 | A map of (A) Chile, (B) the metropolitan region and (C) the Rio Yeso catchment. The automatic weather station Quinta Normal, Yeso Embalse and Termas del Plomo are marked in panels b and c. The hillshade in b is provided by ASTER GDEM and the hillshade/digital elevation model scale in c is that of SnowOFF_2018. Glacier outlines were manually digitised from the Pléiades SnowOFF_2018 orthoimage.
FIGURE 2 | A flowchart of the steps taken to process Pléiades digital elevation models to assess snow depth difference (SDD) in this study. Steps (A-E) are referenced in the methodology text.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579142 3 addressed in previous studies. The flowchart in Figure 2 provides an overview of the processing steps with respect to the Pléiades datasets. The following subsections detail the specific processing steps with reference to this workflow.

Digital Elevation Model Generation
DEMs were generated from triplets of high resolution (0.5 m) imagery using Pléiades 1A and 1B satellites (Marti et al., 2016;Deschamps-Berger et al., 2020;Shaw et al., 2020a) (Figure 2A). Image triplets were taken during four separate acquisition dates ( Table 1) that were either snow-free (hereafter SnowOFF_y, where y refers to the year of interest) or the equivalent snowcovered terrain in late austral winter (hereafter SnowON_y). The SnowON acquisitions were obtained for the 4th and 2nd September in the years 2017 and 2019, respectively. The reference SnowOFF acquisition was obtained January 6, 2018 and was used to generate dDEM maps of snow depth by Shaw et al. (2020a). This scene had minimal snow cover at the highest elevations of catchment (<0.1% of pixels >5,100 m a.s.l. -Supplementary Figure S1) that were ignored in this analysis. A second SnowOFF image set was acquired on April 24, 2019. Unfortunately, this acquisition shortly followed the first small snowfall of the 2019 winter, and thus, the 2018 image triplets (SnowOFF_2018) were used to generate the reference snow-free DEM surface for both years. The SnowOFF_2019 DEM was, however, used to update surface elevation changes over glacierised areas of the catchment (see Correction for Glacier Change).
We applied the semi-global matching stereo algorithm with a binary census transform cost function following Deschamps-Berger et al. (2020) to generate point clouds from Front-Nadir (FN) and Front-Back (FB) image pairs in the stereo process of the Ames Stereo Pipeline [ASP - Shean et al. (2016)]. This function generates two disparity maps from the FN-FB pairs, which are combined in a single point cloud with a joint triangulation process. We opted for a FN-FB set as it offered a larger disparity in point clouds when calculated within the ASP routine, and therefore found as the best trade-off in reducing intersection errors and data gaps due to terrain occlusion. The semi-global matching, binary transform routine in ASP was found by Deschamps-Berger et al. (2020) to reduce the random error and root mean squared error of both snowcovered and stable, snow-free terrain compared to local search algorithms previously implemented for snow depth mapping using Pléiades (Marti et al., 2016;Shaw et al., 2020a). Due to saturation in the SnowON_2017 images (see discussion Limitations of Pléiades for Multi-Annual Snow Depth Derivation), correlation failure of the stereo process in ASP resulted in a loss of ∼17% of the total pixels in the catchment (Shaw et al., 2020a). No saturation was evident in the SnowON_2019 images.
We subsequently applied the point2dem function in ASP in order to map the stereo-generated point cloud onto a rough 50 m grid and provide an adequate absolute georeference for the DEM. This coarse DEM was therefore used as a projection base for the ASP routine mapproject. This routine produces roughly rectified images which were again processed with stereo (Shean et al., 2016). The gridding of the point cloud was this time made at the final resolution of 3 m. We selected this resolution as a balance of accuracy and data gaps in the final DEM product (Marti et al., 2016). This resolution was also considered appropriate by Deschamps-Berger et al. (2020) when comparing Pléiadesderived snow depths to ASO in the Sierra Nevada, United States.

Land Surface Classification
We utilised the ASP mapproject function to ortho-rectify the nadir multispectral image for each acquisition date based upon the respective DEM (Table 1; Figure 2B) and provided a manual training classification of the following surface classes; i) snow, ii) bare ground, iii) snow in shadow and iv) bare ground in shadow. Vegetation and water bodies are almost non-existent and not included in the classification scheme. Each class was digitised manually in ArcMap 10.3 to provide a training set for a random forest classifier in the Orfeo toolbox (Grizonnet et al., 2017). A coauthor of the study provided an independent digitisation of snowcovered/snow-free pixels that were used to evaluate the random forest classifier. Pixels from each class were sampled randomly across the entire space of the catchment. We derive an overall Accuracy Index as the sum of all the true positive scores of a confusion matrix divided by the total number of validation points (523 and 487 in 2017 and 2019, respectively). Using this metric, we found a score of 0.95 and 0.93 for the images of September 2017 and 2019, respectively, indicating a strong performance of the classifier.

Co-Registration of Digital Elevation Models
Using the SnowOFF_2018 DEM as the reference, we relatively coregistered all other DEMs (Table 1; Figure 2C) following the approach of Nuth and Kääb (2011). To horizontally translate and vertically bias correct each DEM, we extract only stable terrain (defined as snow-free, non-glacierized land) for slope angles <45° ( Berthier et al., 2007) using the classified images of each acquisition date (as described in Land Surface Classification). We applied the same translation to each equivalent orthoimage, to retain the same spatial coverage of the surface classification for further analyses. We used the classified bare ground (illuminated) of <45°to remove the median vertical difference of the SnowON DEMs relative to the SnowOFF_2018 reference. The values of the co-registration transforms are given in Supplementary Table S1. Following co-registration of each DEM, we differenced the SnowON_2017 and SnowON_2019 DEMs with the reference SnowOFF_2018 DEM to produce an initial raw difference map (hereafter referred to as a "dDEM") at a 3 m horizontal resolution. Hereafter the naming convention follows that described in Table 2, such that dDEMs refer to raw DEM differences and "SD/SDD" refer to DEM differences with filtered and/or gap-filled data that describe a "snow depth."

Correction for Glacier Change
Due to the relative dates of the SnowON_2017 and SnowOFF_2018 acquisitions, DEM differencing to obtain snow depth also included a signal of glacier mass balance and dynamics (glacier horizontal and vertical motion and subsequent snow and ice melt in the ablation zone). We apply the same corrections as in Shaw et al. (2020a) to account for these effects in the 2017 dDEM. The mean propagated uncertainty for these glaciers was calculated by Shaw et al. (2020a)

Topographic Indices of the Digital Elevation Model
To aid analyses of snow depth differences, we derived topographic indices for the catchment (Supplementary Figure  S2). These indices are similar to those as described in Shaw et al. (2020a), though updated based upon the SnowOFF_2018 DEM ( Table 1) that we derived from the new ASP stereo routine described in Digital Elevation Model Generation. These topographic indices are i) elevation extracted directly from the DEM (hereafter "ELE"), ii) slope angle (°, "SLP"), iii) the topographic position index that determines the relative elevation of a pixel to its surroundings ( (Revuelto et al., 2014) -"TPI"), iv) the catchment aspect (°, "APT"), v) the wind exposure parameter developed by Winstral et al. (2002) ("SX") and, vi) the calculated sky view factor ("SVF"). The SX parameter was calculated using the modal direction of hourly 10 m winds of the closest ERA5 cell for the period 1st April to 3rd September in each year (44.2°in 2017 and 47.4°in 2019). The period was selected as the beginning of the hydrological year (1st April) until the mean Pléiades acquisition date for the 2 years (2nd and 4th September). The ERA5 reanalysis wind vectors cannot represent localised wind conditions in the catchment, though it is a generalised representation of the regional wind direction. A TPI search distance of 60 m was applied following Shaw et al. (2020a).

Confidence Interval of Snow Depth Differences
We evaluated the quality of the inter-annual dDEM difference (d 2 DEM - Table 2) by assessing the empirical cumulative distributions of the stable (snow-free) terrain residuals ( Figure 2D). In doing so, we were able to attribute a range of uncertainty that is derived for 3 m pixels that should in reality be zero (i.e., no snow). We extended this uncertainty range to explore how topographic indices of the basin (Topographic Indices of the Digital Elevation Model) might affect the quality of the d 2 DEM, and thus the assessment of snow depth differences between years (see Snow Depth Comparison). We extracted stable terrain pixels in 10 equal bins of each topographic index and obtained the 25-75% confidence interval of differences in each bin to provide a "strict" assessment of snow depth changes (Snow Depth Comparison), excluding differences that were within the stable terrain confidence intervals. With this, we assumed that the stable pixel residuals are representative of the differences over snow-covered areas. Given that the stable terrain pixels were well distributed in space and across the full range of topographic 2 | Key acronyms to describe elevation models or elevation model differencing to calculate snow depth of inter-annual snow depth differences in this study.

SnowON/OFF_y
The individual digital elevation models generated from the ASP routine (Digital Elevation Model Generation), where "ON" ("OFF") refers to snow-covered ('snow-free) scenes for year "y" indices we believe that this is a fair assumption, though a full test of the relative uncertainties for stable and snow pixels would require a high resolution reference dataset or the entire catchment (e.g., Deschamps-Berger et al., 2020) that was not available to this study, nor would typically be available in studies of this kind.

Snow Depth Comparison
To assess dDEMs of individual years as "snow depth maps" [i.e., an analysis ready "end-product" -Shaw et al. (2020a), Shaw et al. (2020b)], we i) ignore any remaining snow-classified pixels that are negative in the dDEM of either year, ii) set stable terrain differences to zero, and iii) limit extreme outliers of each dDEM for high slope angles following the relation of slope holding capacity given by Bernhardt and Schultz (2010) (<3% of total pixels).
For evaluation of pixel-to-pixel snow depth differences (SDD), we considered: (1) where "SD" is the single year snow depth map (  Figure 2E). For a stricter comparison, we limited SDD by excluding pixels that are over glacierised areas, pixels that are within a buffer area of 6 m (2 pixels) around saturated areas of the SnowON_2017 image (Digital Elevation Model Generation) or pixels that were within the calculated confidence limits of the d 2 DEM over stable terrain for all topographic indices (Confidence Interval of Snow Depth Differences; Supplementary Table S2).
Alternatively, to provide an estimate of a full, catchment-wide SDD, we filled any remaining data gaps of a single year (such as those resulting from image saturation) using a random forest model trained on all available snow depth data and topographic indices, following Shaw et al. (2020a). We considered that mean differences at this full catchment scale is thus the weighted sum of the differences across all elements of analysis: where the proportional area of elements of A (the mean change is zero -no snow in both orthoimages), B (mean change is not significant -within the confidence ranges of topographic indices), C (elements of greater uncertainty -glacierised or gap-filled pixels) and D (the 'strict' SDD where mean change is outside the confidence range) add to 1. Finally, we also investigated the areas of snow absence in each year where the equivalent pixel is snow-covered for the other year (i.e., snow-covered in 2017 and absent in 2019), and we related these differences to the topographic indices above.

Uncertainties at the Catchment Scale
It is not straightforward to analytically estimate the error of the mean snow depth over a catchment as it results from the combination of different errors correlated at different spatial scales (Anderson, 2019). Another, typical approach would be to infer the error from the stable terrain residual. However the median error over stable terrain is zero by correction of the vertical bias (Co-Registration of Digital Elevation Models). We estimated the error on the mean snow depth of the catchment from previous studies which compared Pléiades snow depth with independent measurements. Shaw et al. (2020a) found a 0.22 m bias against terrestrial LiDAR data over a small region of the same catchment while Deschamps-Berger et al. (2020) found a bias ∼0.1 m over a similarly complex terrain of comparable size in the Sierra Nevada, United States. Therefore we considered that σ SD 0.2 is a realistic estimate of the mean error (with the exception of glacierised areas where σ SD is considered as 0.3 m -Correction for Glacier Change). We further assumed that the mean error on SD_2017 and SD_2019 are independent (due to different coregistration vectors, image radiometry and base to height ratios of the winter images - Table 1; Supplementary Table S1) so that the uncertainty of strict SDD can be written as: For the catchment-wide SDD, we follow the structure of Eq. 2 and weight the uncertainties by area: is equal to σ 2 SDD(D) as B and D areas differ only by the significance of the SDD signal. σ 2 SDD(A) was set to zero as we consider that the information about the lack of snow in element A is perfectly known. For C, the sigma value follows Shaw et al. (2020a) who calculated the uncertainty of snow depth corrections over glacier areas due to ice melt and glacier dynamics. The catchment-wide uncertainty can thus be simplified by:

Nivo-Meteorological Data
To relate the quantity and distribution of snow in each year to meteorological conditions during the preceding winter months, we utilised local automatic weather station (AWS) data for the region. We obtained data from the high elevation (2,995 m a.sl.) AWS, Termas del Plomo ( Figure 1C) for the period 1st April to 31st October in both years. The station is maintained by the Chilean Water Directorate and contains information on air temperature (°C), relative humidity (%), wind speed (m s −1 ) and direction (°) and precipitation (mm/h). During 2019, precipitation data are missing for a period of 3 months (May-July). Therefore, we supported this data series with precipitation records from the long-term AWS, Yeso Embalse (2,475 m a.s.l.) for the same period ( Figure 1B). We consider the beginning of the hydrological year (April) until the end of October in order to provide context to the nivo-meteorological conditions preceding and following the image acquisitions. To explore the time series of snow cover in the preceding winter of each year, we extracted the daily catchment fractional snow cover area based upon a 0.2 threshold of the normalised difference snow index (NDSI) from the MODIS (MODerate resolution Imaging Spectroradiometer) MOD10A1 V6. dataset. This dataset was not validated within our study, though provides a generalised view of the snow cover evolution preceding the Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579142 6 Pléiades acquisitions. Finally, to provide a historical context of climatic conditions for the 2 years under investigation, we also obtained air temperature and precipitation records from the long-term (1950-present) (Figure 4). The 2019 winter was both drier and warmer, though typically experienced similar patterns in the speed and direction of the wind, interpreted from both Termas del Plomo AWS and the ERA5 reanalysis ( Figure 4). 2019 experienced a greater quantity of winds from the northwest, though mean wind speeds were similar and no strong differences in the wind patterns were discernible between years. For the period of available data, there were no significant inter-annual differences in the measured incoming shortwave radiation for the lower catchment. Figure 5 shows the stable terrain residuals for individual years and the empirical cumulative distribution of stable terrain residuals between years (d 2 DEM). Inter-annual (2019-2017) stable terrain differences for the entire catchment reveal an uncertainty range of −0.36 to +0.24 m for the 25-75% confidence interval ( Figure 5C), suggesting that within such range there is a 50% chance that SDD are the result of error. For bins of each topographic index (Topographic Indices of the Digital Elevation Model), uncertainty ranges are highly variable (see Supplementary Table S2), with ranges of −1.5 to +0.25 m for slopes >45°and the most convex TPI bin. For shallow slopes (<25°) with high sky view fraction, the range of uncertainty is typically within ±0.2 m. Figure 6 shows the spatial snow depths derived from Pléiades DEM differencing in each year at a horizontal resolution of      Figure 7; Topographic Indices of the Digital Elevation Model). The grey shaded area denotes the range of uncertainty in the stable terrain (d2DEM) for each bin of each index. The right-hand axes represent the ratio of stable to snow pixels for each bin and index (reference of 0 and 1 given by horizontal dashed lines). Low ratios are noted to still have >20,000 stable terrain pixels in total. Arrows for TPI and SX are used to clarify the meaning of positive and negative values for those indices.

Pléiades-Derived Snow Depths
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579142 9 Classification), the calculated snow cover area is considered to be ∼60% in 2017, compared to 50% in 2019. We find mean snow depths of 1.27 ± 0.2 m and 0.92 ± 0.2 m in 2017 and 2019 respectively.
Mean snow depths are related similarly to the topographic indices in each year (Figure 7) whereby sheltered, east facing slopes of 10-30°possess the deepest snow depths. Those sites which represent moderately low sky-view fractions (0.5) are also subject to large inter-annual differences in snow depth, likely as a result of avalanche deposits at the base of steep slopes. In both years, the higher elevation (>4,500 m a.s.l.), steep slopes of a north facing orientation have the shallowest snow depths, as expected due to the effect of slope-dependent holding capacity, wind redistribution and radiative melting/sublimation of snow. The meteorological controls on snow melt/distribution are deemed to be similar based upon the observations from both AWS and reanalysis data sources. However, 2017 again shows a larger variability across all topographic indices.

Quantifying Inter-annual Snow Depth Differences
The mean and standard deviations of SDD are plotted against topographic indices in Figure 8 with the related d 2 DEM stable terrain confidence intervals (grey shaded areas). The smallest confidence intervals are related to flatter and open terrain (Evaluation of Pléiades d Digital Elevation Models; Figure 8B; Supplementary Table S2) and most SDD values lie outside these confidence ranges. The largest SDD values are found for southfacing, concave slopes of slope <30° (Figure 8). Figure 9 shows the pixel-to-pixel SDD for the catchment, again highlighting the mostly negative differences due to dry conditions of the 2019  Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579142 11 austral winter. Total mean catchment-wide SDD ( Figure 9A) is −0.34 ± 0.29 m, with a volume change of −38 ± 20% in 2019, relative to 2017. Figures 9B,C shows the "strict" SDD, based upon removal of glacier and saturation areas and pixel differences within the confidence limits related to topography (Figure 8). With 23% of the total catchment area remaining for this analysis, we find a mean snow depth difference of −0.48 ± 0.28 m, and a 39 ± 15% reduction in 2019 total snow volume compared to 2017. The mean of negative differences is −0.84 m and +0.33 m for positive differences (i.e., for deeper snow depths in 2019). A large proportion of these negative differences appear to be associated with gullies and the base of steep, north-easterly slopes, though positive differences are more sporadic with no clear visible pattern ( Figure 9C).
Negative SDDs are associated more often with shallower slopes ( Figure 10B) that are more sheltered ( Figure 10E) and are concentrated in a bimodal pattern of east and west facing slopes ( Figure 10D). Smaller SDD between years, those that are within the confidence limits of the Pléiades observations, are typically found at lower elevations with steeper, convex slopes than the larger negative differences (not shown), though otherwise follow similar patterns related to topography. Positive SDDs are more common for steep slopes of an east facing orientation, though do not relate well to any of the topographic indices.
We find that in 2019, snow absence is typically more common at lower elevations ( Figure 11A), related to the warmer conditions of the 2019 winter (Figures 3, 4). Absence of snow in 2019 is common for west facing pixels of the catchment ( Figure 11D) and those that are with a higher sky view factor ( Figure 11F), suggesting a possible link to the radiative energy balance of the catchment and the relative timing of winter snowfalls and cloud-free days during the preceding winter.
Absence of snow in 2017 shows a lesser, albeit similar relation to each topographic index, though is most noteworthy for slopes >30°( Figure 11B). The mean of 2017 and 2019 snow depths were 0.68 and 0.70 m, respectively for the areas of snow absence in the alternative year. These areas represent 9.9% and 1.7% of the total catchment area, respectively.

Snow Depth vs. Meteorology of the Central Andes
The recent, prolonged occurrence of drought for the central Chilean Andes is of serious concern for the future water security of the region (Meza et al., 2015;Garreaud et al., 2019). Ayala et al. (2020) estimated that glacier contributions to streamflow of the Maipo River (1955-2099) may have already peaked if climate was to stabilize at the level of the past two decades. Such model findings of the last decades also support the large scale observed patterns in glacier mass balance related to the so-called "mega-drought" (Braun et al., 2019;Dussaillant et al., 2019;Farias-Barahona et al., 2019;Farias-Barahona et al., 2020a;Farias-Barahona et al., 2020b;Masiokas et al., 2020). Nevertheless, the relative contributions of glaciers are governed by the hydrological snow year in the upper catchments of the central Andes, which are typically snow-dominated. For example, estimates of annual streamflow contribution from snowmelt for the Rio Yeso catchment are up to 97% for wet years  and estimated as 54 ± 10% on average for the entire Maipo catchment of 4,843 km 2 between 1955-2016 (Ayala et al., 2020).
The key challenge of all such modelling studies in mountainous regions to date, is the uncertainty in the distribution of solid precipitation in space and time (Ragettli et al., 2014;Ayala et al., 2016;Burger et al., 2019;Ayala et al., 2020;Shaw et al., 2020b), largely governed by the lack of high elevation data and/or the reliability of high elevation snow-precipitation records (e.g., Scaff et al., 2018). Shaw et al. (2020a) addressed this challenge by using Pléiades DEM differencing to obtain an estimate of snow depth for the Rio Yeso catchment. Their study showed that deviations from this observed snow volume could be up to ∼30% if modelling snow depths based upon local meteorological information from the catchment outlet (Termas del Plomo AWS presented here). While this AWS provides one of the highest elevation observations available for the region, calibrating a precipitation gradient from this and other, high elevation weather stations of the region was alone not sufficient to characterise the accumulation and likely re-distribution of snow observed over the >2,400 m elevation range of the catchment.
In this study, we have expanded this research focus to investigate the differences in snow depth between 2 years of the drought period (Figure 3), updating the general workflow based upon recent findings of Deschamps-Berger et al. (2020) and characterising the uncertainties in a robust manner. Our findings suggest that static, elevation-dependent precipitation gradients (such as in the aforementioned studies) would not be sufficiently robust to apply to multiple years of AWS data in the hopes of estimating high elevation snow volumes (Figure 7). For example, we find approximately an 80% reduction of measured 2019 winter precipitation at the Yeso Embalse AWS ( Figure 1B) when compared to the 2017 winter, though we find a difference of −39 ± 15% in the snow volumes of snow covered cells outside the 25-75% confidence limit (23% of the total catchment) and −10% difference in the snow-covered area of the entire catchment.
The consistent patterns of snow-covered area between years, despite sizeable differences in the observed volume of snow ( Figure 6; Supplementary Figure S2) implicates the strong control of both topography (such as exposure or slope) and its interaction with local meteorological conditions (such as incoming shortwave radiation and wind speed/direction) in determining the spatial distribution of snow. Whereas snow volumes have an evident connection to the total precipitation occurring in the basin or in nearby basins (e.g., Termas del Plomo or Yeso Embalse AWSs), the topography acts as a stronger predictor for snow melt, sublimation or mechanical removal (avalanching or wind redistribution), producing very similar inter-annual patterns of snow cover for the end-of-winter. The relationship of snow depth and meteorology can also be inferred from the catchment topography, which combines orographic precipitation differences with elevation up to a threshold (∼4,000 m a.s.l.) followed by the dominant role of ablation or redistribution at high elevation steep and north facing slopes. The mean snow depths at elevations >4,000 m a.s.l. are mostly sustained by concave landforms that Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579142 receive avalanched snow and protect it from both ablation or wind redistribution (TPI, SX and SVF indices). Comparing records of precipitation at the Embalse Yeso and Termas del Plomo AWS in 2017 (not shown), we find that the former could not represent all of the same precipitation events seen by Termas del Plomo AWS in the study catchment . Therefore, we can infer that at least some of the differences in the Yeso Embalse precipitation records and Pléiades observations may be explained by orographic precipitation differences between these stations or by storms originating from Argentina. Additionally, the high uncertainties of volume change at the full catchment scale (−38 ± 20%) and a greater snow compaction of thicker snow packs in 2017 may also partly affect the distinct differences between observed precipitation and snow depth for the catchment.
Unfortunately, Termas del Plomo AWS records were partly missing for the 2019 winter and so differences in timing and magnitude of precipitation events cannot be fully assessed here. Even so, this emphasises the immense difficulty in attempting to model changes of large catchment hydrology using statisticallyderived forcing of precipitation and highlights the potential advantages of satellite stereo DEMs for calibration, validation or initial conditions in glacio-hydrological models (Vögeli et al., 2016;Shaw et al., 2020b). While a greater number of highelevation AWSs would be highly recommendable to support such glacio-hydrological modelling and improve process understanding of snow accumulation and remobilisation of catchments in the Andes , the inferred radiative and wind processes based upon high resolution topography and low-elevation meteorological AWSs are perhaps sufficient to account for the smaller expected differences in snow-covered area. A key challenge remains for linking precipitation events, wind variability and snow depth at the highest reaches of the mountain catchments, especially because the available data for snow depth offers only a snap shot in time ( Figure 4). Nevertheless, the emergence of decimetre accuracy Pléiades and World-View products to estimate snow depth (Marti et al., 2016;McGrath et al., 2019;Shaw et al., 2020a;Deschamps-Berger et al., 2020) may offer new potential in how we can constrain high mountain precipitation, though further steps are required to relate localised in situ meteorological observations to catchment-wide snow volumes on an inter-annual basis.

Limitations of Pléiades for Multi-Annual Snow Depth Derivation
We have presented 2 years of snow depth estimates for almost the same day of year, using the same methodology provided by Deschamps-Berger et al. (2020). The fact that we difference both September DEMs with the same SnowOFF_2018 reference DEM arguably creates a more robust test of the inter-annual differences in snow depth. Nevertheless, the inconsistent acquisition geometries between years means that DEM quality is variable and resulted in data loss. For example, the issues of saturation in the 2017 image acquisition had limited the area of the basin that was comparable between years (by approximately 17% - Shaw et al. (2020a)). This was largely due to non-standardised acquisition parameters set by the tasking, an issue that can be mitigated by specifying the time domain integration (TDI) lines to image the scene. We therefore highlight that there are limitations that are specific to our study (such as saturation or image timing and change detection over glaciers) that are not necessarily representative of the capabilities of Pléiades for snow depth mapping in other domains (Marti et al., 2016). Nevertheless, future tasking of images for this region could benefit from standardisation of acquisition parameters that lead to fewer image saturation problems and data loss.
A further issue with the work presented here is that we are unable to quantify inter-annual variability of shallow snow depths (i.e., small DEM differences) that are particularly prevalent under recent drought conditions, especially for the winter of 2019. Approximately 9% of the total 3 m pixels from the Pléiades acquisition area were within the 25-75% confidence limits calculated from the topographic indices (grey band of Figure 8). This is therefore a sizeable proportion of the catchment that cannot be reliably considered in analyses for inter-annual differences. Nevertheless, evidence of snow-covered area changes from the classification of the high resolution e.g., multispectral Pléiades images allows an analysis of factors that may have affected the presence/absence of snow in a given year ( Figure 11). From this we establish that exposure, slope and aspect control a large proportion of the absences of snow in 2019 ( Figure 11) that can likely be explained by radiative fluxes during the winter. Given that differences in snow-covered areas are therefore typical for steeper slopes, we can assume that snow depths for the year where snow is still present are small given the likely slope holding capacity (Bernhardt and Schulz, 2010). Such smaller differences in spatial snow depth have been shown to have minimal impact upon models of monthly streamflow (e.g., Shaw et al., 2020b) and therefore uncertainty of the reported magnitudes may not be crucial, especially when spatially resampling to coarser resolutions (Deschamps-Berger et al., 2020).
In relation to the relative errors for stable and snow-covered terrain (Figure 8), we assumed that stable terrain uncertainties were representative of the SDD uncertainties, and this can be partially supported by the fact that the stable terrain pixels are well distributed across all of the topographic indices ( Figure 8).
As previously mentioned, studies of this type typically lack a large-scale reference dataset, such as that provided by ASO. Comparing Pléiades and ASO snow depths, Deschamps-Berger et al. (2020) in fact found a higher error for snow-covered pixels than over stable terrain. The authors partly attribute this to the presence of vegetation and to the co-registration process optimized on stable terrain within their study catchment. We highlight that at least the former is not applicable to the Rio Yeso catchment due to the lack of vegetation. Nevertheless, a full assessment of the stable vs. snow errors are not possible given the data availability of this study.
We estimated the mean error from studies which measured error on similar mountainous terrain with independent validation datasets (Deschamps-Berger et al., 2020;Shaw et al., 2020a). We are aware that this neglects the specificity of each study site, though given that the error was estimated over a catchment of similar size and complex topography by Deschamps-Berger et al. (2020) believe that it provides a typical range of errors that are relevant for the Rio Yeso catchment. Finally, we identify that our late winter DEMs aid snow depth estimation only for a brief snapshot in time (Margulis and Fang, 2019;Shaw et al., 2020b). Therefore, these single yearly observations may obscure partly the causes of snow depths in the months preceding and thus it is difficult to build a story of the full winter-time accumulation record (Figure 4). Hedrick et al. (2018) utilised multiple snow depth maps of a single winter season from the aforementioned ASO data in order to update the SWE conditions within an energy balance routine. Their study found a large improvement in SWE estimation if the model had previously been updated with observations of the ASO, thus demonstrating the advantages of a higher temporal resolution for such detailed spatial snow depth datasets. While this could be a recommendation for future work on snow depth derivation and has previously been performed using Pléiades data for glacier mass balance studies (Belart et al., 2017), multiple spaceborne stereo DEMs obtained throughout the season likely could not detect smaller snow accumulation events within the likely uncertainties reported (Marti et al., 2016;McGrath et al., 2019;Deschamps-Berger et al., 2020;Shaw et al., 2020a). Our study has provided a generalised workflow of comparing Pléiades-derived snow depths for two distinct years. However, the calculated levels of uncertainty are both difficult to fully quantify, but also clearly larger than LiDAR-derived snow depths from dedicated surveys (Painter et al., 2016). Accordingly, calculating snow depths and snow depth differences from spaceborne photogrammetry, as presented here, one should carefully plan acquisitions based upon sound knowledge of the local conditions (such as the typical length and intensity of the core winter season and timing of the SWE maxima).

Future Directions
Further years of data, even at single "well-timed" end-of-winter dates could have benefits for modelling seasonal streamflow (Li et al., 2019;Margulis and Fang, 2019;Shaw et al., 2020b), calibration of suitable parameterisations for snow accumulation (e.g., solid precipitation threshold) or elevation gradients of forcing variables, namely precipitation (Ragettli et al., 2014;Ayala et al., 2016;Vögeli et al., 2016;Burger et al., 2019;Ayala et al., 2020). In combination with the increasing resolution of (merged) snow cover products at short revisit periods (e.g., Gascoin et al., 2019), leveraging more, "welltimed" Pléiades imagery may improve our ability to characterize the snow cover distribution and evolution in mountain environments, in particular for data-scarce regions such as the central Andes. We recognise a potential limitation that our comparison years (and possibly future observation years) were during a prolonged and severe drought period (Garreaud et al., 2019), and therefore any parameterisations based upon meteorology or topography using these snow depth datasets may not be applicable to years not under such drought conditions. Nevertheless, the clear contrast of information from in situ observations and high resolution spatial snow depths presented here argue a strong need for further research into forming a link between the high spatial (Pléiades) and temporal (winter AWS records) data series with further years of high resolution data.

CONCLUSION
We present, for the first time, a comparison of 2 years of high resolution, 3 m Pléiades-derived snow depth maps for a high mountain catchment of the central Chilean Andes, building on an initial study for the same catchment. Applying a recently updated methodology for tri-stereo photogrammetry and robustly assessing uncertainties, we compare the pixel-to-pixel and grouped snow depth differences based upon topographic indices (e.g., slope, aspect, sky view) between the late winter of 2017 and 2019.Both years represented drought conditions for the region, though 2019 experienced the second driest winter for the observation record (1950-present), resulting in an 80% reduction of winter (April-September) precipitation compared to 2017 at a high elevation AWS (2,475 m a.s.l.) available in both years. Nevertheless, we find a 10% reduction in snow-covered area for the total catchment and a mean −0.48 ± 0.28 m snow depth difference and 39 ± 15% reduction in total snow volume (2019−2017) for the pixels restricted by stable terrain uncertainty ranges of binned topography indices (23% of the total catchment). We find that the largest differences in snow depth and snow-covered area are related to the aspect, exposure and slope, such that greater depth (and presence) of snow in 2017 was found for west-facing, sheltered pixels of moderate slope (30-40°). Areas of greater snow depth in 2019 appear related to north-easterly slopes of 30-40°, though no clear explanations relating to wind remobilisation or other meteorological phenomena are apparent. Smaller snow volumes of 2019 relate to warmer conditions, lack of winter precipitation and a likely role of the winter energy balance on thin and patchy snowpacks. Nevertheless, we identify a clear need to relate the high spatial (Pléiades) and temporal (AWS) scales within the available datasets in order to address the distinct differences between measured precipitation at lower elevations and observed snow depths at higher elevations. We recommend leveraging further observations of satellite-derived snow depths combined with newly developed snow-cover products to construct a history of the winter snow accumulation regime at such high mountain sites.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: http://doi.org/10.5281/zenodo. 3924571.

AUTHOR CONTRIBUTIONS
TS, CD-B and SG planned the experimental design. SG acquired the Pléiades imagery. CD-B and SG processed the DEMs with input from TS. TS wrote the manuscript with edits, inputs and revisions from CD-B, SG and JM.