Morphological indexes to describe snow-cover patterns in a high-alpine area

Abstract

The spatiotemporal distribution of snow affects hydrological and climatological processes at different scales. Accordingly, quantifying geometric features of snow-cover patterns is important, providing a valuable complement for snow water equivalent (SWE) modelling. This study on satellite-based morphological analysis originally uses two types of geometric indexes: (1) MN, Minkowski numbers (area (MN1), perimeter (MN2), Euler number (MN3)), and (2) CL, average chord length, to describe the morphology of Sentinel-2-derived snow-covered areas (SCAs), within the high-alpine site Zugspitze for a 5 year period. Results indicate that they capture the seasonal variability of snow-cover patterns, particularly during accumulation and ablation. Being to some degree independent from each other, MN2, MN3 and CL provide additional information upon shape, connectivity and length scale of snow cover, compared to most used indexes (e.g. fractional SCA). Correlation values up to +0.7 for MN2, +0.58 for MN3 and +0.46 for CL were observed with selected topographic characteristics, suggesting a close connection between geometric features of snow cover and ground features. Comparing in situ SWE measurements with MN and CL shows a correlation between −0.5 and +0.5. These indexes can hence be applied in combination with in situ data and/or modelling approaches to improve spatially distributed SWE in high-alpine catchments.

1. Introduction

Satellite remote sensing has been widely used for systematic monitoring of spatially distributed snow-covered areas (SCAs), even in near real-time mode, from global to catchment scales. Improved products on seasonal snow-cover mapping are being developing rapidly, and new satellite missions, such as Sentinel-2, are improving spatial and/or temporal resolution of the remotely acquired data. A large number of studies was carried out in the field of deriving SCAs, based on the normalized difference snow index (NDSI) (e.g. Hall and others, 2002; Härer and others, 2018; Gascoin and others, 2019), and fractional snow-covered areas (FSCAs), mainly based on spectral unmixing (e.g. Salomonson and Appel, 2006; Stillinger and others, 2023). Several studies found seasonal changes of (F)SCA to be correlated both with meteorological variables (e.g. Tang and others, 2017) and topography (e.g. Helbig and others, 2015; Poggio and Gimona, 2015; Helbig and others, 2021a, 2021b). The development of assessing snow water equivalent (SWE) using snow-hydrological models in conjunction with remote-sensing products has made significant progress (e.g. Rittge and others, 2016; Premier and others, 2022). Numerous studies have shown a strong correlation between topography and the distribution of snow depth (HS) and SWE during both accumulation and snowmelt periods. There are two main approaches being pursued in this regard. Firstly, researchers are utilizing modelling techniques of varying complexity to assess SWE (Durand and others, 2008; Brown, 2010; Brauchli and others, 2017, Diodato and others, 2022). Secondly, there is a growing interest in employing airborne or terrestrial LiDAR (Grünewald and others, 2013; Melvold and Skaugen, 2013; Helfricht and others, 2015; Skaugen and Melvold, 2019; Weber and others, 2020) as well as very high-resolution drone- and satellite-based stereo-photogrammetric observations (Mendoza and others, 2020; Shaw and others, 2020). Although several approaches to estimate SWE were pursued hitherto, this is still a very challenging task, and currently considered as a most important, not yet fully solved issue in the field of snow hydrology (Dozier and others, 2016).

In this context, it is crucial to improve the existing technologies and to explore novel approaches, e.g. in snow data assimilation methods, as well as to advance the development of accurate and reliable mathematical models that can estimate snow-related variables (e.g. SWE and HS). In this study, we pursue new ways to quantify snow-cover patterns based upon mathematical morphological indexes, i.e. (1) Minkowski numbers (MN) (Mecke and Klaus, 2000), describing the area (MN1), perimeter (MN2) and Euler characteristic (MN3) of a snow-cover pattern, and (2) chord length (CL) distribution (Roberts and Torquato, 1999). Hitherto, we are not aware of previous work using such morphological indexes in this field. We investigate how MN and CL change dynamically over space and time, as a result of snow accumulation and snowmelt processes, thereby mirroring variations in HS and SWE. While here we mainly focus upon the analysis of space–time variations of MN and CL and their correlations to topographic features, our aim is to exploit such tools to complement SWE estimation, as we will outline in more detail in our conclusions. Considering that the morphology and topology of seasonal snow-cover patterns in high-alpine areas are largely controlled by topography, and snow-hydrological and meteorological processes, the main objectives of the study are to (1) evaluate the effectiveness of MN and CL in quantitatively describing the geometric characteristics of snow-cover patterns seasonally, (2) quantify the influence of topographic features upon geometric properties of snow-cover patterns, such as shape, connectivity and scattering and (3) perform a first analysis of the correlation of in situ SWE measurements with MN and CL. In a further step, we hope to enhance the understanding of the processes shaping of the snow-cover patterns in high-alpine regions and to provide additional information for hydrological modelling thereby enhancing our ability to represent and simulate the spatiotemporal occurrence of snow in complex topography.

