Retrieval of Dust Properties From Spectral Snow Reflectance Measurements

Introduction

The spectral reflectance of snow depends on several variables related both to the observation system (e.g., view zenith angle, diffuse or direct radiation etc.) and the physical properties of the surface (Warren, 1982). These latter can be summarized as: snow grain size (Nolin and Dozier, 1993), snow liquid water content (Green et al., 2006), snow density (Dadic et al., 2013) and presence and type of light-absorbing impurities (also referred as light-absorbing particles) (Di Mauro et al., 2015; Dumont et al., 2017; Skiles et al., 2018). Therefore, respective parameters can be derived from snow spectral reflectance measurements obtained both from field spectroscopy, airborne and satellite data, as shown in various studies (Nolin and Dozier, 1993; Green et al., 2006; Negi and Kokhanovsky, 2011; Di Mauro et al., 2015; Seidel et al., 2016). The retrieval of these parameters from remote sensing data is an active field of research (Bohn et al., 2020) and it is useful for assessing the state of the snowpack for hydrological and climatic-related applications. In particular, the development of algorithms for the quantitative retrieval of snow parameters is relevant in the context of current and future space imaging spectroscopy missions.

When snow is mixed with various impurities, its optical properties are modified (Warren and Wiscombe, 1980; Warren, 1982). This process leads to an increase in the absorption of solar radiation and an enhancement of surface melting (Di Mauro et al., 2019). The retrieval of the abundance and optical properties of impurities in snow is still an open issue in the remote sensing of the cryosphere. Furthermore, the timing and intensity of impurities depositions can be directly assimilated in snow hydrology models for runoff prediction (Cluzet et al., 2020).

The determination of snow properties from reflectance measurements has been discussed in a previous study (Kokhanovsky et al., 2018). In this paper, we present a simplified algorithm to retrieve snow grain size, dust mass absorption coefficient and dust mass concentration from spectral diffuse reflectance measurements at three wavelengths located in the visible and near infrared (NIR). The theoretical model is then compared with field spectroscopy data collected in different days from an automated spectrometer in the Alps.

Theory

The spherical albedo r of a semi-infinite plane - parallel snow is defined as

where $W_{refl}$ is the total power scattered into the upward hemisphere from unit area of a snowpack that is being uniformly illuminated by diffuse light from the entire upper hemisphere (say, from overcast sky) and $W_{inc}$ is the total power incident on unit area of the surface (Kokhanovsky et al., 2018). The total power absorbed by the snow uniformly illuminated by diffuse light from the entire upper hemisphere $W_{abs}$ can be found from the following expression:

By definition, it follows that $r=1$ in case there is no absorption of light in snow, which is a good approximation for the clean fresh snow in the visible spectral bands, because of very weak absorption of light by ice in the visible. In near infrared, the absorption of light by ice grains cannot be ignored and the following approximation for the semi-infinite snow spherical albedo can be used (Kokhanovsky and Zege, 2004; Kokhanovsky et al., 2018):

where $l_{tr}=1/\left(1-g\right){\sigma }_{ext}$ is the transport length and $l_{abs}=1/{\sigma }_{abs}$ is the absorption length, $\epsilon =4/\sqrt{3}$, g is the average cosine of scattering angle, ${\sigma }_{ext}$ is the extinction coefficient, ${\sigma }_{abs}$ is the absorption coefficient. Since the concentration of impurities in snow is often low, the value of $l_{tr}$ depends almost solely on the scattering and extinction characteristics of ice grains. The accuracy of Eq. 3 valid for vertically homogeneous plane–parallel semi–infinite snow layers has been studied in Kokhanovsky et al. (2019). An important point is that the transport length for weakly absorbing grains which have sizes much larger than the wavelength λ of the incident light is a very weak function of the wavelength (Kokhanovsky and Zege, 2004). This means that the spectrum of spherical albedo is determined almost solely by the spectrum of the snow absorption coefficient ${\sigma }_{abs}$. This coefficient can be presented in the following way for the case of polluted snow:

where $c_i$ is the volumetric concentration of ice grains in snow, $c_p$ is the volumetric concentration of pollutants in snowpack,

are volumetric absorption coefficients of ice and pollutants, respectively, ${\overline{C}}_{abs,i}$ is the average absorption cross section of ice grains, ${\overline{C}}_{abs,p}$ is the average absorption cross section of impurities in snow, ${\overline{V}}_i$ is the average volume of grains in snow and $\overline{V_p}$ is the average volume of pollutant particles in snow. The impurity absorption coefficient ${\sigma }_{abs}^p=c_p{k}_p$ can be derived from Eqs. 3, 4:

The second term in Eq. 6 can be ignored in the visible. Then it follows:

An important point is that ${\overline{C}}_{abs,i}$ is proportional to the volume ${\overline{V}}_i$ for weakly absorbing ice grains (Kokhanovsky and Zege, 2004):

where ${\gamma }_i\left(\lambda \right)$ is the bulk ice absorption coefficient. The parameter $B$ depends on the shape of ice particles (Kokhanovsky and Macke, 1997; Libois et al., 2014) and real part of ice refractive index, which is a weak function of the wavelength in the visible and near–infrared. If pollutant particles are much smaller as compared to the wavelength of the incident light, then ${\overline{C}}_{abs,p}$ is also proportional to the volume of absorbing particles and it follows:

where ${\gamma }_p\left(\lambda \right)$ is the bulk pollutant (say, soot) absorption coefficient. The parameter D depends on the shape of particles and also on the real part of pollutant refractive index. Therefore, one derives:

where ${\gamma }_i=\frac{4\pi {\chi }_i}{\lambda }, {\gamma }_p=\frac{4\pi {\chi }_p}{\lambda }$, ${\chi }_i$ is the imaginary part of ice refractive index at the wavelength λ, ${\chi }_p$ is the imaginary part of pollutant refractive index at the wavelength λ.

Therefore, it follows for the polluted snow spectral spherical albedo:

or

where $\beta =cD/B,$$c=c_p/c_i$ is the relative pollutant concentration and

is the effective transport length (ETL). One can conclude from Eq. 12 that the soot - polluted snow spherical albedo spectrum in the visible and near infrared is determined by just two constants: $l$ and $\beta ,$ which do not depend on the wavelength λ. Therefore, it is of importance to derive and report these values for different types of snowpacks. Let us assume that pollutants do not influence spectra in the near–infrared. This is a valid assumption at least for mass concentrations of impurity below 1000 ppm (Dang et al., 2015). Then it follows:

where ${\lambda }_{nir}$ is the wavelength in the near–infrared (say, 1020 nm). One can derive for the value of $\beta :$

where ${\lambda }_{vis}$ is the wavelength in the visible (say, 400 nm). We arrive to an important conclusion that information on just two parameters (or measurements at two wavelengths, say 400 and 1020 nm) make it possible to determine the complete spherical albedo spectrum of snow containing impurities, which are much smaller than the wavelength of incident light under assumption that spectra ${\gamma }_i\left(\lambda \right)$ and ${\gamma }_p\left(\lambda \right)$ are known. The parameter $\beta$ is determined by the relative concentration of pollutants c. Therefore, the concentration of c can be derived from $\beta :$

where $\kappa =B/D.$ The value of B depends on the shape of grains and is close to 1.6 for typical snowpacks (Libois et al., 2014). The value of D depends on the type of impurities being close to 1.3 for soot assuming that its refractive index is equal 1.95-0.79i across solar spectrum. One can see that $\kappa$ is close to unity and the parameter $\beta$ can be used to estimate the relative volumetric concentration of soot particles in snow.

Let us consider the physical meaning of the second spectrally invariant parameter $\text{of snow} l$. It follows from Eq. 13:

where we have accounted for the fact that the snow extinction coefficient can be derived from the following geometrical optics equation (Kokhanovsky and Zege, 2004):

where d is the effective ice grain size and ${\overline{ S}}_i$ is the average geometrical cross section of grains in the direction perpendicular to the incident light.

It follows that the value of $l$ is proportional to the effective grain size. Therefore, we can derive the effective grain size from the value of $l$. Namely, it follows:

where

Taking into account that the ratio B/(1-g) is close to 9 (Kokhanovsky et 2019), we derive: $\xi =1/16.$ It should be pointed out that although the snow effective grain size d and concentration of soot c can be derived from the spectral invariants $\beta ,l$ assuming the calibration constants $\left(\kappa , \xi \right),$ the invariants $(\beta ,l)$ derived from the spectra are important characteristics of snow to be reported in the output of snow retrieval algorithms. These parameters can be derived from the measurements itself not assuming a particular shape of ice grains, which is difficult to determine in practice. In case of clean snow, the value of $\beta$ is equal to zero and the spherical albedo spectrum is characterized by just one parameter $l.$ This explains why the model of spherical ice grains can be used to describe the spherical albedo spectrum of clean snow so accurately (Wiscombe and Warren, 1980). Indeed, the effective transport length is the single snow parameter, which can be derived from the measurements. Assuming the wrong shape of particles, one uses wrong constant $\xi$ and derives the biased value of d–although the ETL $d/\xi$, which determines the snow spectral spherical albedo is correct.

