On the Relationship Between Stickiness in DMRT Theory and Physical Parameters of Snowpack: Theoretical Formulation and Experimental Validation With SNOWPACK Snow Model and X-Band SAR Data

This study aims at relating the stickiness parameter (<inline-formula> <tex-math notation="LaTeX">$\tau$ </tex-math></inline-formula>) of the dense media radiative transfer theory in quasi-crystalline approximation of Mie scattering of densely packed sticky spheres (DMRT-QMS), to the physical parameters of the layered snowpack. A relationship has been derived to express <inline-formula> <tex-math notation="LaTeX">$\tau $ </tex-math></inline-formula>, which modulates the attractive contact force between ice spheres, as a function of ice volume fraction (<inline-formula> <tex-math notation="LaTeX">$\phi $ </tex-math></inline-formula>) and coordination number (<inline-formula> <tex-math notation="LaTeX">$n_{c}$ </tex-math></inline-formula>). Since <inline-formula> <tex-math notation="LaTeX">$\tau $ </tex-math></inline-formula> is not a measurable parameter, this is a step forward with respect to what is commonly made in the literature, where <inline-formula> <tex-math notation="LaTeX">$\tau $ </tex-math></inline-formula> is assumed as an arbitrary parameter, generally ranging between 0.1 and 0.3, to fit simulated backscattering data with those measured. As a first validation, DMRT-QMS was integrated with the SNOWPACK model to simulate backscattering at X-band (9.6 GHz) driven by nivo-meteorological data acquired on a test area located in Monti Alti di Ornella, Italy. The simulations were compared with Synthetic Aperture Radar COSMO-SkyMed (CSK) satellite observations. The results show a significant agreement (<inline-formula> <tex-math notation="LaTeX">$R^{2} =0.68$ </tex-math></inline-formula>), although for a limited dataset of eight points in a unique winter season.

values of albedo with respect to natural bodies, therefore 32 reflecting high quantity of solar energy and lowering the 33 surface temperature [4]. Consequently, a decrease in snowpack 34 extension and duration leads to an increase of local warming 35 and accelerate the snow melting and water runoff processes, 36 causing avalanches and floods, in particular flash floods. 37 In fact, several studies predicted that in several mountains 38 of the European Alps, snow water equivalent (SWE) may 39 be reduced significantly by 2100 and that they may become 40 totally snow-free during the summer [5], [6]. This could lead 41 to a point of no return in the management of water resources, 42 energy balance, and winter tourism. Snow is important also in 43 permafrost behavior, due to its thermal insulation characteris-44 tics. As a result of all these aspects, continuous monitoring of 45 snowpack can be crucial and represent an important tool for 46 water resources management [7]. 47 From the physical point of view, snow can be considered 48 as a mixture of air, ice grains, and, if present, liquid water. 49 Since the ice grains occupy between 10% and 40% of the 50 snow volume, snow can be studied as a dense medium with a 51 high fractional volume of scatterers. Consequently, variations 52 in the scatterers inside snowpack cause significant changes 53 of the electromagnetic response. In the microwave range, 54 a variety of radiative transfer model of snowpack has been 55 proposed, such as MEMLS [8] and HUT [9], to simulate 56 the scattering and emission of a layered snowpack. Among 57 them, dense media radiative transfer theory in quasi-crystalline 58 approximation of Mie scattering of densely packed sticky 59 spheres (DMRT-QMS) [10], [11], [12], [13] considering the 60 interactions of the electromagnetic waves scattered from the 61 particles. 62 It is important to parameterize snow microstructure in a 63 radiative transfer model of snowpack in microwave frequen-64 cies. In DMRT-QMS snow representation, grains are con-65 sidered spherical particles of ice. The particles' distribution 66 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ is determined by a free dimensionless parameter called ferent shape of grains present in snow layers, instead of using 117 the real single grain radius a. This parameter, which cannot 118 be measured empirically [33], can be analytically expressed 119 as a function of snow density (ρ) and specific surface area 120 (SSA) with R 0 = 3 × 10 3 /(ρ·SSA) [34]. A scaling factor 121 was sometimes used in order to tune the grains' dimension 122 to obtain a reasonable agreement between experimental and 123 theoretical data. This value is not constant and depends on 124 the considered model and the snow characteristics. The same 125 results can be obtained using grain size in alternative to the 126 equivalent optical grain radius, using a new scaling factor. 127 In [35], it was described an interesting relationship between 128 the grain size scaling factor and τ as ρ varies, in order to 129 obtain the same simulation results.