2. Methods

2.1 Study area

The study area is between Upper Bavaria (Germany) and North Tyrol (Austria), at the Zugspitze Massif, in the European Northern Limestone Alps (Fig. 1a). It has an extent of 29.6 km², a mean altitude of 2,078 m a.s.l. and fully covers in the upper parts the 11.4 km² research catchment Zugspitze. Annual average temperature therein is −4.5°C, and annual precipitation is 2,080 mm (Voigt and others, 2021; Weber and others, 2021). Altitude range spans over more than 2,100 m, encompassing Mount Zugspitze, the highest peak in Germany with 2962 m a.s.l., and the valley floors down to 830 m a.s.l. Large topographic variability and complexity lead to different meteorological conditions over the entire area, originating into distinct spatiotemporal snow-cover patterns during the snow season (Weber and others, 2020), making this site of particular interest snow wise. Most of the area is not or sparsely vegetated. Between ~1,400 and ~2,000 m a.s.l. the area is mainly covered with mountain pines, while below 1,400 m a.s.l. mainly coniferous forests are found. There are three very small glaciers in the upper part of the catchment (see Fig. 1a), i.e. (1) northern Schneeferner, (2) southern Schneeferner and (3) Höllentalferner, with a total surface of 0.34 km² in 2018 (Hagg, 2022). Nevertheless, in 2022 the southern Schneeferner lost its status as glacier due to a significant decline in both thickness and extent during 2018–21, and it is projected that the remaining ice may completely melt within 2024 (Mayer and others, 2022).

Figure 1. Overview of the study site Zugspitze: (a) spatial distribution of altitude and the locations of the glaciers, examples of snow-cover maps of five different days during the season 2019/20, representing typical phases of snow-cover evolution, acquired respectively on (b) 16/10/2019, (c) 18/11/2019, (d) 12/6/2020, (e) 6/8/2020 and (f) 18/9/2020, overlapped on a satellite image of the area under investigation (Google Earth Pro, 2021). The red lines in (b–f) show the outline of the snow cover/glaciers.

The high-alpine seasonal snow-cover dynamic at the study site displays five main stages, as seen from Sentinel-2-derived snow-cover maps (Figs 1b–f; see Section 2.2 regarding satellite data description). Figure 1b shows start of accumulation. The first snowfall events of the season start typically between the end of September and beginning of October in the upper parts of the catchment. Snow cover can appear patchy if some melt occurs in this phase. Figure 1c shows extended and homogeneous SCA. This happens typically in November/December–March in the upper and mid parts, and can reach the valley floors during the high-winter season. Figure 1d displays start of snow ablation. Snowmelt starts predominantly in the mid-elevation range usually during April and May. Snow is mainly patchy in the middle parts. Figure 1e displays the end of snow ablation season. Intense snowmelt with patchy snow in the upper parts occurs in June and July. Figure 1f show snow-free period. August and September are mainly snow free, with occasional snowfall events, with snowpack usually not reaching until autumn. The three glaciated areas are visible with bare ice or firn.

2.2 Satellite-derived snow-cover maps

For the snow-cover mapping exercise in this study, we used Copernicus Sentinel-2 images from the European Space Agency (ESA), available as level 2A product (atmospherically corrected) in the area of investigation, since November 2016. The multispectral images with visible, and near to shortwave infrared bands, were downloaded through the EO Browser of Sentinel Hub (EO Browser, 2022), supported by the ESA. We considered a 5 year time period during November 2016–October 2021, and we selected 129 images (excluding images with cloud cover), distributed in time as shown in Figure 2. The images originally come in a spatial resolution of 10 m in the visible bands, and 20 m in the shortwave infrared bands. However, the Sentinel Hub performs a downscaling procedure as due to WGS84 projection, achieving a value of 6.78 m for the 16-bit TIFF image downloaded within the area of investigation. The identification of SCA was carried out using the NDSI, with a threshold of NDSI ⩾0.42, based upon the Scene Classification algorithm developed by ESA, specifically for Sentinel-2 images (ESA, 2023). Notice that in this study shadows were blind spots. However, more sophisticated algorithms for snow-cover detection in satellite images, e.g. based on machine learning (Barella and others, 2022) could help to extent the investigated area.

