Using Artificial Neural Networks to Couple Satellite C-Band Synthetic Aperture Radar Interferometry and Alpine3D Numerical Model for the Estimation of Snow Cover Extent, Height, and Density

This work presents a new approach for the estimation of snow extent, height, and density in complex orography regions, which combines differential interferometric synthetic-aperture-radar (DInSAR) data and snowpack numerical model data through artificial neural networks (ANNs). The estimation method, subdivided into classification and estimation, is based on two ANNs trained by a DInSAR response model coupled with Alpine3D snow cover numerical model outputs. Auxiliary satellite training data from satellite visible-infrared MODIS imager as well as digital elevation and land cover models are used to discriminate wet and dry snow areas. For snow cover classification the ANN-based estimation methodology is combined with fuzzy-logic and compared with a consolidated decision threshold approach using C-band SAR backscattering information. For snow height (SH) and density estimation, the proposed methodology is compared with an analytical inverse method and two model-based statistical techniques (linear regression and maximum likelihood). The validation is carried out in Central Apennines, a mountainous area in Italy with an extension of about 104 km2 and peaks up to 2912 m, using in situ data collected between December 2018 and February 2019. Results show that the ANN-based technique has a snow cover area classification accuracy of more than 80% when compared MODIS maps. Estimation bias and root mean square error are equal to about 0.5 cm and 20 cm for SH and to 5 kg/m3 and 80 kg/m3 for snow density. As expected, worse results are associated with low DInSAR coherence between two repeat passes and snow melting periods.


I. INTRODUCTION
S EASONAL snow cover is the largest cryospheric component, covering more than 50 million square kilometers of the Earth surface (more than 31% of its land area) every year [1]. Snow cover area (SCA), snow height (SH), and snow density (SD) as well as snow water equivalent (SWE) are the main parameters characterizing the snow accumulation in mountainous regions. Snow cover patterns are governed by the effects of topography, land cover, wind redistribution, solar irradiance, and air temperature. Snow mapping is particularly important in meteorology, hydrology, and climate monitoring applications [2]. The considerable geographical extension of snow cover and its spatial heterogeneity makes it impractical to monitor the above parameters regularly (i.e., with a high spatial and temporal resolution) by means of direct or indirect in situ measurements. Moreover, in the last few decades, a general back-scaling of snow observation networks occurred worldwide [3]. The exploitation of satellite technologies can significantly help systematic snow monitoring [3]. Among the spaceborne remote sensing systems, synthetic aperture radar (SAR) instruments are particularly suitable for the analysis of snow deposits, providing data with spatial resolutions down to some meters, a global scale coverage and a few days revisit time [5], [6]. The centimetric wavelength of radar signal, capable of penetrating into snow layers, makes SAR spaceborne sensors particularly suitable for cryospheric applications [7], [8], [9], [10], [11].
The power and phase information at different frequencies and polarizations can be obtained from the specific processing of SAR images. Such information can be used to retrieve snow parameters and develop active-passive physical-electromagnetic response models. SAR backscattered power strongly depends on the air-snow and snow-ground interface reflections, as well as on the snow volumetric effects [8], [9]. Differential interferometric SAR (DInSAR) techniques have been successfully employed for This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ SWE retrieval, exploiting the phase shift between two sequential images and taking into account the effects of air-snow refraction and the snow-ground reflection [25], [26], [27], [28], [29], [30], [31]. The repeat-pass SAR interferometry can provide a high-resolution mapping of snow depth, but it requires a good coherence between two passes [27], [31]. These limitations for snow detection and mapping may be approached by combining SAR microwave observations with satellite optical spectrometer data at high-moderate spatial resolution [33], [34], [35], [36]. Spaceborne SAR polarimetry can also be exploited since the interaction between the polarized incident wave and the snow anisotropic features causes depolarization effects and specific scattering mechanisms [37], [38], [39], [40], [41]. Satellitebased passive microwave measurements may also be useful in combination with SAR observations, even though with a spatial resolution (tens of kilometers) much worse than the SAR typical resolutions [42], [43].
Most snow remote sensing applications from space have been focused on major mountain systems, such as Sierra Nevada, Alps, and Himalayan regions [14], [15], [29], [31]. The Italian Apennines, have never been considered for this kind of studies by the scientific community, probably due to the high spatial and temporal variability of their snow cover [44]. Nevertheless, the central Apennines region holds a key role for the meteorological dynamics in the Mediterranean area and, moreover, it hosts the southernmost European glacier, the Calderone. The latter is to be considered a glacieret rather than a glacier and its evolution represents a relevant indicator, at least for the medium latitudes, of the impact of climatic changes [45]. In such complex geographical regions, snow cover numerical models, coupled with meteorological weather forecast, can also be exploited to characterize the microphysics and dynamics of the snow cover [46], [47]. Snow cover model simulations can be used to provide auxiliary information in the satellite-based retrieval of the snow cover [48]. On the other hand, meteorological model predictions can be used to mitigate atmospheric water vapor effects on SAR interferograms [49], [50], [51], [52].
In this work, we present a DInSAR retrieval methodology to estimate SCA, height, and density in complex orography regions using Sentinel-1 satellite C-band data. The retrieval methodology is subdivided into classification and estimation subsequent stages. Classification is based on artificial neural networks (ANNs) and fuzzy logic. Estimation is based on an ANN trained by a forward DInSAR response model, coupled with snow cover numerical model outputs. In order to discriminate between wet and dry snow areas, data from satellite visibleinfrared imager are used. A digital elevation model (DEM) and a land cover database are also used to assist in the discrimination process. For snow cover classification, the results obtained from ANN-based retrieval methodology are compared with results from well-established decision threshold approaches. For SH and density estimation, the proposed methodology is compared, in terms of the obtained results, with the well-known analytical inverse method and two model-based statistical techniques (linear regression and maximum likelihood). The validation is carried out in the Italian Central Apennines, using in situ snow data collected between December 2018 and February 2019.
Results show that the ANN-based technique has a relatively good performance in terms of SCA classification accuracy, together with an estimation bias of about 0.5 cm and 5 kg/m 3 for SH and density, respectively. As expected, worse results are associated with low DInSAR coherence between two repeat passes and snow melting periods.
This article is organized as follows. Section II describes data and models including Sentinel 1 SAR data, in situ and auxiliary data as well as snow cover model data and DInSAR forward model. Section III focuses on the retrieval methodology investigating SCA mapping algorithms as well as SH and density inversion techniques. Section IV discusses the results in terms of SCA and height estimates. Section V presents the conclusions, whereas Appendixes give some details on refractive models, the DInSAR sensitivity analysis, winter 2018-2019 conditions, and the density estimation.
II. DATA AND MODELS SAR data are collected from the Sentinel 1 Copernicus mission, whereas ground data are those available from the area of interest (AOI) in central Italy Apennines. Auxiliary data from a land cover dataset, elevation models, and spaceborne optical radiometers are also used. The SAR interferometric model is based on a plane-wave assumption at the air-snow interface neglecting multiple scattering effects. Details are given in the next paragraphs.