130
The surface roughness of snow was deemed macroscop-131 ically negligible, especially at X-band (9.6 GHz), since 132 previous studies showed a small sensitivity of σ 0 with respect 133 to surface roughness at this frequency [36].This assumption 134 is no more valid if the testing area extends over an area of a 135 few meters, so-called small-scale area; in [37], a snow albedo 136 variation between 1% and 3% was observed in midwinter, 137 reaching 10% during the melting season due to local snow 138 roughness. Likewise, in [38], an experiment carried out on 139 a dry snowpack of about 2 m demonstrated that surface 140 roughness and soil moisture under the snowpack cannot be 141 negligible in backscattering computation. The study was car-142 ried out at 5.3 and 13 GHz from bistatic radar. The soil surface 143 roughness contribution under the snow was described analyti-144 cally and experimentally in a recent study using the statistical 145 S-matrix wave propagation in spectral domain (SSWaP-SD) 146 approach [39].

147
This work was devoted at assessing the relationship between 148 τ and the physical parameters of the snowpack, by combining 149 DMRT-QMS and SNOWPACK model simulations, COSMO-150 SkyMed (CSK) X-band synthetic aperture radar (SAR) acqui-151 sitions, and in situ data. Here, SNOWPACK is one of the 152 multilayer snow hydrology models developed by the Swiss 153 Federal Institute for Snow and Avalanche Research, SLF, 154 Davos, Switzerland, for details see Section II-D. First, the τ 155 dependence on physical snow parameters is analyzed based on 156 SNOWPACK simulations, and the optimized τ is estimated by 157 fitting the backscattering simulated by coupling SNOWPACK 158 and DMRT-QMS (σ 0 SP/DMRT ) with the CSK measurements 159 (σ 0 CSK ). The analysis pointed out a significant correlation 160 between τ , ς, and n c . Then, the theoretical relationship 161 between τ , ice volume fraction (ς), and coordination number 162 (n c ) valid in the SHS model is established. The relation-163 ship, which can be considered the main result of this work, 164 is validated in the Monti Alti di Ornella test area, in the 165 eastern Italian Alps, based on SNOWPACK simulations and 166 CSK measurements collected on the area. The input nivo-167 meteorological data for SNOWPACK were provided by the 168 automatic weather station (AWS) installed in the area.

169
This article is organized as follows. Test area, models, 170 and satellite data are described in Section II. The empirical 171 relationship between optimized τ and physical parameters is 172 analyzed in Section III. The theoretical relationship between 173 τ , ς, and n c is established in Section IV and the preliminary 174 validation is described and discussed in Section IV-A.    Table I 216 except ISWR (see Table I). Some data are not fully complete 217 due to temporal gaps. These gaps are usually around 2-3 h. [45]. The fundamental parameters are listed in Table I, together 237 with some ancillary information as ground slope angle [ • ], and 238 date and time combined in the ISO 8601 format. The slope 239 angle was 0 • with an approximation within 5 • after considering 240 the local DEM. Some ancillary information is also necessary 241 and consists of canopies information (if present), snow and 242 soil initial conditions, as well as snow albedo value.

243
If some input parameters are missing, new data can be 244 generated. In one case, the goal is to create new parameters as 245 functions of other variables. In another situation, a statistical 246 temporal interpolation is used when few isolated data failures 247 are present (e.g., a not working sensor or a meteorological 248 parameter that was not measured) [46]. In our case, ISWR 249 was not available and therefore was generated from TA 250 (see Table I) and RH (see Table I) corrected by cloudiness, 251 if available [47].  The scheme in Fig. 2 provides a depiction of the inter-303 action among experimental data, model parameters, models, 304 and model simulations. The AWS data were processed by 305 SNOWPACK. Some of the outputs were used as inputs for 306 DMRT-QMS with the addition of τ , mv, and rms. The simu-307 lated backscatter was compared to the measured one by CSK. 308

310
Some of the nivo-meteorological parameters measured by 311 the AWS were used as inputs to SNOWPACK to estimate 312 the characteristics of the snowpack [66] in the winter season 313 2014-2015. The input data were measured every hour and the 314 output data were generated every 15 min. The simulations of 315 the snow type's distribution inside the snowpack, the presence 316 of liquid water content, SWE, and snow ground temperature 317 (TSG) temporal evolution are shown in Fig. 3.