Figure 2. Sentinel 2 - L2 products monthly availability (excluding images with cloud cover) in the study area for the period between November 2016 and October 2021.

2.3 Snow-cover pattern descriptors: Minkowski numbers and average CL

To quantitatively describe snow-cover patterns, we used two indexes, i.e. (1) MN and (2) CL, which have been recently used to describe morphological features of different binary patterns.

MNs are morphological descriptors related to curvature integrals, which characterize extent, shape and connectivity of spatial patterns (Mecke and Klaus, 2000). This set of morphological descriptors was successfully used to quantify complex binary objects, as well as their temporal dynamics (e.g. Schlüter and Vogel, 2011), and it has been applied in a variety of fields, in particular in soil sciences (e.g. Vogel and others, 2005), but also for snow microstructure analysis (Schleef and others, 2014). For a two-dimensional (2-D) space it is possible to compute three MNs (Parker and others, 2013), i.e. (1) total area of the pattern (MN1), i.e. the total extent of the snow-covered parts, (2) total perimeter of the pattern (MN2), i.e. the total length of the boundary between SCAs and no-snow areas and (3) Euler number (MN3). MN3 is given by the difference between the convex elements and the concave elements of the pattern, giving information on the connectivity between the elements of the pattern (thus the degree of homogeneity and clustering). In plain words, if a pattern is mainly composed of isolated convex objects, MN3 is a positive value, and describes the number of elements in the pattern. If a pattern is a connected network MN3 is a negative value, and the absolute value describes the number of meshes in the network (Vogel and Babel, 2006). MN3 is given as

(1)$${\rm MN}3( X ) = \displaystyle{1 \over {2\pi }}\mathop \smallint \limits_X^\;\left[{\displaystyle{1 \over R}} \right]dc\;\;$$

where X is the pattern, R is the radius of curvature of the circumference of a single object in the pattern and dc is an infinitesimal element of the circumference of the object (Vogel and others, 2005). To compute the MNs for the snow-cover patterns, we used the QuantImPy Python package (Boelens and Tchelepi, 2021). Intuitive examples of the MNs computed on synthetic binary patterns are shown in Figure 3a.

Figure 3. (a) Examples of MNs of different synthetic patterns with white circular elements of different radius, r, belonging to an area of 100 cm². (b) Schematic of CL measurements for cross sections of two-components random patterns. The chords are defined by the intersection of lines with the two-components interface. In the examples, two lines of segments in vertical and horizontal directions are represented.

CL distribution p(z) is a probability function associated with one of the two components of a binary pattern. It has a direct relationship with the morphology of a pattern, and therefore it is a useful tool to describe complex binary structures (Roberts and Torquato, 1999). The CL distribution describes the probability p_i(z)dz that a randomly selected chord identified in the component i (defined as a segment starting and ending at the boundary between the two components, see Fig. 3b) is characterized by a length in the range [z, z + dz] (Roberts and Torquato, 1999). All segments spaced by a given distance and orientated in the same direction are selected; different distance and orientations can be considered, in this case we decided for a distance of one pixel and vertical orientation (i.e. parallel to the north–south axis). The average CL is the first moment of the distribution, which is used in the study as a representative index, since it provides an estimation of the length scale associated with the elements of the component i (Schlüter and Vogel, 2011). The calculation of CL was carried out through the PoreSpy Python package (Gostick and others, 2019), originally designed to extract information from 2-D and three-dimensional (3-D) images of a porous material.

2.4 Topographic descriptors