A. AOI and In Situ Data
The AOI is the central Italy Apennines range between 41.5°-43.5°north latitude and 12.5°-14.5°east longitude. This region, bounded by the Tyrrhenian Sea on the west side and the Adriatic Sea on the east side, is characterized by a complex orography with deep valleys, medium-short rivers, small lakes, and relatively steep mountains [44]. The highest peak is the Corno Grande at 2912 m above sea level in the Gran Sasso Mountain range which hosts the Calderone glacier, the southernmost glacier (indeed, in the last decades transformed into a glacieret) of the European continent.
Within the AOI, there are snow measurement sites, used to validate SAR estimates of snow parameters. These measurements, provided by the Meteomont Avalanche Warning Service (MAWS), are daily acquired within a fenced area (to guarantee the integrity of the snow cover, characterized by requirements about orientation, exposure, and altitude) to allow an institutional snow mantle assessment for reference areas. The measurements provide SH h s , vertical thermal gradient, and density of the snowpack top layer, which we used to derive the mean snowpack density as described in Appendix D. Fig. 1 shows the location of measurement sites available in central Apennines for this study. The selected measurement sites in the AOI are 32, dislocated quite uniformly within the area itself, typically collecting data during the early morning.
In this work, the period of interest is the winter between December 2018 and February 2019. This winter in Central Italy was characterized by a first small snowfall in mid-December 2018, but the high temperatures during the second part of the  month led to a complete melt of the snow cover in most of the measurement sites. The largest snowfalls occurred during January 2019, thanks to the formation of a low-pressure system on the Mediterranean Sea. February 2019 was characterized by warm and dry conditions, which caused the retreat of the snow cover and the complete melting of the snow in most of the measurement sites. After February 2019 only snowfall at high elevations, higher than the measurement sites, were observed. For these reasons, we limit our analyses to the period between December 9, 2018 and February 25, 2019. Weather maps that describe the synoptic conditions of winter 2018-2019 are shown in Appendix C.

B. Sentinel-1 SAR and Auxiliary Data
As mentioned, C-band SAR data in SLC format from Sentinel-1A and Sentinel-1B platforms have been used for obtaining differential interferograms. The Sentinel-1 low-Earthorbit (LEO) satellite acquisition modes include: 1) Strip Map mode (80-km swath, 5 × 5 m 2 spatial resolution); 2) IW swath mode (250-km swath, 5 × 20 m 2 spatial resolution); 3) Extra-Wide swath mode (400-km swath, 20 × 40 m 2 spatial resolution); 4) Wave mode (20-km swath, 5 × 5 m 2 spatial resolution). For the purposes of our study, IW swath products coming from both platforms S1A and S1B are used. The nominal repetitiveness of Sentinel 1 constellation over the AOI is 6 days. Fig. 1 shows the Sentinel 1 swath, characterized by the relative orbit number 22, acquisition time at 05:11 GMT, a 250 km scan width with 5 × 20 m 2 resolution and vertical copolar and crosspolar (VV/VH) polarizations. Table I shows the list   TABLE I  SENTINEL-1 MASTER/SLAVE COUPLES WITHIN DEC. 9, 2018, AND FEB. 25,  2019. ACQUISITION TIME IS 05:11 GMT (RELATIVE ORBIT NUMBER 22) of Sentinel-1 acquisition dates for the SLC-IW images. This set is used to generate 13 interferograms within the period of interest. For all generated interferograms, the master and slave image coupling is based on adjacent dates. Both unwrapped interferometric phase maps and coherence maps are obtained by a standard process of filtering and unwrapping, as better explained later. For each Sentinel-1 SLC-IW image, the backscattering coefficient in both vertical copolar and crosspolar polarizations is also derived. The applicability of DInSAR technique to the snow cover parameter retrieval is strongly limited by the relatively poor repetitiveness of the Sentinel-1 constellation, which is indeed quite high if compared to other satellite SAR missions. Moreover, the DInSAR signal features shows ambiguities over complex orographic terrains and in presence of variable atmospheric conditions between master/slave acquisitions. To overcome as much as possible these SAR observation limits, auxiliary information and data are employed. Table II shows the list of external auxiliary data, used in combination with Sentinel-1 data. The table reports the following.
1) DEM, used to enhance the accuracy of the interferometric data processing and retrieval algorithms [5]. 2) Atmospheric correction model data, used to mitigate Din-SAR interferogram distortions caused by water vapor variable pattern between master and slave image dates [53]. 3) A land cover map, used to assist with dry and wet snow classification [54]. 4) Snow-cover products obtained from the Terra-satellite's "Moderate Resolution Imaging Spectroradiometer" (MODIS) instrument, used to assist in the discrimination of snow-covered areas from snow-free areas [14]. The simulation of fundamental snow cover properties is based on the Alpine3D model simulations. Alpine3D is a Lagrangian numerical model used to simulate high-resolution surface processes, and in particular snow processes in mountainous regions [44], [47]. The Alpine3D model is driven by external data such as measurements at stations or outputs from meteorological models. The Alpine3D model consists of several component modules.
1) The SNOWPACK one-dimensional (1-D) model of vegetation, snow, and soil [47]. 2) A 3-D radiation balance model (which considers shortwave scattering and longwave emission from terrain and tall vegetation). 3) A 3-D drifting snow model solving a diffusion equation for suspended snow and a saltation transport equation. Examples of the use of the Alpine3D model include the analysis of snow cover dynamics for avalanche warning, studies on permafrost development and vegetation changes for the monitoring of climate change evolution, as well as the estimation of soil moisture aimed at obtaining highly accurate meteorological and flood forecasting models. The Alpine3D model, driven by forecasted weather data, has been validated in the AOI [44]. Results showed that Alpine3D is able to reproduce the observed SH and SWE with a bias of 6 cm and 35 mm respectively, with a correlation coefficient of 0.87 for SH and 0.74 for SWE.
The Alpine3D map dataset has been generated with a 3-km spatial resolution at 1-h temporal resolution and has been resampled to a reference resolution of 100 m, which has been adopted in this study to ensure uniformity across all employed spatial data. We deem such a resolution to be sufficient to characterize snow properties for a relatively large AOI, and to be an optimal tradeoff in terms of computational costs [1], [3]. For each pixel of the AOI, the Alpine3D model provides altitude, SH, SD, and snow liquid water fraction as well as the temperature profile. In this work, the generated Alpine3D dataset is basically used to train snow estimation models from DInSAR data.  2. Geometry of DInSAR signal propagation in a snow deposit, characterized by a single uniform snow layer. The yellow ray path refers to the snow absence, whereas the red one to the snow presence inducing refraction. Planar interfaces are assumed, whereas the ground plane is tilted due to the ground local slope where θ is the incidence angle, θ l is the local incidence angle, θ t is the terrain slope angle, and θ s is the Snell refraction angle, Δh s is the relative SH, Δz s is the relative snow thickness. Other symbols are defined in the main text, apart from g that represents the Earth gravitational acceleration vector.