318
The winter season was characterized by some very 319 mild days, i.e., 18/12/2014, 10/01/2015, and 11/02/2015 320 (DD/MM/YYYY, hereafter), in which the energy input in 321 the snowpack was difficult to simulate, also considering the 322 northern exposure of the site. This nivo-meteorological con-323 dition implies that small variations in the control parameters 324 of the SNOWPACK model can determine the possibility of 325 even different snowpack simulations. The one proposed in this 326 article, although characterized by a snowpack a little colder 327 than reality, finds a good match with the snowpack profiles 328 carried out on the ground by the Veneto Regional Agency for 329 Environmental Protection and Prevention (ARPAV), as shown 330 in Fig. 3(d) for HS and SWE. The latter was measured 331 in two different modes: "Stratigraphy SWE" was performed 332 as explained in [67] and "Coring SWE" was performed as 333 explained in [68]. About grain size, some papers have shown 334 a good correspondence between SNOWPACK results and 335 in situ observations, as, for example, in [69]. In our case, 336 the comparison between the in situ measurements and the 337  observes the maximum size of the grains. These differences 345 have also been illustrated in [70]. Therefore, the values of 346 the grain sizes are not comparable, but, given the dimensions 347 observed in the field, those of the model can be considered 348 real and valid.  Table III. 355 We deduced that the soil under the snowpack was not frozen 356 by the TSG values (see Fig. 3) and the knowledge of historical 357 in situ measurements [71]. TSG ranged between 0 • C and 358 −1 • C; only in a few cases, it dropped below −1 • C for 359 very short time, probably due to measurement errors. This 360 is an important assumption because, otherwise, the use of 361 the Oh model would not give correct results. If the soil is 362 completely frozen, i.e., there is no liquid water in it, it should 363 be considered as dry, and therefore, the mv values should be 364 assumed close to zero.

365
A strong variability of σ 0 SP/DMRT as a function of τ can 366 be observed in Fig. 4(a). This diagram demonstrates the 367 importance of a careful choice of τ . A significant vari-368 ability of σ 0 SP/DMRT versus rms and mv is also evident in 369  The SNOWPACK simulations on 10/04/2015 and 378 12/05/2015 show that liquid water was present in the 379 snowpack (see Fig. 5). In April, the snow melting started 380 from the cover top, while in mid-May, the snow melting 381 affected the whole snowpack. In this case, DMRT-QMS 382 cannot predict reliable results. Due to this reason, these 383 last two points were not considered in our simulations. 384 Furthermore, on 10/04/2015, there is noτ that reproduces 385 σ 0 CSK , unless a lower mv value is assumed.

386
Since the values of the snowpack parameters are not con-387 stant but depend on the number of layers, the following 388 implements the following function of ρ: In this section, we derive an analytic relationship among 413 ς, ς, and n c holding in the SHS model. In general, the 414 problem of knowing the position of the scatterers to which 415 apply the scattering rules of DMRT has been addressed 416 in two ways. In the first, which we do not adopt in this 417 work, polydisperse scatterers are considered, i.e., the spheres' 418 diameters are randomly distributed according to a specific 419 probability distribution. In the second, according to the SHS 420 model, monodisperse scatterers are considered. The scatterers 421 are subject to the action of an attractive contact interaction 422 modulated by τ . The effect of τ is the tendency of the 423 scatterers to form clusters, the larger they are, the lower τ is. 424 The case of random distribution is recovered for τ → +∞. 425 In order to compute the scattering from many particles, 426 the classical theory assumes independent scattering and the 427 scattering intensity is set equal to the sum of scattering 428 intensities from each particle [72], [73]. The classical approach 429 ignores the coherent wave interaction among the particles. 430 To consider the collective scattering effects, QCA can be used 431 where k B is Boltzmann's constant and T is the temperature.

444
The potential consists in a hard-core repulsion, which prevents The region above τ perc (ς) is defined as the nonpercolating 465 phase, while the region under τ perc (ς) is defined as the 466 percolating phase. In the percolating phase, particle clusters 467 with average size tending to infinity (system-spanning clusters) 468 exist. This is reflected by the behavior of the average n c of 469 SHS, which is well described by PY, and is the average number 470 of spheres that is in contact with a given one. For given ς and 471 ς, n c is given by [75] 472 n c = 2ςε(ς, τ ) where ε(ς, τ ) is the smallest solution of the quadratic equation 474 A contour plot of the average n c is added in Fig. 8. This 478 indicates that the locus of the percolation points τ perc (ς) in 479 the (ς, τ )-plane, originally identified with the divergence of 480 the mean cluster size, corresponds to n c = 2, as a mean-field 481 condition of minimum connectivity.