The topography of the area was calculated from a 1 m DEM (LDBV, 2022) resampled to 6.78 m (resolution of snow-cover maps). We considered topographic descriptors of altitude, aspect, slope and curvature. Altitude is a main driver for snow accumulation and ablation: because of the temperature gradient, at the highest altitudes snow-cover accumulation starts earlier in the season, and lasts longer, often leading to a deeper snowpack. Aspect becomes especially significant during snow ablation, and needs to be considered due to both different effects of solar radiation, i.e. north-facing slopes are more shaded than south-facing slopes, and different effects of wind, i.e. leeward slopes are more subjected to snow accumulation then windward slopes (Lòpez-Moreno and others, 2013). The area of interest includes cells with aspect in almost all directions. We considered two components of aspect, i.e. (1) along the north–south direction and (2) along the east–west direction. This means their values increase monotonically with exposition to solar radiation. Accordingly, such components of aspect are more clearly related to snow processes: higher values indicate lower exposure to solar radiation, and often deeper and more extended snowpack, while lower values indicate vice versa. Slope and curvature mainly affect snow deposition due to gravitational transport. Slope is the inclination of the ground to the horizontal, and it is related to the local effect of gravity (Rowbotham and others, 2005), affecting processes such as deposition, and lateral snow transport (Bernhardt and Schulz, 2010; Bombelli and others, 2021); steeper slopes have generally thinner snow cover. Curvature is the convexity/concavity of the ground, with implications in the acceleration/deceleration, and convergence/divergence of materials like snow, affecting erosion and snow deposition (Hanzer and others, 2016). Convex features were attributed to positive values up to +1 and concave features to negative values down to −1. Histograms with the frequency distributions of each topographic feature are shown in Figure 4.

Figure 4. Frequency distribution of the topographic features in the study area: (a) altitude, (b) aspect (classes are defined as follows: north: −22.5° to 22.5°, north–east: 22.5–67.5°, east: 67.5–112.5°, south–east: 112.5–157.5°, south: 157.5–202.5°, south–west: 202.5–247.5°, west: 247.5–292.5°, north–West: 292.5–337.5°), (c) slope and (d) curvature.

2.5 Method for spatial correlation of morphology of snow-cover pattern and topographic descriptors

To assess the correlation level between snow-cover features as quantified with MNs and CL (Section 2.3) and the chosen topographic characteristics (Section 2.4), we calculated the spatial distribution of the morphological indexes and the topographical features within the entire area of investigation. To do so, we developed an algorithm to compute MNs and CL across a moving window upon the snow-cover maps. The moving window needs to be (1) large enough to cover a snow-cover pattern, representative of local geometry, and topology of snow distribution, and still (2) small enough to be represented by spatially averaged topographic indexes. To identify (a range for an) optimal window extent, a qualitative sensitivity analysis was performed. We computed the MNs in multiple areas within the catchment, using different side lengths ranging between 50 m and 1 km, during October 2020–September 2021. We then visually compared the resulting time series in order to identify the smallest window size that could still capture a representative snow pattern. Based on this analysis, we established that, for the purpose of this study, a window size of 250 m × 250 m was the minimum extent required to identify most significant temporal variations in the spatial distribution of the snow pattern morphologic indexes. This size was effective in capturing changes, both during the onset of winter accumulation and during the melting phase in spring/summer, while also effectively describing the local topographic features. Smaller window extents were found to capture poorly the spatiotemporal variability of snow-cover patterns. On the contrary, larger window sizes, although capturing well the variability of the indexes in time, were less representative of the underlying topography. Thus, using larger windows result in lower (estimated) correlation values between the indexes and topography. We believe that a window size of 250 m × 250 m is an optimal choice for our study, since it balances capturing significant temporal variations in snow-cover patterns, while adequately describing local topographic features. It is important to highlight that the optimal window size is influenced by the topographic characteristics of individual mountain regions. In areas with higher topographic variability, a smaller window size is required to represent the ground sufficiently, whereas in regions with lower topographic variability, a larger window size is beneficial to encompass less frequent changes of snow-cover patterns. In this regard, for future studies, the choice of the moving window size could be optimized using statistical algorithms, which take into consideration the prevailing topography and the magnitude of the geometric features of the snow-cover patterns. The outputs of the moving window approach are four raster files, displaying the spatial distribution of the four snow-cover pattern descriptors, MN1–MN3 and CL, for each of the snow-cover images. To link snow-cover patterns to topography, the same moving window frame was applied to compute a spatial average of the topographic descriptors. Examples of the spatial distribution of the morphological indexes, together with the distribution of the spatially averaged topographic descriptors are shown in Figure 5. Spatial correlation was quantified through Spearman rank correlation coefficient (ρ). Spearman's ρ is a non-parametric measure, and it quantifies strength and direction of a (monotonic) relationship between two variables, without an a priori functional assumption (e.g. a linear relationship).

Figure 5. (a) SCA of 21/09/2017 (top) and 16/10/2017 (bottom) and the corresponding spatial distribution of the morphological indexes: SCA (MN1), perimeter of the snow cover (MN2), Euler characteristic of the snow cover (MN3) and average chord length of the snow cover (CL). White areas regarding MN1, MN2, MN3 and CL represent no calculated morphological value either due to no or full snow coverage in the moving window frame. (b) Spatial distribution of the averaged topographic descriptors in the area of study. A moving window algorithm with a window extent of 250 m × 250 m was applied to compute the morphological indexes (a) and to average the topographic descriptors (b).