C. SAR Differential-Phase Forward Model
The snow estimation methodology from SAR differential interferometry, based on satellite repeat passes above the same scene, takes advantage of the nearly linear dependence existing between the SAR signal path delay and SH h s [7], [11]. Even though widely explained in literature, we briefly summarize here the principles to introduce the basic notation. Table III lists the main definitions, acronyms, and symbols, used throughout the text. Fig. 2 illustrates how the SAR signal, supposed to be locally approximated by a plane wavefront, is affected by the snow refraction causing a path delay with respect to the same surface without snow. Note that the incidence angle θ is the angle between the normal of the flat earth (ellipsoid) and the satellite pointing, whereas the local incidence angle θ l is the angle between the normal of the topography and the satellite pointing. Depending on the relief within the scene and the topography slope angle θ t , these two angles can be significantly different (if θ t = 0 in Fig. 2 then θ l = θ). Each interface in Fig. 2 is considered as a specular plane, for simplicity, and possible multiple scattering effects within the snow layer (due to ice grains, ice lenses, and/or vertical refrozen structures larger or comparable with SAR central wavelength) are disregarded [27], [29].
In Fig. 2, we exemplify the problem geometry for one planeparallel uniform non-absorbing isotropic snow layer above a surface with a local slope [25], [29]. The snow complex permittivity can be expressed by [13] with ε Rs and ε Js the real and imaginary part of ε s (with j the imaginary unit), ε rs the relative snow permittivity with ε Rrs and ε Jrs the real and imaginary part respectively, and ε 0 the vacuum permittivity. The adopted snow permittivity models, mainly depending on the dry SD ρ ds and liquid water fraction f w as highlighted in (1) to include wet snow, are reported for completeness in Appendix A [55], [57].
Note that the SD ρ s differs from the dry SD ρ ds due to the presence of liquid water fraction. Moreover, snowflakes are usually compressed by weight after ground deposition (snow settling) and some days later, recrystallization may form ice grains that are spherical at first, vertically extended at later phases, and which typically increase in size over time (snow metamorphism) [55], [56]. This snow microstructures typically cause macroscopic anisotropy that reflects into polarization-dependent effective dielectric constant of the air-ice mixture, perturbed by liquid water from the melting process. As a result, the snow complex permittivity in (1) can be polarization dependent and snow can be modeled as an air medium populated by aligned spheroidal ice particles using the Maxwell-Garnet theory [26].
The point P on the ground is covered by a snow layer whose relative height is (see Fig. 2) with h s the SH at time t 2 (>t 1 ) and h s0 the reference surface (ground surface at time t 1 ) height, whereas the quantity Δz s in (2) is the relative snow thickness, measured along the terrain normal. In (2), θ t represents an important feature in complex orography. If R s is the distance between the SAR sensor and snow surface, the plane wave fronts of the signal will be refracted at point Q and, covering the distance ΔR s , will reach the point P.
In absence of snow, no refraction would occur and propagation distance in air to point P is R a +ΔR a . If n s = √ ε Rrs is the snow refractive index real part, the geometric path difference ΔR between the two conditions (i.e., snow presence and absence) is determined by being ΔR s = Δz s /cosθ s and ΔR a = ΔR s cos(θ l -θ s ). By applying trigonometry and using the Snell law sin θ = n s sin θ s [13], we obtain Changes in snow layer thickness Δz s by a positive (due to snowfall) or negative (due to snow ablation) amount reflect into a change of the SAR wave front two-way differential phase at pp co-polarization given by ΔΦ pps = 2 (2πΔR/λ). Using (2), it holds for the differential phase due to snow Δh s (5) being k = 2π/λ with f the SAR operating frequency and λ = v s /f 0 the SAR wavelength in the snow layer propagating at phase velocity v s = c/n s with c the vacuum light velocity. For Sentinel-1 DInSAR, p = v using the copolar returns at vertical (v) polarization. Note that for an anisotropic snow layer, due to its polarization-dependent effective permittivity, ΔΦ vvs is not equal to ΔΦ hhs thus causing a copolar differential phase shift (or phase difference) Δ Φ cos = ΔΦ vvs − ΔΦ hhs , fairly well correlated at X band to the h s [26]. In the following text, we will keep the notation ΔΦ pps for generality.
The previous relation (5) constitutes the basic forward model that can be used to translate the dry or wet snow thickness into DInSAR differential phase ΔΦ s . For a dry snow uniform isotropic layer, using the permittivity model (A.1) illustrated in Appendix A, it can be shown that (3) can be approximated with an error less than 3% (for θ less than 50°) by [27] ΔΦ pps ∼ = ka c b c +θ being the terms a c = 1 (average optimal value) and b c = 1.59, whereas the relative SWE height Δh swe is defined as follows: with Δz swe the relative SWE thickness (see Fig. 2). The last term in (7) holds only for a uniform snow layer. Moreover, the approximate expression in (6), in terms of Δh swe , holds also for a non-uniform snow layer with a variable density below the air-snow interface, i.e., [27] (8) where N is the number of the different snow ith layers. For a wet snow layer, ΔΦ s becomes a function of the snow water fraction so that (6) must be replaced with the more general expression (5) [29], [40]. As an example, since n s is a function of ρ ds and f w , using the Wiesmann-Mätzler snow permittivity model in Appendix A, Fig. 3 shows how ΔΦ s is affected by different values of ρ ds and f w for two incidence angles θ = 20°a nd θ = 40°, assuming θ t = 0 and setting the frequency at 5.6 GHz, the one of Sentinel-1 SAR. ΔΦ pps is less sensitive to h s variation for low snow densities and lower incidence angles, even though the impact of the latter is relatively low. The increase in liquid water fraction enhances the sensitivity in all cases. For Δh s larger than a few centimeters, we expect a SAR phase folding issue at C band (that is, ΔΦ pps >180°).
The analytical dependence between ΔΦ pps and Δh s in (5) can also be used for carrying out an uncertainty analysis to evaluate how errors on θ l , ρ ds , and f w would influence ΔΦ pps . Using the first-order error propagation theory and assuming the statistical independence of all uncertainties, the DInSAR differential phase standard deviation σ ΔΦ pps can be expressed by the square root of the error variance sum weighted by the respective square partial derivative being σ θ l , σ ρ ds , and σ f w the standard deviations of θ l , ρ ds , and f w , respectively. The three partial-derivative terms of (9) are analytically derived and discussed in Appendix B using, as an upper bound, the linear-mixing model for the snow relative permittivity described in Appendix A.