We have assumed that impurity particles are much smaller than the wavelength of incident light. This is often not the case. Let us assume now that the volumetric absorption coefficient of impurities can be presented as (Kokhanovsky et al., 2018):

where $\alpha$ is the absorption Angström parameter and $k_0$$=k_p\left({\lambda }_0\right)$. We shall assume that ${\lambda }_0=1\mu m$. Then it follows instead of Eq. 11:

or

where

and $l$ is determined by Eq. 13 (Dalzell and Sarofim, 1969). The parameters $l$ and $\beta$ have the same meaning as discussed above except the value of D must be substituted by$k_0$. The parameter $\alpha$ describes the spectral properties of the absorption coefficient of pollutants in the visible. It is equal to one for soot. However, generally speaking, it can take various values depending on the type of pollutants, their size distributions and spectral optical constants. Therefore, we conclude that the snow spectral spherical albedo is determined by just three spectral invariants: $\alpha ,\beta ,l$. Therefore, each snow sample must be characterized by these constants. This makes it possible to determine the snow spectral spherical albedo at any wavelength in the visible and near infrared. Also, we can derive the snow broadband albedo, snow plane albedo $r_p$ and snow reflection function $R.$ In particular, it follows (Kokhanovsky and Zege, 2004; Kokhanovsky et al., 2018; Kokhanovsky et al., 2019):

where $R_0$ is the reflectance of a nonabsorbing semi–infinite snow layer, $\psi =u\left({\mu }_0\right)u\left(\mu \right)/R_0$, $u\left({\mu }_0\right)$ is the so–called escape function of a semi-infinite snow layer (Kokhanovsky et al., 2019), ${\mu }_0$ is the cosine of the solar zenith angle (SZA), $\mu$ is the viewing zenith angle. The values of $R_0$ and u can be derived using various parameterizations (Kokhanovsky et al., 2019) or from the measurements itself (Kokhanovsky et al., 2018). The equations presented above can be therefore used in context of multispectral sensors and can be easily optimized for hyperspectral space sensors. In this work we analyze the ground spectral measurements of the snow plane albedo, which is the spherical albedo power the escape function (see Eq. 26). We present simple approximate equation for the escape function in Appendix 1.

Determination of Polluted Snow Properties From Spectral Reflectance Measurements

The theory presented above has been used to process the spectral reflectance measurements collected in the field at the Torgnon experimental site (Valle d’Aosta, Italy, 45.845012N, 7.578517E) on May 16, 17, 18, 2018, in the presence of dust - loaded snow. Mineral dust deposited on seasonal snow at this site has probably a Saharan origin (Di Mauro et al., 2019), as reported by other studies in the Alps (Dumont et al., 2017).

In Figure 1, we show an example of a picture of the site acquired on the 17th of May (2018) from a webcam (namely Phenocam Network) installed at the site (see Di Mauro et al., 2019 for further details).

FIGURE 1. The photo of the Torgnon site in Italian Alps (45.845012N, 7.578517E) taken on May 17, 2018. Dust is visible on the surface of snow.

This experimental site is equipped with several instruments for snow modeling and snow hydrology (Colombo et al., 2019; Di Mauro et al., 2019). The site is seasonally covered by snow from October to May and it is located at the altitude of 2,160 m in the Italian European Alps. Since November 2017, a RoX spectrometer (JB Hyperspectral Devices, Düsseldorf, Germany) is installed at the site, and it continuously acquires spectral data each 3 min. The spectrometer collects incident and reflected radiation in the wavelengths between 400 and 900 nm. This makes it possible to derive the spectral plane (or spherical for the overcast sky) albedo. The instrument features a full width at half maximum of 1.5 nm, and a spectral sampling of 0.65 nm (JB Hyperspectral devices). Albedo is calculated as the ratio between the reflected and incident fluxes. Following previous studies (Picard et al., 2016a; Dumont et al., 2017), the automatic system is equipped with a cosine diffusor both for the upwelling and downwelling fluxes in order to measure the Bi-Hemispherical Reflectance factor (BHR) (here referred as spectral plane albedo). For this study, 30 measurements collected at the minimum of the solar zenith angle have been averaged. In order to remove noise from data, spectra were smoothed with a Savitzky-Golay filter by using a third order polynomial and a kernel of 99 data points. Albedo spectra were corrected for a constant slope (10 degrees) and aspect (230 degrees) of the snow surface, as proposed in Weiser et al. (2016). Slope and aspect values were calculated from a digital surface model obtained from a drone survey of the snow-covered area. The measurements above 865 nm appear still noisy and, therefore, are not present in the subsequent analysis.