3. Results and discussion

3.1 Seasonal evolution of the morphological descriptors of the snow-cover patterns

Time series of MN and CL of the entire image extent of the Zugspitze region for the 5 year period November 2016–October 2021 show a marked annual periodic behaviour (Fig. 6).

Figure 6. Time series of MN and CL and the corresponding correlation values with examples of topographic features: (a) MN1 and correlation with altitude, (b) MN2 and correlation with slope, (c) MN3 and correlation with curvature, (d) CL and correlation with N/S and E/W components of aspect. Data have been computed on 129 snow-cover maps produced from Sentinel-2 satellite images of the Zugspitze research catchment acquired in a time period between November 2016 and October 2021.

3.1.1 Minkowski numbers (MN1–3)

Peaks of MN1 were detected mainly between November and March, and minimum values were achieved between August and September, describing the transition from (almost) fully SCAs to (almost) snow-free conditions (Fig. 6a). The increase from minima to maxima is faster than the decline from maxima to minima, in line with longer snow ablation at the end of the snow season, than SCA expansion at the beginning. MN2, i.e. perimeter, shows mostly low values when high values of MN1 occur, indicating an extended and homogeneous snow-cover distribution, but even when low values of MN1 occur, indicating a small amount of snow in the catchment (Fig. 6b). This means MN2 shows lower values, when the area is mainly homogeneously snow-covered or almost snow-free, and high values during the ablation period, and occasionally also in summer or autumn, with first snow accumulation at the highest altitudes, when swift changes in space of the snow cover occurs. Comparing the curves of MN1 and MN2, MN2 has larger interannual fluctuations. MN3, i.e. Euler number, shows a periodicity almost symmetrical to MN1 (Fig. 6c). When snow is quite homogenously spread over the entire area (mainly during the accumulation period), the patterns tend to be more connected, showing negative MN3 values. However, small fluctuations during the winter period are more prominent in MN3 compared to MN1. In general, when snow cover decreases, the patterns tend to be sparser and disconnected, showing positive MN3 values. High positive values resist until autumn, and as soon as most of the area is covered by snow, MN3 suddenly changes to negative values. MN3 also allows us to identify summer or early autumn occasional snowfall events. The visible change between a minimum and a peak in September 2017, for example, shows an early autumn extended, yet short-lasting (the area is almost snow-free again until end of October) snow accumulation. This event is also visible in a peak in MN2 and a minimum in MN1.

3.1.2 Average CL

CL shows a steep increase during accumulation, and a slower dynamic during ablation (Fig. 6d). Rather short peaks are specifically observed at times, when the conditions move from less/no-snow conditions to large accumulation. This is mainly the case at the beginning of the winter period, but it can also occur in spring, with a cold period, as it was for example the case in March/April 2021. Lowest values are observed mainly between July and September, when SCA achieves the least values.

Overall, MN and CL seem to quantitatively describe the temporal variability of snow-cover patterns well, and their seasonal and annual periodicity. Although MN1 alone gives already a good depiction of the spatiotemporal distribution of snow, MN2–3 and CL represent different features, and identify complementary morphologic information, so it is reasonable to consider them jointly for snow-cover pattern description. In fact, these indexes not only give information about the extent of the SCA, but they also quantify its shape, connectivity and length scale of snow cover. This is clearly visible in Figure 6 where MN1, representing the total extent of SCA achieves a plateau, usually between December and April (Fig. 6a), while the other indexes continue to exhibit more variability (Figs 6b–d). This shows that MN2, MN3 and CL help in capturing additional information of the geometric features of the snow-cover patterns.

To investigate how independent the information provided by each index is, we computed mutual linear correlation coefficients (Table 1). The results show that all indexes display a certain level of independence, and therefore each one carries additional information content, which might have some link to the controlling topographic features. Especially, MN2 is the least correlated compared to the others. For example, two patterns with a similar MN1 could be characterized by a totally different MN2, i.e. if one of them is more homogeneous and less dispersed, while the other one is more scattered. Similarly, two patterns with a similar MN2 could be characterized by very different MN3 values, if one of them consists of multiple separated small elements, while the other consists of few, but bigger branched elements. Compared to MN, CL seems to be less sensitive in the characterization of the ablation phase. However, it shows distinct peaks when the most parts of the site moves from no/less snow to a larger homogeneous snow extent, which is the case mainly in the early winter season, and this is therefore an added value further to the use of MN.