III. METHODOLOGY
The previous section was dedicated to the description of data and DInSAR forward models. Starting from that analysis, this section will be devoted to the illustration of the proposed inversion approaches that will be applied in the next sections.

A. Retrieval Methodology
The adopted processing methods include DInSAR standard processing schemes, inversion methods to obtain snow properties estimates from SAR sensor measurements, decision tree, and validation steps. Three main branches are visible in Fig. 4.
For the retrieval processes, we have opted for two separate ANNs as inversion methods [60]: one for the snow classification and one for the SH and density estimation. Indeed, the diagram of Fig. 4 can be easily adapted, with minor modifications, to represent other retrieval methods, for instance, based on analytical and statistical approaches as discussed later. The overall ANN-based retrieval scheme produces for the whole area.
3) SD or SWE estimated maps. Each step is discussed in what follows. To properly handle data fusion, Fig. 4 shows the SAR-based integrated approach for the ANN retrieval of SCA, type, height, and density, coupling Alpine3D model and auxiliary data. Input sources include the following.
1) Data from satellite SAR.
2) Data from auxiliary remote sensing observations such as spaceborne optical and near-infrared sensors (these data are used for the training of the ANN, which performs the snow cover classification and also to validate results). 3) Data products derived from satellite data such as DEMs and land cover classification maps (these data are used within the ANN to improve the precision of the results). 4) Data from numerical models for snow cover evolution simulations (these data are used for the physical-analytical snow model, which requires the knowledge of some data not retrievable with satellites like SD and liquid water fraction). 5) Ground-based data (e.g., in situ measurements: these data are used to validate results). As mentioned, all data listed in Tables II and III have been reprojected and bilinearly interpolated onto a common resample grid with a linear resolution of about 100 m [1], [3]. As discussed in [27], the quality of the measured interferometric coherence γ ppsm between two consecutive acquisitions is crucial for the DInSAR method application. Incoherent gaps between interferograms tend to degrade the overall DInSAR signal-to-noise ratio. The possible coherence loss on snowcovered areas may be caused by the variation of the liquid water fraction due to positive temperatures, melting processes, and/or rainfall. Coherence loss may also occur when strong winds lead to a snow erosion and redistribution. Another possible cause of the coherence loss is the snow structure collapse by gravity or external forces over weak or thin snow layers. Finally, the vegetation cover above the snow cover can also degrade the DIn-SAR coherence due to the snow weight on the plant branches. Using DInSAR lower frequencies, the coherence temporal loss is slower, the penetration depth is larger due to less impact of liquid water fraction and phase-wrapping errors are usually smaller, but SH and density retrieval is less accurate due to longer wavelengths. In this respect, the use of C-band SAR data may be a good compromise.
In DInSAR, the measured differential phase ΔΦ ppsm at snow surface is computed from the measured absolute phase Φ ppsm of 2 repeated passes at time t 1 (master image) and t 2 (slave image), that is [14] where Φ ppsm is the two-way propagation delay due to the radar wave refraction in the snow cover, referred to snow-free conditions (see Fig. 2).
In the flow diagram of Fig. 4, master and slave SAR SLC measurements are processed according to a consolidated scheme. Assuming a sufficiently high coherence between two SAR acquisitions over the same snow-covered scene, the measured snow interferometric phase ΔΦ ppsm (with the subscript m standing for measured in the following text) at pp polarization is obtained for the (x,y) pixel from the measured interferometric phase ΔΦ ppm by (e.g., [12], [13]): where ΔΦ ppf lat and ΔΦ pptopo are the differential phase differences due to changes of the relative distance between satellite and target for flat Earth and for topography, respectively, whereas ΔΦ ppatm is generated from changes in atmospheric propagation due to ionosphere and troposphere, and ΔΦ ppnoise is differential phase noise.
2) Orbit file application and geocoding.
3) Enhanced spectral diversity and interferogram generation. 4) TOPSAR debursting. 5) Topographic phase removal. 6) Goldstein phase and multi-look filter. 7) Masking and phase unwrapping. 8) Terrain correction. The phase unwrapping is based on the Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping algorithm, which is an effective and efficient two-dimensional phase unwrapping algorithm [63]. The atmospheric phase screen correction, derived from the Generic Atmospheric Correction Online Service for InSAR (GACOS) products, has been applied to the unwrapped phase values [53]. The SAR local incidence angle information, required in (5), is obtained as a subproduct of the interferogram terrain correction process. The final output contains Sentinel-1 unwrapped phase interferograms, coherence, and local incidence angle maps. SAR SLC data are also processed to obtain the SAR backscattering copolar coefficients.