482
Starting from (6) and (7), we can derive ς as a function of 483 n c and ς with the condition 486 n c < 4(3ς(1 + ς/2)) 1/2 (1 − ς) (10) 487 see Fig. 7 for a plot. The last condition is a consequence of the 488 fact that only the smallest solution of (7) is acceptable, and it 489 is obtained by setting equal to zero the partial derivative with 490 respect to n c in (9).

492
Following the suggestion given by Löwe and Picard [26], 493 we can identify n c and ς of the SHS model with those derived 494 from SNOWPACK simulations. In this way, due to (9), τ is no 495 longer a free parameter of DMRT-QMS, but a consequence of 496 the snowpack microstructure simulated by SNOWPACK. This 497 model uses four primary independent parameters to describe 498 the microstructure of the snowpack: sphericity, dendricity, 499 grain size, and bond size [42]. All other microstructure para-500 meters of SNOWPACK can be derived from other parameters. 501 In particular, the coordination number n c is derived from snow 502 density ρ by means of (2) [42]. This is the average number 503 of bonds per grain and it determines the interconnectivity of 504 the ice matrix.

505
Before considering the pair (ς, n c ) generated by SNOW-506 PACK as valid arguments for determining ς, it is necessary 507 to verify that n c is actually the smaller of the two values 508 associated with the same τ -value. This is equivalent to fulfill 509 (10). Another condition to be satisfied is (8) with ε replaced 510 by (6). Fig. 7 shows that both conditions are fulfilled for two 511 examples of SNOWPACK data.

512
The identification between n c generated by SNOWPACK 513 with that one of SHS does not guarantee that the snow-  Table V shows an example of input data for DMRT-QMS 536 with ς calculated according to (9). We note that τ mainly  In Section III, we showed that, using the Oh model inte-   (9) are not reasonable for the soil type of the test area. The 566 simulations with the best agreement with σ 0 CSK were for 567 mv = 0.1-0.15 cm 3 /cm 3 . the average value of data determined experimentally by those 588 authors. 589 Leppänen et al. [69] compared the predictions of SNOW-590 PACK with measurements in a taiga snowpack. The results 591 show that SNOWPACK tends, in some cases, to underestimate 592 the grain size G by a factor ranging between 1 and 2. 593 Considering that the larger the diameter of the scatterers D 594 of DMRT-QMS, the higher the simulated backscattering, and 595 the results of Section III on the dependence of backscattering 596 on mv, we can deduce that an underestimation of G can lead 597 to an overestimation of the mv values assumed in this work. 598 However, we cannot verify whether the results in [69] can also 599 be applied to the snow type of our test area located in Monti 600 Alti di Ornella [80]. suggestions.

647
of the Mediterranean region (OPTIMED-WATER) in the frame of FP7 of the study of avalanche dynamics. At the same time, 1100 he began training in the field of avalanche forecast-1101 ing. He is currently a Senior Avalanche Forecaster at 1102 the Arabba Avalanche Center and AINEVA, Trento, 1103 Italy. In 1996, it was the first participation in the scientific expedition to 1104 Antarctica in the context of the PNRA, followed by the 1998 expedition 1105 again in the field of international projects related to remote sensing. In 2005, 1106 the third scientific expedition was on the Antarctic plateau for the study of 1107 mega snow dunes. In Artico, at the Italian base Dirigibile Italia, Ny-Ålesund 1108 (Svalbard), Italy, he continued his scientific research in the field of snowpack 1109 classification with multispectral signatures, participating in the creation of the 1110 SISpec and Snowcrystals.it database. Since 2004, he has been involved in the 1111 operational development of the SNOWPACK model of the snowpack with 1112 the SLF, Davos, Switzerland, and AlpSolut s.r.l., Livigno, Italy. Since 2015, 1113 he has been dealing with snow water equivalent of estimating water resources 1114 for large basins (PO, Piave, Brenta) through the use of MODIS images and 1115 ground monitoring networks. He has been a lecturer in snow, avalanche, and 1116 snow slope stability courses in the courses of AINEVA, of the Italian Alpine 1117 Club and since 2002, in professional training courses for ski instructors and 1118 mountain guides. Since 2000, he has been a member of the Editorial Board of 1119 "Neve e Valanghe" magazine. From 2005 to 2013, he was the representative of 1120 Italy in the international meetings of EAWS. He is the author or the coauthor 1121 of 50 articles published in peer-reviewed international journals, 150 articles 1122 in Italian trade journals, and over 700 citations.