Table 1. Linear correlation coefficient between each index, computed on the entire images’ extent

3.2 Spatial correlation of morphologic snow-cover pattern and topographic descriptors

A spatial correlation analysis between the morphological indexes MN and CL and topographic descriptors (see Section 2.5) was performed, with the aim of getting an insight on how topography affects snow-cover processes described by morphologic indexes. Figure 6 shows the correlation between MN (a–c) and CL (d) with selected topographic descriptors. Correlation between the morphologic snow pattern descriptors and all topographic features is shown in Figure 8 in the Appendix.

3.2.1 Area (MN1)

Correlations between MN1 and all topographic features display the highest values, and the most marked periodic behaviour, when compared to the other indexes. The results confirm the tendency of snow to be present for a longer period in higher and flatter areas, as well as for northern and eastern expositions. In Figure 6a the correlation of altitude and MN1 shows values up to 0.8, for less/no-snow conditions, but lower values, close to 0, for high snow coverage. Interestingly, the correlation with altitude shows even slightly negative values, down to −0.25, during the high-winter season, possibly caused by extended/homogeneous snow cover at low and mid altitudes (flatter and more wind shielded ground), and more heterogeneous/fragmented snow cover at high altitudes (steeper and more wind exposed ground) (see Fig. 1d).

3.2.2 Perimeter (MN2)

Correlation between MN2 and topographic features shows a marked periodic behaviour, characterized by a wide range of both positive and negative values (see Fig. 6b). This suggests opposite effects of the same topographic feature on MN2 as varying in time, possibly caused by the intrinsically non-monotonic behaviour of this index. MN2 shows (1) negative correlation with altitude and aspect (both facilitating snow deposition/accumulation), achieving values down to −0.23 and −0.36, respectively, and positive correlation with slope and curvature (facilitating lateral snow transport), achieving values up to +0.60 and +0.40, respectively, from December to March (when there are extended/homogeneous snow-cover patterns). MN2 also shows (2) positive correlation with altitude, up to +0.71 and aspect, up to +0.68, and negative correlation with slope, down to −0.76 and curvature, down to −0.31, during June–September (when there are limited/sparse SCAs). This likely happens because within an extended SCAs, the biggest patches with no snow (high MN2) are placed at the lowest altitude and aspect, and where highest values of slope and curvature (convex features) occur, while at the highest altitude and aspect, and at lowest values of slope and curvature (concave features), snow tends to accumulate and/or last longer, and therefore areas without snow cover are smaller (low values MN2). On the contrary, in a sparse SCA, the biggest patches of snow (still high MN2) are placed at the highest altitude and aspect, and at the lowest slopes and curvatures (concave features), whereas at low values of altitude and aspect and high values of slope and curvature (convex features) snow is absent, or distributed in small and sparse patches (low MN2). A visual explanation of how MN2 is affected by slope and by aspect is given in Figures 7a and b (top panel), respectively.

Figure 7. (a) Comparison of spatial distribution of slope values in the area of study (top left) and SCAs (in red, top right, bottom left and bottom right) during different times of the season: first snowfall events and accumulation (10/11/2018), extended snow cover (13/2/2019) and melting phase (23/7/2019). The snow-cover maps explain why snow patterns show highest MN2 values at high values of slope (light-shaded areas) when there is an extended snow cover (13/2/2019), and highest MN2 values at low values of slope (dark-shaded areas) when there is a limited snow cover (10/11/2019, 23/7/2019). (b) Top: Schematic of why snow patterns show highest MN2 values at low values of aspect (south/west-facing slopes) when there is an extended snow cover, and highest MN2 values at high values of aspect (north/east-facing slopes) when there is a limited snow cover. Bottom: Schematic of why snow patterns show negative values of MN3 (prevalently concave features) on convex ground when there is an extended snow cover, and positive values of MN3 (prevalently convex features) on concave ground when there is a limited snow cover.

3.2.3 Euler characteristic (MN3)