B. SCA Classification
For the SCA mapping, we need to address two separate objectives.
1) Identifying the presence or absence of snow (primary binary mask). 2) Discriminating the wet snow from dry snow within the SCA (secondary binary mask). When using only SAR backscattering data, the Nagler-Rott empirical approach can be adopted [10], [28]. This technique uses the SAR backscattering coefficient σ 0 to discriminate between wet and dry snow. For each pixel (x,y), two backscattering measurements are considered: 1) σ 0 vvmref (x,y), the copolar vertically polarized data from a reference SAR image; 2) σ 0 vvm (x,y), the copolar vertically polarized data from the SAR image during the snow cover period.
According to the Nagler-Rott method, the copolar ratio r vv (x,y) indicates the presence of wet snow in case its value is below a threshold, and the presence of dry snow otherwise [10] where r 0wet is the wet-snow threshold (in dB), very often empirically set to about −3 dB. A further extension of (12) is considering also the cross-polar backscattering coefficient σ 0 vvm (x, y) and θ l so that the combined ratio r c is given by [28] where W is an empirically derived angular weighting tree function [28]. To identify the wet snow cover, the ratio r c (x, y, θ l ) is typically set below the threshold of −2 dB.
A limitation of the Nagler-Rott approach is the exploitation of SAR backscattering data only. From Sentinel-1 DInSAR acquisitions, we can also extract the vv-copolar differential phase difference ΔΦ vvsm and coherence γ vvsm as well as integrate auxiliary data available from Table II. To perform an effective data fusion, we have designed an ANN and combined it with a fuzzy logic algorithm to obtain snow/no-snow and dry-snow/wet-snow cover maps.
The proposed Snow Cover Classification ANN (SCC-ANN) employs a feed-forward backpropagation architecture with a Bayesian regularization. The SCC-ANN is composed by a single hidden layer with 14 neurons using a sigmoid activation function. Formally speaking, we can write for each pixel (x,y) wherex scc (x, y) is the estimated snow cover classification parameter vector, whereas y mSAR (x, y) and y ap (x, y) are the input data parameter vectors extracted for each pixel from Sentinel-1 SAR observations and a priori auxiliary data in Table II, respectively. Specifically, 1) The input measurement vector y mSAR (x, y) includes: i) DinSAR coherence γ vvsm (x,y); ii) SAR copolar backscattering coefficient σ 0 vvm (x,y) for master and slave images; iii) SAR backscattering cross-polar coefficient σ 0 vhm (x,y) for master and slave images; iv) SAR backscattering copolar coefficient σ 0 vvmref (x,y) at the reference (nosnow) date; v) SAR backscattering cross-polar coefficient σ 0 vhmref (x,y) at the reference (no-snow) date.
2) The input vector y ap (x, y) includes: i) local incidence angle θ l (x,y); ii) SRTM DEM data h topo (x,y); iii) CORINE land cover labels L c (x,y).
3) The output vectorx scc (x, y) includes two snow parameters: i) the fractional snow cover (FSC) f SC (x,y), that is the percentage of snow cover for each pixel; ii) the volumetric liquid water fraction f w (x,y), that is the percentage of liquid water in a unit volume of snow. To train SCC-ANN, the output variables f SC (x,y) and f w (x,y) have been extracted from MODIS and Alpine3D respectively, after a bilinear interpolation and resampling onto a common grid at 100-m resolution. By selecting more than 10.000 samples within the area/period of interest and dividing the whole set into a validation (15%) and test (15%), the SCC-ANN shows relatively good results, as discussed later.
Based on the SCC-ANN estimatedf SC (x, y), we can perform the first task to identify the presence/absence of snow (primary binary mask). This can be accomplished by a simple FSC threshold algorithm, like the one of Nagler-Rott in (12), that isf where f 0snow is an arbitrary threshold (in %), set to an optimal value during the training stage (typically between 5% and 25%). The SCA is finally computed by summing up all the snow pixels satisfying snow presence in (15). We can perform the second task to discriminate the dry/wet snow (secondary binary mask) using the SCC-ANN estimated f SC (x,y) andf w (x, y). Instead of relying on a fixed threshold for the fraction of water, we have adopted a fuzzy logic approach also involving the snow surface temperature T s (x, y) as a third parameter coming from the Alpine3D simulations. Fuzzy logic is a form of many-valued logic in which the truth value of variables may be any real number between 0 and 1, extending de facto the Boolean logic [62]. Fuzzy logic models or sets are an extremely convenient and suitable form of representing uncertain information. Two membership functions M w and M t are defined for the snow water fraction and snow surface temperature, respectively: where in our case the coefficients are w 1 = 0.1%, w 2 = 0.8%, If the inference function I is defined through a multiplicative formula by means of then the following defuzzification rule is applied to the inference function: where the threshold I 0 is usually set to 0.5, equal to a uniform probability of the 2 fuzzy states.

C. SH and Density Estimation
Once the SCA and wet/snow mask is obtained, we can approach the estimation of Δh s (x, y) for each image pixel from SAR and DInSAR observations. In this respect, satellite optical data and products are not useful due to the high snow and cloud extinction at those wavelengths. Microwave remote observations are the only ones that can provide a physically based estimate of the snow depth.
The artificial Neural Network Inversion (NNI) method is based on the use of an ANN, as already done in Section III-B. The snow depth estimation operation is performed by a dedicated neural network, named SnowPack Estimation ANN (SPE-ANN). The architecture of the SPE-ANN is like the SCC-ANN one, using a Bayesian regularization with 1 hidden layer and 14 neurons.
Formally speaking, we can write for each pixel (x,y) wherex spe (x, y) is the estimated snow cover parameter vector, whereas z mSAR (x, y) and z ap (x, y) are the input data parameter vectors extracted for each pixel from SAR observations and a priori integrative data in Table II, respectively. Specifically, for Sentinel-1 SAR data 1) input vector z mSAR (x, y) includes: i) DinSAR coherence γ vvs (x,y); ii) DInSAR co-polar unwrapped differential phase ΔΦ vvs (x, y); iii) DInSAR cross-polar unwrapped differential phase ΔΦ vhs (x, y). 2) input vector z ap (x, y) includes: i) local incidence angle θ l (x,y); ii) atmospheric water vapor correction maps ΔΦ atm (x, y), derived from the GACOS algorithm; iii) SRTM DEM data h topo (x,y); iv) CORINE land cover labels L c (x,y). 3) output vectorx spe (x, y) includes two snow parameters: i) snow relative height Δĥ s (x, y); ii) SDρ s (x, y). To train SPE-ANN, the output variables Δĥ s (x, y) and ρ ds (x, y) have been extracted from Alpine3D, after a bilinear interpolation and resampling onto a common grid at 100-m resolution.
To intercompare NNI with other estimation approaches, we have implemented the following inversion algorithms: A brief description of the these inversion algorithms is given hereafter.
1) The PAI algorithm is obtained by analytically inverting the DInSAR equation in (5), that is: The PAI inversion algorithm requires that θ l and θ t are known, but more importantly an assumption on ρ ds and f w to estimate ε Rrs in (21) (see Appendix A). The latter parameters are not usually known, and their arbitrary assignment constitutes the main limitation of this inversion method.
2) The LRI algorithms are based on the statistical linearregression approach, described in the following equation: where the regression coefficients a h0 and a h1 are derived from the DInSAR forward model simulations. The latter can be based on: 1) a MonteCarlo random approach, where the input model parameters (i.e., h s , ρ ds , f w in (5) using (A.2)) are uniformly varied between a minimum and a maximum value (the algorithm is then named LRI-MC); 2) Alpine3D numerical outputs, where the same input model parameters are provided by the Alpine3D simulations described in Section II-D (the algorithm is then named LRI-A3).
The main difference between LRI-MC and LRI-A3 is that the input snow parameters are randomly selected in the first case and physically constrained in the second case. This implies that there are Monte Carlo triplets that do not have necessarily a physical meaning.
3) The MLI algorithms are statistical approaches where the error probability density function is maximized, that is under a Gaussian hypothesis its negative argument is minimized [61]. If a forward model is available, minimizing the error (or objective) function implies to minimize the square difference between the forward model observable ΔΦ pps [derived from (5)] and the corresponding measured one ΔΦ ppsm . Formally speaking, using the DInSAR observable ΔΦ pps (x, y), we can write the estimation through the objective function as where argmin is argument-of-the-minimum operator restituting the elements Δĥ s andρ ds of the domain at which the function values are minimized. With respect to PAI and LRI algorithms, the MLI technique can provide a statistical estimate of both SH and density. As for LRI, we can distinguish again between a Monte Carlo random and an Alpine3D model approach, defining MLI-MC and MLI-A3 respectively. For computational purposes, the minimization can be carried out by searching in a pregenerated database and setting a confidence interval to provide an associate uncertainty.
For h s estimate, we can further investigate a method combining MLI and Alpine3D simulations. This combined approach, called MLI-A3C, can use Alpine3D model values of SH whenever the DInSAR coherence ρ vvs is below a given threshold (e.g., 0.3). This approach should improve the accuracy of the estimated Δĥ s when the DinSAR coherence is too low to provide reliable estimates of SH. The combination between MLI estimates and Alpine3D model values can be linearly weighted by the coherence γ vvm as follows:

IV. RESULTS
After introducing the retrieval methodologies, in this section, we will illustrate the numerical results obtained from the forward and inverse models.

A. Snow Cover Model Simulations
An example of a SH variation map, simulated from the Alpine3D model between Dec. 27, 2018 and Jan. 2, 2019, is shown in Fig. 5. The reason for the snow depletion in many high-altitude areas is due to the lack of snowfall and temperature raising during the considered week. This h s variation map is a useful representation when compared with DInSAR interferograms. Fig. 6 shows bi-dimensional histogram of dry SD and volumetric liquid water fraction, generated by the Alpine3D snow cover model within the area (see Fig. 1) and period of interest. SD values range from 50 till 400 kg/m 3 , whereas f w spans from 0% up to 6% denoting relatively dry samples, probably due to the exclusion of the last winter months. Note that a minimum threshold of h s equal to 1 cm is applied to plot data in Fig. 6 (smaller values are considered not numerically significant). The DInSAR forward model in (5) can be coupled with the numerical outputs of the Alpine3D snow cover model for the area and period of interest. Fig. 7 shows the simulated ΔΦ pps values using the Alpine3D dataset of Δh s , f w and ρ s . The snow cover model constrains the variability of the DInSAR differential phase as a function of the simulated variability of the physical snow parameters. Moreover, the linear correlation between ΔΦ pps and Δh s is apparent from Fig. 7. This correlation can be used as a precious a priori information within the inversion process.

B. SH Estimation Expected Errors
A numerical intercomparison of the snow relative height estimation algorithms can be carried out on simulated data. Fig. 8 shows the comparison between the estimated and "true" SH using the Monte Carlo database for training (50% data) and test (50% data). NNI shows a standard deviation error of about 5.7% (with about zero mean) and an estimate correlation of 0.98. The PAI technique with an average SD of 250 kg/m 3 shows a better estimate correlation (about 0.96) than using an extreme value of 500 kg/m 3 (0.92, not shown). The two inversion algorithms LRI and MLI appear to be equivalently performing in terms of estimate correlation (about 0.96), mean error, and error standard deviation. These results confirm that the PAI approach is prone to be more biased and inaccurate than other retrieval techniques.
The same analysis can be repeated using a realistic snow scenario, provided by the Alpine3D model, instead of the Monte Carlo random one, as shown in Fig. 9.

C. SCA Classification
The dry/wet snow classification map, estimated by the proposed SCC-ANN method, is shown in Fig. 10. With the purpose of providing an independent reference, we have performed the same classification using the Nagler-Rott algorithm, described in (13), as shown in Fig. 11. A partial, but noticeable, correspondence for the wet snow areas can be observed by comparing the two maps, with a predominance of wet snow when applying the SCC-ANN retrieval.    By applying the NNI technique to data of Figs. 10, 12 shows the Δh s map, estimated on Jan. 2, 2019. For the considered couple of dates, a validation has been performed comparing results with the 32 in situ station measurements acquired on Jan. 2, 2019. The comparison indicates a correlation of 0.77 for samples whose DInSAR coherence was found to be above 0.4 (corresponding to a subset of 9 samples on 32 total samples available for that date). SCA, estimated from the proposed SCC-ANN method, can be validated by using MODIS satellite products expressed in Fig. 10. Example of dry/wet snow cover map estimated by SCC-ANN on Jan. 2, 2019. The AOI is the same as in Fig. 1 and the snow area extent is the one derived from SCC-ANN algorithm. Fig. 11. Example of dry/wet snow cover map estimated on Jan. 2, 2019 with the algorithm proposed by Nagler and Rott (2016). Magenta color indicates wet snow, whereas cyan indicates dry snow. The AOI is the same as in Fig. 1 and the snow area extent is the one derived from SCC-ANN algorithm.
terms of SCF at about 500-m resolution. This intercomparison is quite cumbersome since MODIS-derived snow cover, even though post-processed to avoid cloud coverage, shows some residual errors that can lead to a significant underestimation of the SCA. These effects may be dependent on the SCF percentage threshold. Within these limitations, we have performed an intercomparison between MODIS SCF and SCC-ANN SCA estimates using 50% of the available dataset. Apart from probability of detection (POD) and false alarm ratio (FAR), we can also use the classification accuracy ACC derived from the confusion (or contingency) matrix elements, as follows [60]: where TN represents the number of True Negatives, TP represents the number of True Positives, FN represents the number of False Negatives, and FP represents the number of False Positives. Using a MODIS SCF threshold of 1% as a reference, numerical results show that, for the SCC-ANN neural network approach, POD is 0.451, FAR is 0.04 and ACC is 0.78. When increasing the MODIS SCF threshold to 20% and 50%, POD becomes, respectively, 0.55 and 0.60, FAR is 0.11 and 0.24 and ACC 0.83 and 0.87. As mentioned, the POD should be read with caution, but the accuracy is quite high and the FAR is relatively low.