The field data provides the plane albedo. The spherical albedo is estimated from the plane albedo using Eq. 26. Namely, it follows:

The parameter α has been derived from Eq. 24 using BHR measurements at two wavelengths in the visible (${\lambda }_1=410nm, {\lambda }_2=500nm$), which were selected because available in many space sensors:

where $z=\ln r\left({\lambda }_2\right)/\ln r\left({\lambda }_1\right)$. It follows for the product $b=\beta l$ (see Eq. 24) under assumption that absorption of light by grains in the visible can be neglected as compared to light absorption by dust grains:

The parameter $l$ is derived from Eq. 24 at the wavelength ${\lambda }_3=865nm$:

Then it follows:

The relative volumetric concentration of dust particles can be estimated from Eq. 25:

The relative mass concentration $c_m$ is related to the value of c by the following equation:

where ${\rho }_d=2.65g/c{m}^3$ is the density of dust and ${\rho }_i=0.917g/c{m}^3$ is the density of ice. Therefore, it follows:

One can see that one needs the value of $k_0$ to derive the relative mass concentration $c_m.$ The value of $k_0$ can be estimated from the derived value of α. In particular, we have used the following relationship between $k_0$ and α derived using the electromagnetic scattering theory (see Supplementary Appendix 2):

where $b_0=10.916, b_1=-2.0831, b_2=0.5441$ and $k_0$ is in $m{m}^{-1}.$ The derived dust mass absorption coefficient (MAC), that is ${\sigma }_{abs}^m=k_0/{\rho }_d{\left(\lambda /{\lambda }_0\right)}^{-\alpha }$ for three days is presented in Figure 2, where we also show the MAC measured experimentally (Caponi et al., 2017) in a chamber for the dust aerosol from Sahara (Algeria, particles with diameters smaller than 10.6 µm). The close agreement of experimental and derived (for May 18, 2018) dust MAC has been found. Also absorption Angström exponents for the case May 18, 2018 and that measured in the chamber (2.5, see Caponi et al. (2017)) almost coincide. We have also derived the effective diameter of dust grains $d_{ef}$ (see Supplementary Appendix 2).

FIGURE 2. The derived dust mass absorption coefficients for the three days. Symbols represent experimental measurement of the mass absorption coefficient in the CESAM simulation chamber (Caponi et al., 2017).

The solar zenith angle at the time of measurements needed for the determination of spherical plane albedo and results of calculation of various snow parameters for May 16, 17, 18, 2018 at the Torgnon site are given in Table 1. They are based on theory presented above and also on the plane albedo measurements performed at three wavelengths: 410, 500, and 865 nm.The plane albedo was modelled using the following simple analytical equation derived from Eqs 24, 26:

with the parameters shown in Table 1. The function $u\left({\mu }_0\right)$ has been calculated as described in Appendix 1; Eq. 36 can be used to estimate spectral plane albedo at any other solar zenith angle, if the measurement at a particular SZA known (under assumption that snow properties have not been changed during the corresponding time interval).

TABLE 1. The parameters of dust–loaded snow derived at three dates at the Torgnon site.

We show the spectral albedo measured on site and the results of calculations using Eq. 36 with the parameters (α, β, $l$) specified in Table 1 for each date in Figure 3. The spectral imaginary part of ice refractive index has been taken according (Warren and Brand, 2008) with the updates for the visible part of electromagnetic spectrum proposed in Picard et al., 2016b. It follows that Eq. 36 can indeed be used to represent the spectral plane albedo of dust–loaded snow. This opens an avenue for the determination of various snow intrinsic characteristics based on analytical Eq. 36 as discussed above. Unfortunately, concentration and size distribution of dust at the site are not available for May 2018. In a previous work (Di Mauro et al., 2019), we showed that dust dominates the impurity content in snow during the melting season. We can reasonably state that the derived concentration and size of dust particles are close to those reported in previous work (Di Mauro et al., 2015; Dumont et al., 2017; Di Mauro et al., 2019) dealing with same type of impurities. The derived values of absorption Angström parameter are close to those reported in Caponi et al. (2017).

FIGURE 3. The spectral snow albedo: measurements–solid lines, calculations using Eq. 36–dashed lines. Dates refer to the acquisitions of the albedo spectra. Arrows show the wavelengths used in the inversion procedure (410, 500, and 865 nm).