Correlations between MN3 and topographic features show lower values, more frequent changes and a less marked periodic behaviour compared to the other indexes, which makes the interpretation more difficult. Negative correlation is observed against topographic features facilitating snow deposition, accumulation and later ablation, as they lead to more connected snow-cover patterns (negative MN3). Further, positive correlation is observed with topographic variables inhibiting snow accumulation and facilitating faster ablation, because they lead to more fractionated snow-cover patterns (positive MN3). As far as curvature is concerned, correlation shows positive values (1) during winter (extended and homogeneous snow cover), achieving values up to +0.34, suggesting a predominant effect of convex ground features in snow accumulation and redistribution processes, which leads to a well-connected pattern, and (2) during summer (limited and sparse snow cover) achieving values up to +0.24, suggesting a predominant effect of concave ground features, leading to a more poorly connected pattern (Fig. 6c). Both combinations of convex features/connected snow-cover patterns, and concave features/disconnected snow-cover patterns are characterized by a positive correlation. A schematic of the effects of curvature on MN3 in different seasons is shown in Figure 7b (bottom panel).

3.2.4 Average CL

Correlation between CL and aspect also does not show a clear periodic behaviour (Fig. 6d), however some observations can still be done, i.e. (1) the north–south component appears to show lower values of correlation and a less marked periodicity than the east–west component, suggesting that, in the area of investigation, the north–south component is less relevant to snow-related processes (Fig. 6e). (2) The east–west component appears to show positive correlation values up to +0.38 starting from the ablation period in spring, up to the first snowfall events in late autumn, and correlation values close to zero during winter (Fig. 6d). As CL is an average value, its use may entail some loss of information, especially for snow patterns with wide CL distribution, this is why the periodicity of its dependence against topography is less marked, when compared to MN.

In summary, the index showing the highest correlation values with topography is MN1, the second is MN2, MN3 third, followed by CL. It was possible to identify periods of the year where the effect of topography on snow process of accumulation, transport and melting is at least partly captured by the correlation with the indexes. All of them show the highest correlation values during the initial seasonal accumulation of snow (approximately September–December) and during ablation (approximately April–July), excluding correlation between MN3 and curvature. This is because (1) during the high-winter period, snow cover is mostly homogeneous and reaches mid- and low areas, and therefore in this period topography does not strongly affect the geometric patterns of the SCA given the high values of HS that dampen its effect; and (2) during the ablation phase, the snow cover gets more and more patchy, and in this period topography has a large impact on when and where the snowpack melts. The same phenomenon is visible during the initial season of accumulation, when topography controls the areas where snow is transported and accumulated. Further on, (3) during summertime, the snow-cover extent is zero, or very little, including the outlines of the three glaciers. Clearly, topography has only a minor impact in this period time, except for some occasional summer snowfalls. Finally, we believe it is important to note the large size of the data samples, which exceed 640,000 elements (number of moving windows used for calculation). This means that the obtained correlation values are well estimated and characterized by high-statistical significance.

3.3 Temporal correlation of morphologic snow-cover pattern and SWE

Since MN and CL showed a certain degree of correlation against chosen topographic features, and because topographic features are known to control snow processes such as accumulation and transport, thereby affecting SWE distribution, we performed a correlation test between the morphological indexes and corresponding values of in situ SWE measurements. To do so, we used SWE observations provided by an automatic weather station operated by the Bavarian Avalanche Warning Service (LWD) which is situated within the site at 2,420 m a.s.l. For this comparison, we focused on snow-cover patterns within a 500 m × 500 m area centred on the LWD station. The reason for selecting a larger area compared to previous analyses was due to the relatively homogeneous topography in this particular part of the catchment. Opting for a smaller area, such as 250 m × 250 m, would have led to snow-cover patterns dominated by either full snow coverage or no-snow cover at all. Consequently, the morphological indexes would not have identified interesting geometric features of snow-cover patterns and would have shown limited correlation with topographic features. By considering a larger extent of snow cover, we aimed to adequately represent the seasonal variations of snow cover and to ensure meaningful correlations. To measure temporal correlation, we employed the Spearman rank correlation coefficient between the morphological indexes MN and CL, and daily average values of SWE. Additionally, we examined the correlation against changes in SWE along 2 (previous) weeks, i.e. increase/decrease (observed during accumulation and melt phase, respectively), with respect to the date of each available snow-cover map. The results are given in Table 2. Therein, it is shown that some of the indexes actually display some degree of correlation/information, with/about the value of SWE, and/or its change in time. In general, the temporal correlations of HS and MN and CL would be similar to the ones performed with SWE. However, as snow density changes over season, e.g. the snowpack gets more compacted due to wind and radiation influences as well as rain-on-snow events, the connection of SWE and the morphological indexes might be represented differently in the accumulation and ablation phases, which should be considered in future research. This facet is to be further explored, but we preliminary observe that the use of MN and CL, and/or of their dynamics in time, may provide a complement in SWE estimation models, when information on the ground is not available.