D. SH and Density Estimation
Considering the whole period from Dec. 9, 2018 till Feb. 25, 2019, we have at disposal 211 in situ measurements of SHs and densities within the AOI, shown in Fig. 1. In order to validate results, the SCC-ANN and SPE-ANN neural network retrieval methods can be applied to retrieve Δh s and ρ s in correspondence of each ground snow-measurement station.
To quantify the comparison, 4 further error statistical indexes are used to evaluate the overall performance during the considered season. The estimate error can be defined as where Δh sm (x, y, t) and Δĥ s are the measured and estimated relative SHs, respectively, at the location (x,y) and time t. The error metric consists of the error bias b δ , the root mean square where the angle brackets represent the ensemble averaging over time (number of interferograms) and space (snow measurement sites), where σ hm and σ h are the standard deviations of Δh sm and Δĥ s , respectively. Moreover, the index of agreement (IA) parameter can be useful to provide a comprehensive assessment of the degree of similarity between estimates and measurements,  Table II), obtained as a mean of all available in situ data and from all the snow-height retrieval algorithms compared in Fig. 13  and is defined as [64] i δ = Δĥ s (x, y, t)−Δh sm (x, y, t) 2 Δĥ s (x, y, t)− Δh sm (x, y, t) +|Δh sm (x, y, t)| 2 (28) being 0≤i δ ≤ 1 where 0 indicates no agreement and 1 indicates complete match.
Using the in situ validation set, Fig. 13 shows a bar diagram of the performance indexes b δ , r δ , ρ δ , and i δ for the Δh s estimation, according to different inversion algorithms. The diagram includes six of the inversion approaches, described in Section III and compared in this study, that is, SPE-ANN (neural network) against PAI (physical analytical), LRI (linear regression), MLI-MC (maximum likelihood trained with Monte Carlo), MLI-A3 (maximum likelihood trained with Alpine3D), and MLI-A3C (maximum likelihood trained with Alpine3D and combined with Alpine3D). Fig. 13 shows that, for the h s retrieval, the NNI estimation method provides the best performances for all considered four error indexes with a bias of about 0.7 cm, RMSE of about 17 cm, correlation of about 0.5, and IA of about 0.6. The second-best approach is the combined MLI-A3C followed by MLI-A3 and MLI-MC with the PAI inversion algorithm showing the worst estimation capability with a bias of about −2 cm, RMSE of about 30 cm, correlation of about 0.05, and IA of about 0.4. The behavior of PAI is indeed expected to be the worst due to the assumptions that have to be necessarily made on the SD, as discussed previously.
Results for the SD are similarly shown in Fig. 14, even though only three algorithms, providing SD retrieval, are compared (namely, MLI-MC, MLI-A3, SPE-ANN). For the SD estimate, the results obtained from the SPE-ANN neural network approach are again better than maximum-likelihood, with an IA of about 0.85. Fig. 15 finally shows the average time-cumulated snow as a function of available overpass dates, obtained as a mean of all available in situ data and from all the SH and density retrieval algorithms compared in Figs. 13 and 14. Relative SH retrievals are converted into snow absolute height by starting from the zero value at the initial date on December 9, 2018. Standard deviations of all estimates at each site at a given date for SH and density are also shown in Fig. 15.
These plots visually confirm what has already been discussed about the capability of ANN-based approaches to better follow the accumulation and ablation of the snow deposit due to snowfall and melting and/or mobilization factors, respectively. In terms of estimate correlation with in situ measurements, the SPE-ANN technique shows a value equal to 0.49 using all 211 total available samples (i.e., when no subselection based on a coherence threshold is performed) that becomes about 0.55 for the sample subset whose DinSAR measured coherence is above 0.4 (corresponding to a subset of 42 samples on 211 samples totally available).

V. CONCLUSION
This work has been devoted to developing a new approach for the estimation of snow cover extent, SH, and SWE, which combines C-band SAR Sentinel-1 satellite data coupled with numerical snow cover model data using ANNs. The AOI is within the Italian central Apennines, which is novel to this kind of studies. Besides the SAR C-band data, auxiliary data sources adopted for this study include optical sensor data (MODIS), atmospheric correction data (GACOS), land cover data (CORINE), and DEMs.
An ANN-based approach has been developed for both wet and dry SCA classification and estimation. Such approach has been compared with other physical and statistical inversion approaches. For the SCA classification, the proposed ANNbased method has been combined with a fuzzy-logic scheme to discriminate between dry and wet snow. For SH and density estimation, the proposed ANN-based method has been compared with a physical-analytical method, linear regressionbased methods, and maximum likelihood methods. Synthetic forward model simulations obtained using a uniform probabilistic distribution (Monte Carlo database), as well as using snow cover model simulations (Alpine3D database), have been used for training purposes in the SH and density estimation methods.
The AOI, in the Italian central Apennines, has been analyzed with a spatial resolution of 100 m using three months of the 2018-2019 winter. SCA retrieval has been validated using MODIS satellite snow cover maps, obtaining accuracies from 0.78 up to 0.87 when discriminating between snow-covered areas and snow-free areas. Maps of dry and wet snow cover obtained from the SCC-ANN have been compared with analogous maps obtained from the Nagler-Rott classification methodology, noting a substantial similarity in SCAs, with more wet snow subareas detected when using the SCC-ANN retrieval. SH estimates have been validated using measurements from more than 200 in situ observations. Imposing a DInSAR coherence larger than 0.4 (corresponding to 42 samples out of a total of 211 samples), we have obtained an estimate correlation value of 0.55. The ANN-based overall error has a bias and a root mean square equal to about 0.5 and 20 cm for SH and to 5 and 80 kg/m 3 for SD. Worse results mainly depend on the low DInSAR coherence values, as well as on snow melting effects and meteorologically driven changes of snow state and accumulation.
The main advantages of a DInSAR-based retrieval of snow cover are the achievable spatial resolution, the availability of data in all weather conditions, and the capability of snow layer probing (depending on the SAR frequency). On the other side, the main DInSAR drawbacks consist of relatively small swaths and long revisit periods, causing a coherence loss in the interferograms and, ultimately, unreliable estimates.
However, the growing interest in SAR satellite missions and the current trend in deploying increasingly large satellite constellations could gradually reduce this limitation up to acceptable levels for snow cover retrieval applications.
Future developments of this work will include the extension of the proposed approach to other SAR missions, possibly not only at C band, but also at L and X band. By adopting a multifrequency and multimission SAR approach, we might overcome, totally or in part, some of the limitations we encountered in the present work.
In general, especially during the melting periods, snow is not dry due to the presence of a liquid water fraction. A simplified frequency-dependent model of snow permittivity is based on the linear mixing (LM) approach that, for the (wet) snow, can be formulated as follows [12]: where f w is the volumetric liquid water fraction (%), f is the frequency, and ε Rrw is the liquid water relative permittivity following the Debye model [13]. At C band, the pure water value of ε Rrw is about 66.7 (see later). The LM model of ε Rrs can be compared with the one from Wiesmann and Mätzler (WM) whose general formulation applies also in wet snow conditions. The WM model exploits the water relaxation spectrum so that the frequency-dependent snow relative permittivity ε Rrs , for small values of f w , can be decomposed in four additive terms, one of them being the expression for dry snow in (A.1), that is [55] ε (W M) Rrs 3) being Re the real part operator with ε Rra , ε Rrb , and ε Rrc the Debye frequency-dependent complex terms. The latter ε Rrk (k = a,b,c) is expressed by where the Debye parameters ε ∞k , ε sk , and f 0k are depending on the depolarization factors. Details can be found in (53)-(57) of the WM-cited paper [55]. Another well-known frequency-dependent snow permittivity model is the one proposed by Sihvola and Tiuri (ST) [58], [59]. Based on measurements around 1 GHz scaled to higher frequencies by using the Debye formulas, they proposed for the snow relative permittivity real part ε Rrs the following formula [59]: being c 0 = 1, c 1 = 1.7, c 2 = 0.7, c 3 = 0.1, c 4 = 0.8 with w 0 = 4.9, w 1 = 82.8, and f 0 = 8.84 GHz related to a relaxation time of 18 ps [note that in (A.6) ρ ds is expressed in g/cm 3 and f in GHz]. Fig. 16 shows a comparison between the LM simplified model and the more accurate WM and ST models for an average (250 kg/m 3 ) and higher (500 kg/m 3 ) dry SD. All models show a substantial linear dependence on small values of f w (between 0 and 4%), being the LM providing values of about 60% larger than WM and ST for water fraction larger than 3%. For low values of f w (i.e., <1%), the overestimation of LM relative permittivity reduces to less than 25%.
To evaluate the WM and ST approximations of snow relative permittivity, we can use the Alpine3D snow cover model simulations which provide, among other snow parameters, both  using the Alpine3D outputs for the area and period of interest. In this realistic scenario, dry SD values range from 50 to 400 kg/m 3 , whereas the water fraction ranges from 0% to 5% for most samples. The discrepancy average and standard deviation between the two snow permittivity models in (A.3) and (A.5) are 0.048 and 0.098, respectively, that is ST and WM models differ less than 3% on average. The overestimation of LM permittivity model with respect to WM and ST models is (not shown) 0.28 and 0.33 on average, respectively, that is about 16% and 18% on average.