Regarding the data collected on 18th May, we noticed a mismatch between observed and modeled spectra between 650 and 750 nm. Since this spectral range is particularly sensitive to chlorophyll concentration, the mismatch can be due to other impurities not accounted for in our study, such as cryospheric algae (Di Mauro et al., 2020). Further analysis is needed to decouple the effect of organic and inorganic impurities in snow at this experimental site in the European Alps.

During the three days considered in this study, we observe a steady increase of dust load (11.7–106.9 ppm) in surface snow. This is explained by the resurfacing of impurity layers during the melting season (Di Mauro et al., 2019; Skiles and Painter, 2019; Tuzet et al., 2020), a well-known process in the Alpine region. The reinforcing albedo feedback induced by the concentration of impurities in surface snow strongly accelerates the melting and represent an important application for remote sensing and snow hydrology studies.

The possibility of directly estimating dust concentration and snow grain size represents a breakthrough in the remote sensing of the cryosphere. In fact, this will allow to map the spatial and temporal variability of impurities in snow at global scale and to evaluate the role of impurities on the melting processes of snow and ice.

Conclusion

In this paper, we present a technique for the simultaneous determination of the size of ice grains and size of dust grains and also mass concentration of pollutants in snow. The technique to determine the mass absorption coefficient of impurity particles deposited in snow is also proposed. The technique is valid for the case of plane–parallel semi-infinite vertically homogeneous snow layers illuminated by a plane parallel direct solar light (or by the diffuse light from the overcast sky). The proposed method has been tested on field spectroscopy data and can be extended to account multiple pollutants in snowpack and to remote sensing data with different spectral resolution. Also, one can account for the fact that in reality the measured albedo is a mixture of white sky and blue sky albedo determined by the diffuse to direct light ratio. Although no validation data are available for this study, the modeling results obtained with the proposed model are comparable with previous literature focusing on the impact of dust in snow. The validation of the proposed technique will be a subject of our future study. This will lead to the simple and robust technique needed for the quantitative estimations of various polluted snow intrinsic characteristics from multispectral and hyperspectral data. The technique could be also applied to the analysis of snow reflectance spectra $R\left(\lambda \right)$ (see Eq. 26). The assumption (see Eq. 14) related to the strength of absorption by pollutants (in the near infrared) can be dropped and other ways for the estimation of the triplet (α, β, $l$) from Eq. 36 can be used (e.g., optimal estimation technique based on all available spectral channels).

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

AK: theory, conceptualization. BDM, RC, and RG: experimental measurements. AK, BDM, RC, and RG: writing of the paper. All authors reviewed and approved the final version.

Conflict of Interest

Author AK was employed by company Vitrociset Belgium SPRL, A Leonardo Company.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The data shown in Figure 2 by symbols have been supplied by Adriana Rocha-Lima. This is acknowledged with many thanks. AK acknowledges the support of the ESRIN/ESA contract 4000125043-ESA/AO/1-9101/17/I-NB EO SCIENCE FOR SOCIETY. BDM, RG, and RC acknowledge the support of the MUSICA project (Multiband Ultrawide SpectroImager for Cryosphere Analysis, contract: 2019-NAZ-0045) funded by the Italian Space Agency (ASI). BD, RG, and RC would like also to thank Tommaso Julitta (JB Hyperspectral Devices) for the help in pre-processing the Rox data, and Edoardo Cremonese and the personnel of the Environmental Protection Agency of Aosta Valley (ARPA-VDA) for maintaining the experimental site in Torgnon (AO, Italy). The hyperspectral data was processed from raw digital numbers to calibrated radiances using the FieldSpectroscopyDP (https://github.com/tommasojulitta/FieldSpectroscopyDP) and FieldSpectroscopyCC (https://github.com/tommasojulitta/FieldSpectroscopyCC) packages in R (R Core Team, 2017).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2021.644551/full#supplementary-material

Appendix 1

Escape Function

We shall assume that the snow phase function is described by the Henyey – Greenstein phase function with the asymmetry parameter g=0.75 (Henyey and Greenstein, 1941; Kokhanovsky et al., 2019). The escape function has been numerically calculated for this specific case in Yanovitskij (1997). The corresponding dependence is shown in Figure A1.1 together with approximation derived by us from the fit of numerical values:

FIGURE A1.1. The dependence of the escape function on the cosine of the solar zenith angle for the Henyey – Greenstein phase function at g=0.75 calculated using exact radiative transfer theory (symbols), and Eq. (A1.1) (solid line) and Eq. (A1.2) (dashed line).

It follows from Fig. A1.1 that a linear approximation for the escape function (Sobolev, 1972; Kokhanovsky et al., 2019):

can be also used (at μ₀>0.2).

References