Table 2. Temporal correlation between morphological indexes MN and CL and SWE, corresponding to an area of 500 m × 500 m in the research catchment

Values refer to Spearman rank correlation coefficient.

4. Conclusions and outlook

In this study, we used satellite images covering a 5 year period at the high-alpine Zugspitze region to verify that some morphologic indexes, namely MN1 (area), MN2 (perimeter), MN3 (Euler characteristic) and CL (chord length distribution) can (1) quantitatively describe different facets of seasonal snow-cover pattern evolution and their annual periodicity, and (2) show seasonally an overall medium-to-high correlation against selected topographic features. Most noteworthy results are found during snow accumulation in winter, and during snow ablation in spring/early summer. Therein, HS is limited, and the effect of topography on snow distribution is evident, leading to heterogeneous and unevenly distributed snow-cover patterns. During deep winter, fewer changes of the morphologic indexes and of their correlations against topography occur, because snow homogeneously covers the area, and likely deep snowpack masks the effect of topography on snow distribution processes. A similar effect regarding MN and CL can be observed in summer, given that almost no snow is present.

This study is meant to be a preliminary investigation of the potential of selected morphologic indexes to depict snow-cover patterns and needs further refinement, and investigation in a next step. Besides their value as descriptors of snow processes, we believe that the presented morphological indexes could improve snow-hydrological monitoring in alpine complex terrain. MN2, MN3 and CL indeed provide additional information as compared to MN1 or FSCA, which describe both the areas. This means, MN and CL indexes not only give information about the extent of SCAs, but they also quantify its shape, connectivity and length scale of snow cover. More information is given for example when achieves a plateau, which is the case usually between December and April, while the other indexes continue to exhibit more variability. Such information might be a valuable complement for SWE estimation, in addition to other commonly used indexes (e.g. FSCA). Also, the indexes demonstrate independence from each other, implying that they contain distinct pieces of information about snow-cover patterns. All of the indexes further exhibit a certain degree of correlation with some selected topographic features, i.e. they allow us to describe the geometric characteristics of snow patterns influenced by the underlying topography. Moreover, some degree of correlation was found against SWE, which proves that there is a connection between the geometric properties of snow-cover patterns, and stored water as snow on the ground.

Potential open paths for further investigation can be identified as follows.

1. (1) Given that one limitation of this study is the fact that the temporal resolution of the chosen cloud-free Sentinel-2-derived snow-cover maps might be too coarse to detect all relevant changes of the snow-cover patterns, particularly at the beginning of winter and during the melting phase, this should be improved. A logical next step would be to examine how the results would be affected by looking at higher temporal resolution, using e.g. MODIS satellite data or webcam images, the latter even being available to resolve a sub-daily resolution. In addition, by using higher (webcam) or lower (MODIS) spatial resolution sensors, this approach could be tested on different spatial scales, to verify whether other resolutions show meaningful results, especially compared to topography.

2. (2) Given that this study showed that it is possible to quantitatively describe snow-cover patterns jointly with the selected morphological indexes, we believe that this approach could be applied in combination with other in situ data and/or modelling approaches, e.g. to improve spatially distributed SWE and/or HS estimations on a finer spatial scale, which are still characterized by high uncertainties. We expect that complementing snow-hydrological models and monitoring setups with information on morphological indexes, e.g. through statistical modelling (Grünewald and others, 2013) and/or artificial intelligence framework, could improve HS and SWE derivation in the future. Some examples of application could be (i) to incorporate MN2, MN3 and CL as independent variables in a linear regression model for real-time SWE estimation (such as the ones applied in Yang and others, 2022), thereby combining satellite-derived information on snow covered area extent (FSCA) and morphology (MN and CL) with in situ snow pillow SWE measurements, physiographic data and historical SWE patterns, (ii) to use combined information on SCA, MN and CL, together with information on local topography, as constraints for 2-D or 3-D interpolation of ground observations of SWE, obtained from a network of snow pillows (Fassnacht and others, 2003; Molotch and others, 2004) and (iii) to exploit information on timing and geometry of snow-cover wasting/disappearance provided by MN and CL for a more accurate backward calculation of melting rates, to trace accumulated SWE back to the last significant snowfall (Molotch and others, 2008). In a next step, we aim to investigate whether adding MN and CL to information on in situ snow measurements, and snowpack modelling may improve spatially distributed SWE estimations by using different correlation analysis techniques (e.g. multivariate regression models, principal component analysis).