B. Differential Phase Uncertainty Analysis
Starting from (5), we can use the LM snow permittivity model for a dry SD less than 400 kg/m 3 in (A.3) to derive the following expression for the DInSAR differential phase, in the case of dry snow, for θ t = 0°Δ (B.1) For a wet snow condition, we can employ the LM model expressed in (A.2), a choice that may represent an upper bound of snow permittivity values (see Fig. 16).
We can then generalize the previous expression for dry snow as Using (B.2), the partial derivatives of the three terms of (9) in Section II-C are the following: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (1−f w )(1+a 1 ρ ds +a 2 ρ 3 ds ) + ε Rrw f w − sin 2 (θ) ∂ΔΦ pps ∂f w = −kΔh s × 1+ a 1 ρ ds +a 2 ρ 3 ds − ε Rrw (1−f w )(1+a 1 ρ ds +a 2 ρ 3 ds ) + ε Rrw f w − sin 2 (θ) (B.3) Substituting (B.3) into (9), we get the ΔΦ pps uncertainty expression as a function of the key parameters θ l , ρ ds , and f w . Fig. 18 shows the sensitivity of ΔΦ pps as a function of Δh s for lower and higher values of θ l , f w , and ρ ds . The curves are linear Fig. 18. Uncertainty of DInSAR differential phase ΔΦ pps versus snow relative height, using (9) varying the liquid water fraction (top panels with low value f w = 0.1%, bottom panels with large value f w = 4%) and local incidence angle (left panels with lower value θ l = 20°, right panels with higher value θ l = 40°). In (9), the standard deviations are set to σ θ l = 3°, σ fw =0.1 f w , σ ρ ds = 100 kg/m 3 . and symmetric with respect to Δh s due to its proportionality from (B.3). For a given incidence angle, a larger liquid water fraction increases the differential phase uncertainty, especially for lower SD. For a given SD, the differential phase is larger for larger incidence angles. For larger water fraction, the effect of the dry SD can be considered negligible, especially for smaller snow thickness. The largest uncertainty of the differential phase is expected for highly wet deep snow layers and larger incidence angles. Indeed, during December 2018, the Central Apennines were mainly influenced by a high-pressure system which caused a positive temperature anomaly and a deficit of precipitations. January 2019 instead was dominated by a low-pressure system which brought low temperatures and a positive precipitation anomaly. During February, again, a high-pressure system caused a positive temperature anomaly and a small deficit of precipitations.

C. Synoptic View of Winter 2018-2019
Thus, monthly anomalies compensate each other and winter 2018-2019 results to be in line with the average of past winters.

D. Estimating Mean Snowpack Density From Snowpack Top Layer Density
The density data collected daily by the Meteomont service unfortunately refers only to the snowpack top layer, thus they are not representative of the mean snowpack density.
To overcome this problem and to be able to use the snowpack top density observations for validation, we used an independent dataset of snowpack vertical profiles, still provided by the Meteomont service, to find a relation between the top and mean snowpack density values.
The new dataset consists of 138 vertical stratigraphies of the snowpack properties done during winters 2018-2019, 2019-2020, and 2020-2021, which include the SD of each layer identified in the snowpack. The stratigraphies are done weekly, and since just a few of them occurred on the same day of SAR acquisitions, we could not use them for validation.
Thus, for each vertical profile, we calculated the mean density and we extracted the density of the top layer. Then we compared top and mean density obtaining a correlation coefficient equal where ρ mean and ρ top are the mean and the top snowpack densities, respectively. Substituting to ρ top the density value coming from the daily Meteomont dataset, we were able to estimate the corresponding ρ mean , and use it to validate the ANN estimations.

ACKNOWLEDGMENT
The authors wish to thank the ESA Research Service and Support (RSS) group from Frascati (Rome, Italy) for the valuable help in improving the SAR processing chain. The European Space Agency and its Network of Resources (NoR) are also acknowledged for sponsoring and providing high-performance computing resources. We are grateful to the Italian Meteomont Avalanche Warning Service (MAWS) (https://meteomont.carabinieri.it/home) for providing daily snow data within the area of interest. The authors wish to commemorate one of the co-authors of this manuscript, Prof. Frank Silvio Marzano, who suddenly passed away on May 2022 after the research work presented herein had already been completed. His support was fundamental in the organization of the study and in the writing of this article. His contribution was-and still is-inspiring in many ways.
Paolo Tuccella received the Laurea degree in physics