The two-layered radiative transfer model for snow reflectance and its application to remote sensing of the Antarctic snow surface from space

The two-LAyered snow Radiative Transfer (LART) model has been proposed for snow remote sensing applications. It is based on analytical approximations of the radiative transfer theory. The geometrical optics approximation has been used to derive the local snow optical parameters, such as the probability of photon absorption by ice grains and the average cosine of single light scattering in a given direction in a snowpack. The application of the model to the selected area in Antarctica has shown that the technique is capable of retrieving the snow grain size both in the upper and lower snow layers, with grains larger in the lower snow layer as one might expect due to the metamorphism processes. Such a conclusion is confirmed by ground measurements of the vertical snow grain size variability in Antarctica.

1 Introduction

It has been shown in the past that multispectral measurements can be used to retrieve the vertical variation of cloud and snow properties (Li et al., 2011; Kokhanovsky and Rozanov, 2012). In this paper, we propose a fast and robust technique for the determination of the vertical variation of ice grain size in snow using hyperspectral satellite measurements, which opens the way for the large-scale assessment of the vertical inhomogeneity of various snow surfaces. The first retrievals of snow vertical inhomogeneity using EnMAP hyperspectral measurements of top-of-atmosphere (TOA) radiance are performed. EnMAP (Guanter et al., 2015; Storch et al., 2023) is the satellite mission for the assessment of various environmental parameters based on hyperspectral reflected solar light measurements in the spectral range of 410–2,450 nm with a spatial resolution of 30 m. The solution to the inverse radiative transfer problem is based on the two-LAyered snow Radiative Transfer (LART) model developed by us for the simulation of both bottom-of-atmosphere (BOA) and top-of-atmosphere spectral reflectance for the cases of underlying snow surfaces. The unique feature of this model is that it makes it possible to model the snow reflectance using simple analytical relationships between the snow nadir reflectance and ice grain sizes, which facilitates the solution of the inverse problem. Only case 1 snow is considered (dry and clean snow). However, the model can be easily generalized to the case of polluted snow. Other limitations of the forward model used in the retrieval of snow properties are listed as follows:

- the semi-infinite snow layer with no contribution from the underlying surface (1);

- the plane-parallel snow layer (the absence of 3D effects, sastrugi, surface inclination, and horizontal variability of snow properties) (2);

- the close-packed effects are ignored (3);

- the shape of particles is assumed ahead of retrievals (Koch fractals of the second generation) (4);

- only two snow layers are assumed, with the smaller ice grains at the top (5);

- the size distribution of ice crystals is not accounted for (6);

- the atmospheric aerosol optical thickness is small (7);

- the clear sky atmosphere without cloudy media and fogs between the instrument and snow surface under observation is assumed (8).

Some of the assumptions given above are not critical for the proposed theory (e.g., 1, 5, and 6) and can be removed if needed.

The remainder of this paper is structured as follows: in Section 2, we introduce the radiative transfer forward model and discuss the dependence of local and global snow characteristics on the snow microstructure and the optical thickness of the upper snow layer. In Section 3, we propose a semi-analytical approach for the inverse problem solution. In Section 4, we discuss the application of the developed technique for snow remote sensing using the EnMAP L1C radiance data.

2 Forward model

2.1 Theory

In this work, we shall use the shortwave infrared (SWIR) EnMAP channels, where ice absorbs incoming solar radiation and atmospheric contribution (outside gaseous absorption bands) can be neglected for a clean Antarctic atmosphere. Therefore, we assume that the top-of-atmosphere reflectance in the selected channels coincides with the snow reflectance. It is supposed that the snow is composed of a semi-infinite layer underneath and a relatively thin layer of fresh snow at the top. Therefore, the snow reflectance at the cosine of the solar zenith angle $\xi$, the cosine of the viewing zenith angle $\eta$, and the relative azimuthal angle $\varphi$ are modeled as (Sobolev, 1975; Liou, 2002; Kokhanovsky, 2006)

$R\left(\xi ,\eta ,\varphi \right)=R_1\left(\xi ,\eta ,\varphi ,\tau \right)+\frac{t_1^{\downarrow }\left(\xi ,\tau \right)t_1^{\uparrow }\left(\eta ,\tau \right)r_{2\infty }}{1-r_1\left(\tau \right)r_{2\infty }},(1)$

where $R_1$ is the reflectivity of the upper layer with the snow optical thickness (SOT) $\tau$ for the case of a black underlying surface, $t_1^{\downarrow }\left(\xi ,\tau \right)$ is the snow transmittance for the upper layer for the light propagation from the snow surface down, $t_1^{\uparrow }\left(\eta ,\tau \right)$ is the same quantity for the light propagation from the second snow layer up to the snow surface, $r_1\left(\tau \right)$ is the spherical albedo of the upper layer, and $r_{2\infty }$ is the spherical albedo of the lower semi-infinite snow layer. It is assumed that the upper snow layer is homogeneous in the vertical direction, and therefore, $\tau =k_{ext}L$, where $k_{ext}$ is the extinction coefficient and L is the geometrical thickness of the upper layer. The lower snow layer is considered as Lambertian reflector.

The value of $r_{2\infty }$ is given by the very accurate van de Hulst approximation (van de Hulst, 1974):

$r_{2\infty }=\frac{\left(1-as\right)\left(1-s\right)}{\left(1+bs\right)},(2)$

where a = 0.139, b = 1.17, and

$s=\sqrt{\frac{\left(1-{\omega }_{02}\right)}{\left(1-g_2{\omega }_{02}\right)}}(3)$

is the similarity parameter, ${\omega }_{02}$ is the single scattering albedo (SSA), and $g_2$ is the average cosine of the single light scattering angle in the lower snow layer. We shall assume that the upper snow layer is optically thick. In the case of the optically thin layer, the horizontal variability of snow optical thickness starts to play a role, and this effect cannot be accounted for in the framework of the 1D radiative transfer theory used in this work. For the case of optically thick layers, there is a possibility to represent $t_1^{\downarrow }\left(\xi ,\tau \right)$ and $t_1^{\uparrow }\left(\eta ,\tau \right)$ in the following form (Sobolev, 1975; Zege et al., 1991; Kokhanovsky, 2006; Kokhanovsky, 2021):

$t_1^{\downarrow }\left(\xi ,\tau \right)=t_1\left(\tau \right)u\left(\xi \right),(4)$

$t_1^{\uparrow }\left(\eta ,\tau \right)=t_1\left(\tau \right)u\left(\eta \right),(5)$

where $u\left(\xi \right)$ is the escape function and $t_1\left(\tau \right)$ is the diffuse transmission coefficient of the upper layer. The reflection function of the upper layer is represented as (Sobolev, 1975; Zege et al., 1991)

$R_1\left(\xi ,\eta ,\varphi ,\tau \right)=R_{1\infty }\left(\xi ,\eta ,\mathit{\varphi }\right)-\Delta \left(\tau \right),(6)$

where $R_{1\infty }\left(\xi ,\eta ,\varphi \right)$ is the reflection function of a semi-infinite layer with local optical properties of the upper layer. In this work, we shall calculate the function $\Delta \left(\tau \right)$ and coefficients $r_1\left(\tau \right)$ and $t_1\left(\tau \right)$ in the framework of the weak absorption approximation (Zege et al., 1991; Kokhanovsky, 2021). Such an approach is possible because snow can be considered a weakly absorbing turbid medium in the visible part of the infrared region of the electromagnetic spectrum (till the wavelength approximately equal to 1300 nm). At larger wavelengths, snow is a strongly absorbing medium, and the parameter $\Delta \left(\tau \right)\ll 1$ and its accurate calculation are not required. In the case of strong absorption, $R_1\left(\xi ,\eta ,\varphi ,\tau \right)\approx R_{1\infty }\left(\xi ,\eta ,\varphi \right)$. Clearly, the second term in Eq. 6 is of importance only for weakly absorbing media; otherwise, this term is close to zero, and its accurate evaluation is not needed. Therefore, we shall use the following approximations for the functions $\Delta \left(\tau \right),$ $r_1\left(\tau \right)$, and $t_1\left(\tau \right)$ (Zege et al., 1991):

$\Delta \left(\tau \right)=t_1\left(\tau \right)\exp \left(-x_1-y_1\right)u\left(\xi \right)u\left(\eta \right),(7)$

$r_1=\frac{sh{x}_1}{sh\left(x_1+y_1\right)},t_1=\frac{sh{y}_1}{sh\left(x_1+y_1\right)},(8)$

where

$x_1={\kappa }_1\tau ,y_1=4q_1{\kappa }_1,{\kappa }_1=\sqrt{3\left(1-{\omega }_{01}\right)\left(1-g_1\right)},q_1=1/3\left(1-g_1\right),(9)$

where ${\omega }_{01}$ is the single scattering albedo and $g_1$ is the average cosine of the single light scattering angle in the upper snow layer. The parameter ${\kappa }_1$ is called the diffusion exponent in the radiative transfer theory. It gives the rate of the light attenuation in the upper layer with the optical thickness (exp (-${\kappa }_1\tau$)) at large values of $\tau$. The value of $y_1$ is proportional to the similarity parameter $s_1$ for the upper layer as ${\omega }_0\to 1$ ($y_1=4s_1/\sqrt{3}$). From Eq. 8, $r_1\to e^{-y_1}$ and $t_1\to 0$ as $\tau \to \infty$ as it should be (Zege et al., 1991). The escape function is derived as follows (Kokhanovsky et al., 2021):

$u\left(\xi \right)=\frac{3}{5}\xi +\frac{1+\sqrt{\xi }}{3}.(10)$

It follows that $R\left(\xi ,\eta ,\varphi \right)\to R_{1\infty }\left(\xi ,\eta ,\varphi \right)$ as $\tau \to \infty$ as it should be. In addition, the same asymptotic solution can be derived as ${\omega }_{02}\to {\omega }_{01},g_2\to g_1$, and

$r_{2\infty }=\exp \left(-y_2\right),(11)$

where $y_2=4q_2{\kappa }_2,{\kappa }_2=\sqrt{3\left(1-{\omega }_{02}\right)\left(1-g_2\right)},\text{and\hspace{0.5em}}{q}_2=1/3\left(1-g_2\right).$ It should be pointed out that approximate Eqs 2, 11 give almost the same results at small values of the probability of photon absorption (PPA) $\beta =1-{\omega }_0$, which is the case for the weak absorption approximation considered here. In particular it follows, from Eq. 11 as $\beta \to 0$ that $r_{2\infty }\to 1-Ds$, where $D=4/\sqrt{3}$. It follows from Eq. 2, in the same case, that $r_{2\infty }\to 1-D^{*}s$, where $D^{*}=2.309$. It should be pointed out that $D\approx D^{*}.$

One can see that the equations given above make it possible to reduce the calculation of the reflectance of a two-layered snow to the calculation of the reflection function of a semi-infinite snow layer $R_{1\infty }\left(\xi ,\eta ,\varphi \right)$ with the parameters of the upper layer. The function $R_{1\infty }\left(\xi ,\eta ,\varphi \right)$ can be calculated and tabulated using various radiative transfer codes (Yanovitskij, 1997; Mishchenko et al., 1999). In this work, we take into account that EnMAP measurements are performed close to nadir directions, and therefore, we can use the following relationship between the values of $R_{\infty }$ in the nadir direction and the value of spherical albedo $r_{1\infty }$:

$R_{1\infty }=a_0+a_1{r}_{1\infty }+a_2{r}_{1\infty }^2.(12)$

This simple equation has been derived by Kokhanovsky et al. (2024) under the assumption that $\beta <1/2$, which is characteristic for snow media in the optical range of the electromagnetic spectrum, where EnMAP performs the measurements. The parameters $a_0,a_1,a_2$ are given in Kokhanovsky et al. (2024) as functions of the solar zenith angle for the case of the assumed single light scattering angular distribution in the snow (Henyey–Greenstein phase function (Henyey and Greenstein, 1941) with the value of the parameter g = 0.75), and $r_{1\infty }$ is derived using Eq. 2 with the similarity parameter for the upper snow layer.

The proposed LART forward model (Eqs 1–10, (Yanovitskij, 1997)) can be used to simulate the spectral nadir reflectance of snow layers. It is import to compare the results of calculations using the proposed simple analytical forward model and the numerical solution of the integro-differential radiative transfer equation (RTE) (Sobolev, 1975; Liou, 2002) for solar light reflectance. The numerical solution of the RTE was found using the PYDOME forward radiative transfer model, a Python implementation of the discrete ordinate method with a matrix exponential (DOM-ME) approach. One expands the RTE into a cosine Fourier series and transforms the light radiance in the angular domain into its discrete counterpart by selecting a discrete set of Gaussian points and weights in the framework of the DOM-ME approach. Therefore, the integral term of the RTE is changed to a simple summation operation, leading to a system of differential equations solvable via the matrix exponential method. To enhance the efficiency of solving the eigenvalue problem, a scaling transformation, as described by Efremenko et al. (2013), is applied. The theoretical foundation of this method is detailed by Efremenko et al. (2017). PYDOME operates in two modes. In the first mode, layer equations for each atmospheric layer are derived and integrated into the system. By applying boundary conditions at the top and bottom of the atmosphere and ensuring radiance continuity at layer boundaries, the system yields the unknown radiance components. In the second mode, the problem is approached through the Riccati matrix equation to determine reflection and transmission matrices (Chuprov et al., 2022). The solution for the multilayered medium is obtained using the matrix operator method, which computes the reflection and transmission matrices for the entire medium. The model incorporates several phase function truncation techniques, including the delta-M (Wiscombe, 1977) and delta-fit (Hu et al., 2000) methods and TMS (Nakajima and Tanaka, 1988) correction. An optional modification for the small-angle approximation of the discrete ordinate method is available to bypass phase function truncation (Afanas’ev et al., 2020). For this study, the first mode was utilized, augmented with delta-M and TMS correction methods. The number of discrete ordinates was increased in powers of 2 until the results from two consecutive iterations differed by less than 0.0001%. A total of 64 discrete ordinates were utilized. The necessary number of azimuthal harmonics was determined based on the Cauchy convergence criterion, which had to be satisfied twice to ensure the accuracy and reliability of the computational results. The dependence of reflectance on the upper layer optical thickness for various values of the SSA in both layers is shown in Figures 1A, B. The calculations have been performed using both the LART forward model and the numerical PYDOME RTE solution. Except for the case of SSA equal to 0.9999, the semi-infinite snow limit is achieved at the upper snow optical thickness in the range of 10–40, depending on the level of light absorption in the upper layer. The cases with a single scattering albedo larger than 0.9999 occur in the visible range (Kokhanovsky, 2021); the spectral range not used in the inverse problem solution is discussed below. Figure 1C shows that the accuracy of our forward model is excellent. The relative errors are smaller than 5% at solar zenith angles larger than 78.5°, which is comparable to the accuracy of the EnMAP measurements. This is of great importance for the simplification of both forward and inverse radiative transfer problems in layered snow optics.

Figure 1

Figure 1. (A) Dependence of the reflectance of a two-layered snow system on the optical thickness of the upper layer calculated using the approximation given by Eq. 1 and the numerical solution of the radiative transfer equation for the nadir observation and the solar zenith angle 60°. It is assumed that the lower layer is semi-infinite. The phase function is the same in both layers (Henyey–Greenstein phase function with the average cosine of the scattering angle equal to 0.75). The single scattering albedo values in the upper (first number) and lower (second number) snow layers are given in the figure for each curve. (B) The same as in (A), except for the single scattering albedo of the lower layer, which is equal to 0.99. (C) The dependence of the error of the approximation given by Eq. 1 on the cosine of the solar zenith angle at the nadir observation for the several values of the upper layer optical thickness. All other parameters are the same as in (A). The spread along the vertical is due to seven combinations of the single scattering albedo values for the upper and lower layers given in (A).

2.2 Solar light reflectance for two-layered snow surfaces

2.2.1 Local optical properties of snow

To model the snow spectral reflectance using the formulas given above, one needs to calculate the values of the probability of photon absorption β, the average cosine of scattering angle g, and the extinction coefficient $K_{ext}$ as functions of the snow grain diameter d and the wavelength $\lambda$. This can be done either using the Maxwell electromagnetic theory or the geometrical optics approximation (Liou, 2002; Kokhanovsky, 2021). In both cases, the numerical calculations require a long time, which is not practical as far as remote sensing problems are concerned. Therefore, we use the parameterization of local optical properties derived using the geometrical optics approach for fractal ice grains (Kokhanovsky, 2021; Kokhanovsky et al., 2024). For the dry and clean snow considered in this paper, the spectral dependencies of the probability of photon absorption (PPA) β and the average cosine of scattering angle for the single light scattering event g depend on the spectral refractive index of ice, the size of ice grains, and their shape. We use the spectral refractive index of ice given by Warren and Brand (2008). The following geometrical optics equations for the randomly oriented identical fractal ice grains are used for the modeling of functions $\beta \left(\lambda \right)$ and $g\left(\lambda \right)$ (Kokhanovsky et al., 2024):

$\beta =\frac{1}{2}\left(1-\rho \right)\left(1-e^{-\sigma c}\right),(13)$

$g=g_{\infty }-\left(g_{\infty }-g_0\right)e^{-\epsilon c},(14)$

where $c=\alpha d$ is the attenuation parameter, $\alpha =\frac{4\pi \chi }{\lambda }$ is the bulk ice absorption coefficient, $\chi$ is the imaginary part of the ice refractive index at the wavelength λ, and $d$ is the effective grain diameter defined as $d=3V/S$, where V is the ice grain volume and S is the projected area of the particle on the plane perpendicular to the incident light beam (Kokhanovsky, 2021). For spherical particles, S = $\mathrm{\Sigma }/2$, where $\mathrm{\Sigma }$ is the surface area of particles. The parameters $\sigma$ and $\epsilon$ depend on the shape of ice grains and not on their size (Kokhanovsky, 2021). As $c\to \infty ,$

$\beta \left(c\to \infty \right)=\frac{1}{2}\left(1-\rho \right),g\left(c\to \infty \right)=g_{\infty }.(15)$

One can see that the parameters $\rho$ and $g_{\infty }$ describe the probability of light absorption and the average cosine of a light scattering angle for strongly absorbing particles. For such particles, one can neglect the light transmitted through the scattering particles, and the single light scattering pattern is determined exclusively by light diffraction and reflection events. In the case of randomly oriented identical nonspherical particles (say, spheroids), the light scattering pattern due to reflection coincides with that for spheres (Kokhanovsky, 2021). Therefore, the parameters $\rho$ and $g_{\infty }$ can be found assuming the spherical shape of particles. This is not the case for the parameter

$g\left(c\to 0\right)=g_0,(16)$

which strongly depends on the shape of particles and, therefore, on the type of snow. The parameters $\sigma$ and $\epsilon$ can be found using geometrical optics calculations of the parameters:

$\varphi =\frac{2\beta }{1-\rho }\text{\hspace{0.17em}and\hspace{0.5em}}\psi =\frac{g-g_0}{g_{\infty }-g_0}(17)$

as $c\to 0$. Both parameters are proportional to the value of c, and at small values of c,

$\sigma =\frac{\varphi }{c},\epsilon =\frac{\psi }{c}.(18)$

The pair (σ,ε) depends on the shape of the ice grains. In this work, we shall assume that σ = 0.8571 and ε = 0.9045 as for fractal ice grains (Kokhanovsky et al., 2024).

The parameterizations for the values of $g_{\infty }$ and ρ in terms of the real part of the ice refractive index n are (Kokhanovsky et al., 2024)

$g_{\infty }=1.008-0.11\left(\mathit{n}-1\right)\text{\hspace{0.17em}and\hspace{0.17em}}\mathrm{\rho }=0.0123+0.1622\left(\mathit{n}-1\right).(19)$

The value of $g\left(c\to 0\right)=g_0$ has been calculated using the following approximation (Kokhanovsky et al., 2024):

$g_0=0.9919-0.769\left(\mathit{n}-1\right),(20)$

which is valid for fractal ice grains.

Eqs 13, 14 are valid for arbitrary identical nonspherical randomly oriented particles, which are much larger than the wavelength of incident light ($x\to \infty ,2x\left(n-1\right)\to \infty ,x=\pi d/\lambda$ (Kokhanovsky et al., 2024)). The particles in snow can be partially oriented, and they are characterized by statistical distributions of sizes and shapes. Therefore, Eqs 13, 14 require averaging with respect to the poorly known size/shape/orientation distributions. We do not perform such averaging in this work. Instead, we use the optical diameter d, which represents the snow microstructure via Eqs 13, 14.

It should be pointed out that the extinction coefficient $k_{ext}$ of turbid media with larger particles does not depend on the wavelength and can be found using the well-known equation (Kokhanovsky, 2021)

$k_{ext}=3c_v{d}^{-1},(21)$

where $c_v$ is the volumetric concentration of scatterers. The value of $c_v$ for snow can be derived if the snow density ${\rho }_s$ is known:

$c_v=\frac{{\rho }_s}{{\rho }_i},(22)$

where ${\rho }_i$ = 0.917 $g/c{m}^3$ is the bulk ice density at the standard atmospheric pressure and a temperature of zero degrees. The average value of $c_v$ on the East Antarctica plateau is 0.355 $g/c{m}^3$ (Weinhart et al., 2020). Therefore,

$k_{ext}=A{d}^{-1},(23)$

where $A=1.16$. Assuming that the upper snow layer is homogeneous, we can derive

$\tau =A\frac{L}{d},(24)$

where L is the geometrical thickness of the layer. One can see that optical thickness approximately coincides with the number of ice grain layers on the snow thickness L. In addition, the snow optical thickness can be used to estimate the upper layer snow geometrical thickness; that is, from Eq. 24: $L=\tau d/A$.

The results of calculations of the local optical properties and also the attenuation parameter c and spherical albedo of a semi-infinite snow layer are given in Figure 2. One can see that the respective parameters vary considerably depending on the size of the particles. The same is true for the spherical albedo. This points out the fact that the determination of ice grain size is possible from optical measurements if the remote sensing channels are chosen in such a way that they provide maximal sensitivity to the snow ice grain size. Usually, one uses wavelengths 1,026 and 1,235 nm (Kokhanovsky, 2021) for the respective retrievals (see Figure 2F).

Figure 2

Figure 2. Dependence of SSA (A), PPA (B), average cosine of scattering angle (C), attenuation parameter (D), similarity parameter (E), and spherical albedo (F) on the wavelength for various effective diameters of snow grains.

2.2.2 Spectral nadir reflectance

The snow spectral nadir reflectance is shown in Figure 3. The calculations have been performed using the approximate forward model presented above at various values of the parameters ($d_1,d_2,\tau$). For the case shown in Figure 3A, it has been assumed that the effective diameter of grains in the lower layer is 0.2 mm and that it is two times smaller in the upper fresh snow layer. The spectral refractive index of ice as presented in Warren and Brand (2008) has been used. Furthermore, the reflectance increases with the optical thickness of the upper layer. This is related to the fact that particles in the upper layer are smaller and more reflective in the NIR. Figure 3B corresponds to the case of two times larger grain sizes in both layers. The spectral dependence is similar to that presented in Figure 3A, except that the dependence of reflectance on the upper layer optical thickness in the SWIR (e.g., at 2,233 nm) is almost absent. The dependence of the reflectance on the lower layer ice grain size is studied in Figure 4A at $d_1=0.1mm\text{\hspace{0.17em}and\hspace{0.17em}}\tau =10$. It is observed that the lower layer ice grain size does not influence reflectance at longer wavelengths, a feature that will be used in the retrieval procedure later on. Figure 4B presents the dependence of reflectance on the upper layer grain size at $d_2=0.6mm\text{\hspace{0.17em}and\hspace{0.17em}}\tau =10$. One can see that the reflectance is sensitive to the value of $d_1$ in a broad spectral range, which is not the case for the diameter $d_2$.

Figure 3

Figure 3. (A) Spectral dependence of the reflectance calculated using the LART model (Eqs 1–10 (Yanovitskij, 1997; Kokhanovsky et al., 2024; Weinhart et al., 2020)) for several values of snow optical thickness. The results for the semi-infinite snow with the parameters of the upper and lower layers are also given. The observation is in the nadir direction, and the solar zenith angle is equal to 60°. (B) The same as in (A), except the diameters of the upper and lower snow layers are two times larger.

Figure 4

Figure 4. (A) The same as in Fig.3, except for a snow optical thickness of 10, the diameter of ice grains in the upper layer is 0.1 mm and variable diameters of ice grains in the lower layer are in the 0.2–0.6 mm range, as shown in the figure. The calculations for the case of a homogeneous layer of ice grains with a diameter of 0.6 mm are given by the lower green line. (B) The same as in (A), except that the diameter of ice grains in the upper layer is equal to 0.1, 0.2, 0.3, and 0.4 mm and the diameter of ice grains in the lower layer is equal to 0.6 mm. The calculations for the case of a homogeneous layer of ice grains with a diameter of 0.6 mm are given by the lower green line.

3 Remote sensing of a vertically inhomogeneous snowpack using backscattered solar light

3.1 Theory

The task of the inverse problem is to derive the parameters of the model specified above (namely, $d_1$, $d_2$, and $\tau$) from the hyperspectral solar light reflectance in the spectral range of 1,000–2,450 nm, where the influence of atmospheric scattering is weak. In particular, we select the EnMAP wavelengths ${\lambda }_1$ = 1,026 nm, ${\lambda }_2$ = 1,235 nm, and ${\lambda }_3=$ 2,233 nm to derive the snow parameters. To reduce the influence of the assumed angular distribution of the single light scattering by snow crystals, we modify the coefficient $a_2$ in Eq. 12. In particular, from Eq. 12, we have $R_{1\infty }\left(r=1\right)=a_0+a_1+a_2$, where $R_{1\infty 0}\equiv R_{1\infty }\left(r=1\right)$ is the reflectance for a nonabsorbing semi-infinite light scattering layer with the assumed Henyey–Greenstein phase function at g = 0.75. To reduce the dependence of equations on the assumed single light scattering law, we propose the following representation for the coefficient $a_2$:

$a_2=R_{1\infty 0}-a_1-a_0.(25)$

Now, instead of the assumed three coefficients, we have just two ($a_0{\text{and}a}_1$), and the third coefficient ($a_2$) is calculated using Eq. 26, where $R_{1\infty 0}$ can be derived from the EnMAP data, as proposed by Kokhanovsky et al. (2019):

$R_{1\infty 0}=R_{meas}^{\gamma }\left({\lambda }_0\right)R_{meas}^{1-\gamma }\left({\lambda }_1\right),(26)$

where

$\gamma =1/\left(1-b\right),b=\sqrt{\frac{\alpha \left({\lambda }_0\right)}{\alpha \left({\lambda }_1\right)},}(27)$

where ${\lambda }_0=864nm,$ ${\lambda }_1=1026nm,\alpha \left(\lambda \right)$ is the bulk ice absorption coefficient at the wavelength $\lambda$, and $R_{meas}\left(\lambda \right)$ is the EnMAP reflectance at the wavelength $\lambda$. An important point is that Eq. 26 does not contain any unknown parameters and can be readily used for the estimation of the parameter $R_{1\infty 0}$ needed for the calculation of the coefficient $a_2,$ which is used in the retrieval procedure described below.

Furthermore, we note that the solar light at the wavelength ${\lambda }_3=$ 2,233 nm is almost not influenced by atmospheric scattering and absorption processes, and it hardly penetrates the second snow layer due to the strong absorption of bulk ice at this wavelength. Therefore, we can assume that $R_{meas}\left({\lambda }_3\right)\approx R_{1\infty }$ (see Figure 4A). This makes it possible to derive the value of $r_{1\infty }$ analytically using the EnMAP measurements at a single wavelength. For this, one needs to solve the quadratic equation given below, which is derived from Eq. 12:

$R_{meas}\left({\lambda }_3\right)-a_0-a_1{r}_{1\infty }\left({\lambda }_3\right)-a_2{r}_{1\infty }^2\left({\lambda }_3\right)=0.(28)$

The answer is

$r_{1\infty }=\frac{\sqrt{a_1^2-4a_2\left(a_0-R_{meas}\left({\lambda }_3\right)\right)}-a_1}{2a_2}.(29)$

On the other hand, the following is derived using Eq. 2 for the wavelength ${\lambda }_3$:

$s\left({\lambda }_3\right)=\frac{\psi -\sqrt{{\psi }^2-4a\left(1-r_s\left({\lambda }_3\right)\right)}}{2a},(30)$

where $\psi =1+a+b{r}_{1\infty }\left({\lambda }_3\right)$. Finally, using Eq. 3 and Eqs 13, 14, we can derive the value of the diameter of ice grains in the upper layer from the value of $s\left({\lambda }_3\right)$ under the assumption that the parameter ε = σ:

$d_1=\frac{\ln W}{\sigma \alpha \left({\lambda }_3\right)},(31)$

where

$W=\frac{2A}{\sqrt{B^2-4AC}-B},(32)$

$A=pv,B=\left(1-p\right)v-\left(1-w\right)p+p/s^2\left({\lambda }_3\right),C=w+\left(1-w\right)p+p/s^2\left({\lambda }_3\right),(33)$

$v=g_{\infty }-g_0,p=\left(1-\rho \right)/2,w=1-g_{\infty }.(34)$

In the next step, we find the values of $d_2$ and $\tau$, which minimize the function

$F={\displaystyle \sum _{m=1}^3}{\left|R_{sim}\left({\lambda }_m\right)-R_{meas}\left({\lambda }_m\right)\right|}^2,(35)$

where $R_{sim}\left({\lambda }_m\right)$ is the simulated reflectance (see Eq. 1), $R_{meas}\left({\lambda }_m\right)$ is the measured reflectance at the wavelength ${\lambda }_m$ specified above, and $d_1$ is calculated using Eq. 31. Because the derived equations are analytical, we use the simple technique of minimization based on multiple runs of the model in the range of the diameters $d_2$ between $d_{-}=d_1\text{ and }{d}_{+}=3mm$, where we take into account that the lower snow layers are characterized by the size of grains larger than those in the upper layer due to the snow metamorphism processes. The optical thickness of the upper layer is sought in the range of 1–40. The snow layers with lower values of $\tau$ are horizontally inhomogeneous, and Eq. 1 cannot be applied. The snow layers with optical thicknesses of 40 can be considered semi-infinite in the spectral range under study (see Figure 1). Therefore, for such cases, the values of $d_2$ and $\tau$ cannot be derived, and therefore, $d_2=d_1$ and $\tau \to \infty$ (a single snow layer case).

After finding the pair $\left(d_2,\tau \right)$, we use the second iteration to improve the previously derived value of $d_1$ not using the approximation of a semi-infinite layer. Instead, we use the following transcendent equation at the wavelength $\lambda$ = ${\lambda }_3$ to derive $d_1$:

$R_{meas}\left(\lambda \right)-R_{1\infty }\left(\lambda \right)+t_1\left(\lambda \right)\exp \left(-x_1\left(\lambda \right)-y_1\left(\lambda \right)\right)u\left(\xi \right)u\left(\eta \right)-\frac{t_1^2\left(\lambda \right)u\left(\xi \right)u\left(\eta \right)r_{2\infty }\left(\lambda \right)}{1-r_1\left(\lambda \right)r_{2\infty }\left(\lambda \right)}=0.(36)$

This equation follows directly from Eqs 1, 7. The values of $\left(d_2,\tau \right)$ derived from the minimization of the function given by Eq. 35 are used in Eq. 36. The sought parameter $d_1$ is derived numerically using the Brent method (Brent, 1973) under the assumption that the lower and upper boundaries for the search of the final value of $d_1$ are 1.5 times lower and larger than the value of $d_1$ derived in the first iteration. Afterward, we derive the improved values of $d_2$ and $\tau$ solving Eq. 36 at wavelengths $\lambda$ = ${\lambda }_1$ and $\lambda$ = ${\lambda }_2$, respectively. In the case of the $d_2$ determination, the Brent method (Brent, 1973) is used. As far as the optical thickness of the upper layer is concerned, it can be found in Eq. 36 analytically. Indeed, this equation can be re-written in the following form with respect to the value of $z=\exp \left({\kappa }_1\tau \right)$:

$A{z}^2-Bz+C=0,(37)$

where

$A=\left(1-r_{2\infty }{e}^{-y_1}\right)p_1{e}^{2y_1},(38)$

$B=2\left(1-r_{2\infty }ch{y}_1\right)p_1+2\left(r_{2\infty }-e^{y_1}\right)p_{2}sh{y}_1+4p_{3}s{h}^2{y}_1,(39)$

$C=\left(1-r_{2\infty }{e}^{y_1}\right)p_1{e}^{-2y_1}+2\left(r_{2\infty }-e^{-y_1}\right)p_{2}sh{y}_1,(40)$

$p_1=R_{meas}\left({\lambda }_1\right)-R_{1\infty },p_2=u\left(\xi \right)u\left(\eta \right)e^{-y_1},p_3=u\left(\xi \right)u\left(\eta \right)R_{2\infty }.(41)$

Eq 37 has the following physical solution:

$z=\frac{B-\sqrt{B^2-4AC}}{2A}.(42)$

Therefore, the optical thickness of the upper snow layer is

$\tau =\frac{\ln z}{2{\kappa }_1\left({\lambda }_2\right)}.(43)$

In the final step, we make the last iteration to find the updated value of $d_2$ from Eq. 36 at $\lambda$ = ${\lambda }_1$ and use the values of $d_1$ and $\tau$ derived at previous steps of the retrieval procedure.

The determination of the parameters $d_1$, $d_2$, and $\tau$, as discussed above, makes it possible to derive the bottom-of-atmosphere reflectance using Eq. 1 at any EnMAP channel. Therefore, the atmospheric correction can be performed in a fast, robust, and accurate way. The spectral spherical albedo r and plane albedo $r_p$ (Kokhanovsky, 2021) can be estimated using the integration of Eq. 1. After integration, (Kokhanovsky, 2021)

$r_p=r_{p1}\left(\tau \right)+\frac{t_1^2\left(\tau \right)u\left(\xi \right)r_{2\infty }}{1-r_1\left(\tau \right)r_{2\infty }},(44)$

$r=r_1\left(\tau \right)+\frac{t_1^2\left(\tau \right)r_{2\infty }}{1-r_1\left(\tau \right)r_{2\infty }}.(45)$

If needed, the plane and spherical broadband albedo can also be derived by integrating Eqs 44, 45 with respect to the wavelength (Kokhanovsky, 2021). In the case of layered snow, this procedure will lead to more accurate values of the spectral reflectance above 1,400 nm, and therefore, more accurate values of the broadband albedo will be derived. The same is true for the snow absorptance at diffuse $\left(a=1-r\right)$ and direct $\left(a_p=1-r_p\right)$ light illumination conditions.

3.2 Retrieval of the layered snow characteristics in the vicinity of the Concordia research station

We have applied the technique described above to the 30 × 30 km area around the Concordia Research Station (75.09978 S, 123.332196 E) in Antarctica. The EnMAP measurements were performed on 21 December 2023 (00:47:46 UTC, 08:47:46 UTC local time). The station is located on the Antarctic Plateau at a height of 3.2 km above sea level. The EnMAP L1C radiance data have been used. The radiance has been transformed into reflectance, as explained by Kokhanovsky et al. (2023).

The meteorological parameters at the ground derived from a Vaisala WXT520 weather station installed above the Physics shelter, Concordia Station (http://refir.fi.ino.it/rtDomeC/), are presented in Figure 5A. In particular, one can see that the air temperature is $-{33}^{\circ }C,$ the wind speed is quite low (5m/c), and its direction is close to 250° at the moment of measurements. The average wind direction can also be assessed using EnMAP data. We have selected the EnMAP channel 2,233 nm and identified all pixels within the plume with reflectance values larger than 0.22, which is larger than the reflectance for snow pixels not covered by the plume. The location of these pixels is shown in Figure 5B. They represent the plume in the direction of 4° from the OX axis (the OX axis corresponds to the 270° wind direction). Therefore, we derive that the average wind direction during the plume formation’s existence was 274°, which is close to the data shown in Figure 5A (at 8:47 local time).

Figure 5

Figure 5. (A) Meteorological parameters at the ground. (B) Wind direction estimated from EnMAP measurements. The coefficient $\mathrm{\Upsilon }$ = 0.07 shown in this figure gives the value of $\tan \left(\zeta \right),$ where $\zeta$ is the deviation from the OX axis (the exact west–east wind direction). We find that $\zeta =4^{\circ }$ and the average wind direction during the plume existence time is ${274}^{\circ }$, which coincides with the plume direction in the EnMAP imagery shown in Figure 6C.

The measured spectral top-of-atmosphere reflectance at the three channels used in the retrievals is shown in Figure 6. The station is clearly seen in all maps given in Figure 6. The track leading to the station is visible in the maps at 1,026 and 1,235 nm. This is not the case at 2,233 nm, which could be related to the fact that the track is covered by the fresh snow layer and that the radiation at this wavelength hardly penetrates the snow depth at which the track is located.

Figure 6

Figure 6. (A) Reflectance R (1,026 nm) map. (B) Reflectance R (1,235 nm) map (C). Reflectance R (2,233 nm) map.

The measured and fitting spectral reflectances (for the assumed single semi-infinite snow layer and a two-layered snowpack) at the location (74.9759 S, 123.0508 E, upper left corner in Figure 6A) are shown in Figure 7A. The calculations for the single semi-infinite snow layer have been done using Eq. 12. The value of the diameter of ice grains $d_{\infty }$ has been derived using Eq. 12 at a wavelength of 1,026 nm (Kokhanovsky et al., 2024). It has been found that $d_{\infty }=0.28mm.$ The parameters of the two-layered snow system for the studied case found by us solving the radiative transfer inverse problem, as discussed above, were $d_1=0.14mm$ and $d_2=0.39mm$. It is derived that $d_{\infty }\approx \left(d_1+d_2\right)/2$ for the case studied. One can see that the one-layer semi-infinite snow model with just two free parameters ($d_{\infty },R_{1\infty 0}$) is capable of describing the measured EnMAP reflectance (outside gaseous absorption bands) quite well. In particular, this model can be used to predict surface snow reflectance at wavelengths smaller than 1,200 nm very accurately. Therefore, the atmospheric correction can be performed in an easy way, avoiding searching for information on atmospheric constituents at the location studied. This is quite an important result as far as the atmospheric correction of VNIR EnMAP channels over snow is concerned. The two-layered snow and one-layer snow models produce very close results in the VNIR region of the electromagnetic spectrum. This is related to the fact that the upper snow layer is thin and that radiation at this spectral range penetrates deeper into the snow with very small influences due to the thin layer with other ice grain microstructure at the top. The situation is radically different in the SWIR (e.g., see the 1,500–2,500 nm spectral range in Figure 6A), where radiation hardly penetrates the snow and only the upper snow layer influences the solar light reflectance. As we have found using optical measurements, the diameter of grains in this upper layer is two times smaller than the value of $d_{\infty }$. This leads to an underestimation of the reflectance in the SWIR in the framework of the single-layer snow model if the two-layered snow is present on the ground during the measurements. This feature has already been observed for a long time in the ground measurements of solar light reflectance from snow surfaces in Antarctica. In particular, the ground observations and modeling of Antarctic snow surface albedo performed by Grenfell et al. (1994) show that the albedo spectra in a broad spectral range can be fitted using the two-layered snow model, which confirms our findings. The derived diameters of particles were 0.06 and 0.2 mm in the upper and lower snow layers, respectively, which are approximately two times smaller than those found in our findings. This difference could be related to the different microstructure of snow during respective ground and satellite measurements. Grenfell et al. (1994) utilized the assumption of spherical ice grains, which is different from our more realistic assumption of fractal ice particles. This could also contribute to the found differences found. The derived value of the upper snow layer geometrical thickness was 0.25 mm, which is also approximately two times smaller than that in our case. Grenfell et al. (1994) did not provide the value of $\tau$. However, taking into account that $\tau \approx L/d$ with both L and d being two times smaller in Grenfell et al. (1994) than those in our case, we arrive at the value of $\tau \approx 4.2$, which is similar to the result we derived. One can see that the SWIR reflectance is sensitive to the presence of the extremely thin snow layer at the top of the old snow.

Figure 7

Figure 7. (A) Measured snow spectral reflectance in the nadir direction. The results of numerical simulations according to the developed forward model for single- and two-layered snow are also shown. The solar zenith angle is 56°. (B) Measured snow spectral reflectance in the nadir direction. The results of numerical simulations according to the developed forward model for the TOA and BOA spectral reflectances in the case of a two-layered snowpack are also shown. The solar zenith angle is 56°. (C). The temporal variation of the PWV at the Concordia station. The vertical line corresponds to the time of measurement. The detected cloud optical thickness (cloud OT) is below 0.1. The lower panel shows the accuracy of fit used for the determination of the PWV. (D). The variation in the total ozone derived from ground measurements (UV-RAD) and ERA5 (https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5) at the Concordia Research Station. ERA5 is the fifth generation of the European Center for Medium-range Weather Forecasts (ECMWF) atmospheric reanalysis of the global climate, covering the period from January 1940 to the present. The rapid fluctuations of the TOC at the time of EnMAP measurements (December 21, 00:47:46 UTC) are clearly seen.

We also simulated the TOA EnMAP spectrum using the modified version of the radiative transfer model SNOWTRAN (Kokhanovsky et al., 2024). In particular, the LART model described in this work has been incorporated into the SNOWTRAN model. The results are given in Figure 7B, where we used the values of the precipitable water vapor PWV = 0.6 mm measured on the ground using the technique described by Kokhanovsky et al. (2024) (see Fig. 7c). The total ozone column was assumed to be equal to 360DU, which is close to the ground-measured value of total ozone (338DU). The radiometer used for the ground measurements of the total ozone is described by Kokhanovsky et al. (2024). We have found that the assumed value of the total ozone column (7% higher value as measured on the ground) better describes EnMAP observations, which is possibly due to some differences in the location of the EnMAP ground scene with a spectrum shown in Figure 7A and the location of the instrument. The fluctuations in the temporal variation of the total ozone (see Figure 7D) and errors in ground measurements could also contribute to this difference.

Figure 7B shows that the measured and simulated TOA spectra are close to each other. This leads to the conclusion that the radiative transfer in the atmosphere-underlying surface system over the area around the Concordia Research Station is well-understood. The accurate performance of the EnMAP mission over bright surfaces is confirmed. Both TOC and PWV can be retrieved from ENMAP measurements, as demonstrated by Kokhanovsky et al. (2024) and confirmed by Figure 7B. The aerosol load in the area close to the station is very low due to the absence of atmospheric pollution sources and the high elevation of the station. Light scattering by molecules is far more important at this location than that by atmospheric aerosols. We have assumed that the aerosol optical thickness at 550 nm is 0.01 and the Angström exponent is equal to 1.6, which is typical for the aerosol over the Concordia Research Station (Six et al., 2005).

The maps of the retrieved snow characteristics are given in Figure 8. The respective statistical characteristics of the retrieved parameters are presented in Table 1. The average value of $d_{\infty }$ was 0.33 mm for the area shown in Figure 8. One can see that the diameter of ice grains in the lower snow layer is approximately two times larger than the diameter of ice crystals in the thin upper layer, which has a thickness of approximately 0.5 mm. The snow grain sizes in the upper snow layer are somewhat larger in the lower right corner of the respective map, with smaller crystals in the plume of the Concordia station, which could be expected on physical grounds due to the smaller size of crystals in the plume. Generally speaking, the spatial variation of the snow parameters is quite low, as has already been reported using ground observations (Gay et al., 2002). In particular, Gay et al. (2002) analyzed more than 500 samples taken in the interior of the Antarctic ice sheet and found that the diameters of particles are in the range of 0.2–0.4 mm, which is similar to our findings based on the analysis of the hyperspectral satellite imagery. The statistical distributions of the derived grain diameters are given in Figure 9. One can see that the spatial variation of snow parameters in the interior of the ice sheet is quite low. We present the spatial distribution of the derived upper layer snow geometrical and optical thicknesses in Figure 10. The average upper layer optical thickness is 3.16, and the estimated geometrical thickness of the upper layer is 0.51 mm. One can see that the upper layer is very thin and consists of three to four layers of tiny ice crystals. We may conclude that the aged snow is covered by a thin veil of tiny ice crystals at the top. This is consistent with the data presented by Turner et al. (2019), where it is stated that the highest daily precipitation at Dome C is just 2.77 mm with a 90th percentile of 0.39 mm. The thin layer of crystals is laid down and then begins to age and metamorphose with Sun and wind exposure. Gallet et al. (2011) also found that the thin layer of crystals (h = 0.25 mm) on the top of aged snow is consistent with the spectral snow reflectance observations in the range of 1,600–2,500 nm (although the crystals in their simulation were much smaller (d = 0.03 mm) than the sizes of crystals in the upper layer found by us (d = 0.16 mm, see Table 1). It should be pointed out that the values of d < 0.1 mm are not typical for the Antarctic Plateau (Sato et al., 1981; Gay et al., 2002). We should stress that the presence of the upper thin layer of ice crystals can be detected if measurements in the spectral range of 1,600–2,500 nm are available. Otherwise, the signal is dominated by the semi-infinite snow layer underneath.

Figure 8

Figure 8. Ice grain size in the upper (upper picture) and lower (lower picture) snow layers.

Table 1

Table 1. Average value and the coefficient of variance (CV = standard deviation/mean) of various snow characteristics. The last three columns show the average values of the spectral reflectances at wavelengths 1,026, 1,235, and 2,233 nm.

Figure 9

Figure 9. Statistical distributions of the derived snow diameters in the upper and lower snow layers.

Figure 10

Figure 10. Geometrical (upper panel) and optical (lower panel) upper snow layer thickness.

The estimation of the snow parameters can be performed only in the absence of clouds and pollution sources on the ground. Therefore, the derived results at the location of the plume from the station do not represent the properties of the underlying surface. In particular, Figure 8 (upper panel) shows that the sizes of particles in the exhaust plume are smaller than those of the particles on the ground as one might expect.

In addition to the snow grain size, we have derived the snow-specific surface area $\mathrm{\Sigma }$ related to the diameter of ice grains using the following simple relationship (Kokhanovsky, 2021):

$\mathrm{\Sigma }=\frac{6}{{\rho }_{i}d}.(46)$

The spatial distribution of the specific surface area in the upper layer is shown in Figure 11. The average value is 40 $m^2/kg$, and the coefficient of variance is 8%. Gallet et al. (2011) found using ground measurements that $\mathrm{\Sigma }$ is in the range of 24–55 $m^2/kg$ around Dome C, which is consistent with the retrieval of this parameter using EnMAP spaceborne observations (see Figure 11).

Figure 11

Figure 11. Snow-specific surface area in the upper snow layer.

It has been reported (Inoue et al., 2024) that the snow-specific surface area along the traverse road from the coast to Dome Fuji, Antarctica, was $45\pm 11m^2/kg$ on the traverse distance of 500–1,066 km from the coast. We have found that the snow-specific surface area is $40\pm 3m^2/kg$ in the 30 × 30 km ground scene close to the Concordia Research Station, which is 1,100 km from the coast and has a similar elevation as Dome Fuji. Therefore, we can conclude that the variability of the specific surface area of the surface snow on the vast Antarctic Plateau is rather low, which is related to the unique environmental conditions in the area characterized by dry air, rare storms, and low temperatures.

4 Conclusion

We have proposed the fast and simple, fully analytical LART model for the calculation of the BOA and TOA spectral reflectance over layered snow surfaces in the VNIR and SWIR regions of the electromagnetic spectrum. It has been assumed that the snow is an ideally horizontal surface without sastrugi and any inclination. The close-packed media effects have been neglected. The ice particles have been modeled as Koch fractals of the second generation (Kokhanovsky, 2021). We have accounted for the light absorption by atmospheric gases such as molecular oxygen, water vapor, and ozone. The vertical inhomogeneity of the snow surface has been accounted for by assuming a two-layered snow model. We have supposed that the upper layer is optically thick and the lower layer is semi-infinite. The derived technique has been applied to EnMAP data, which show the capability of detecting thin fresh snow layers overlying aged snow using the capability of the EnMAP to measure the solar light reflectance from snow in the SWIR. Therefore, in this work, we have confirmed the possible existence of a thin upper snow layer with tiny ice crystals on the Antarctic plateau postulated by Grenfell et al. (1994) and further studied by Turner et al. (2019). This layer, being geometrically thin, has important impacts on the radiative budget in the SWIR region of the electromagnetic spectrum, leading to a larger reflectance of electromagnetic radiation from snow in the SWIR than in the case of a homogeneous layer with the grain size of the aged snow underneath.

The determination of snow characteristics using only VNIR measurements may lead to biases in the retrieved snow grain sizes and, therefore, snow albedo, which is an essential climate variable. The derived technique can be used in the improvement of spaceborne spectral and broadband albedo determination using satellite measurements, as well as for the detection and characterization of fresh snowfall events using EnMAP and other hyperspectral datasets. The EnMAP data can also be used to retrieve the values of TOC, PWV, and the direction of the wind on the ground in cases where strong aerosol emission sources are present on the ground.

We believe that the model proposed in this work will be useful for the determination of snow parameters using other hyperspectral and multispectral instrumentation, including ground-based optical spectrometers and radiometers. The technique can be easily generalized to study other optically thick turbid media with large scatterers, including soil and desert surfaces.

Future work will focus on the validation of the derived parameters, the estimation of uncertainties in retrievals, and the runs of the model for various regions and seasons.

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: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing–original draft, and writing–review and editing. MB: data curation, writing–review and editing. KS: project administration, resources, software, supervision, writing–review and editing. DE: formal analysis, investigation, software, writing–original draft, and writing–review and editing. BP: conceptualization, data curation, formal analysis, investigation, writing–original draft, and writing–review and editing. GB: conceptualization, data curation, formal analysis, investigation, methodology, visualization, writing–original draft, and writing–review and editing. RS: methodology, writing–original draft. SC: project administration, resources, supervision, methodology, investigation, writing–original draft.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article. The research was funded by the EnMAP mission science program funded within the German Earth observation program by the DLR Space Agency with resources from the German Federal Ministry for Economic Affairs and Climate Action (grant numbers 50EE1923 and 50EE2401).

Acknowledgments

The authors thank DLR for providing EnMAP data and ECMWF for providing ERA5 re-analysis data.

Conflict of interest

The 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.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Afanas’ev, V. P., Basov, A. Y., Budak, V. P., Efremenko, D. S., and Kokhanovsky, A. A. (2020). Analysis of the discrete theory of radiative transfer in the coupled “ocean–atmosphere” system: current status, problems and development prospects. J. Mar. Sci. Eng. 8 (3), 202. doi:10.3390/jmse8030202

CrossRef Full Text | Google Scholar

Brent, R. P. (1973). Algorithms for Minimization without Derivatives. Englewood Cliffs, NJ: Prentice-Hall.

Google Scholar

Chuprov, I., Konstantinov, D., Efremenko, D., Zemlyakov, V., and Gao, J. (2022). Solution of the radiative transfer equation for vertically inhomogeneous media by numerical integration solvers: comparative analysis. Light and Eng., 21–30. doi:10.33383/2022-030

CrossRef Full Text | Google Scholar

Efremenko, D., Doicu, A., Loyola, D., and Trautmann, T. (2013). Acceleration techniques for the discrete ordinate method. J. Quantitative Spectrosc. Radiat. Transf. 114, 73–81. doi:10.1016/j.jqsrt.2012.08.014

CrossRef Full Text | Google Scholar

Efremenko, D., Molina García, V., Gimeno García, S., and Doicu, A. (2017). A review of the matrix-exponential formalism in radiative transfer. J. Quantitative Spectrosc. Radiat. Transf. 196, 17–45. doi:10.1016/j.jqsrt.2017.02.015

CrossRef Full Text | Google Scholar

Gallet, J.-C., Domine, F., Arnaud, L., Picard, G., and Savarino, J. (2011). Vertical profile of the specific surface area and density of the snow at Dome C and on a transect to Dumont D'Urville, Antarctica – albedo calculations and comparison to remote sensing products. Cryosphere 5, 631–649. doi:10.5194/tc-5-631-2011

CrossRef Full Text | Google Scholar

Gay, M., Fily, M., Genthon, C., Frezzotti, M., Oerter, H., and Winther, J.-G. (2002). Snow grain size measurements in Antarctica. J. Glaciol. 48, 527–535. doi:10.3189/172756502781831016

CrossRef Full Text | Google Scholar

Grenfell, T. C., Warren, S. G., and Mullen, P. C. (1994). Reflection of solar radiation by the Antarctic snow surface at ultraviolet, visible, and near-infrared wavelengths. J. Geophys. Res. 99 (D9), 18669–18684. doi:10.1029/94jd01484

CrossRef Full Text | Google Scholar

Guanter, L., Kaufmann, H., Segl, K., Foerster, S., Rogass, C., Chabrillat, S., et al. (2015). The EnMAP spaceborne imaging spectroscopy mission for Earth Observation. Remote Sens. 7, 8830–8857. doi:10.3390/rs70708830

CrossRef Full Text | Google Scholar

Henyey, L. G., and Greenstein, J. L. (1941). Diffuse radiation in the Galaxy. Astrophysical J. 93, 70–83. doi:10.1086/144246

CrossRef Full Text | Google Scholar

Hu, Y.-X., Wielicki, B., Lin, B., Gibson, B., Tsay, S.-C., Stamnes, K., et al. (2000). δ-Fit: a fast and accurate treatment of particle scattering phase functions with weighted singular-value decomposition least-squares fitting. J. Quantitative Spectrosc. Radiat. Transf. 65 (4), 681–690. doi:10.1016/s0022-4073(99)00147-8

CrossRef Full Text | Google Scholar

Inoue, R., Aoki, T., Fujita, S., Tsutaki, S., Motoyama, H., Nakazawa, F., et al. (2024). Spatiotemporal variation in the specific surface area of surface snow measured along the traverse route from the coast to Dome Fuji, Antarctica. EGUsphere. doi:10.5194/egusphere-2024-769

CrossRef Full Text | Google Scholar

Kokhanovsky, A., Brell, M., Segl, K., and Chabrillat, S. (2024). SNOWTRAN: a fast radiative transfer model for polar hyperspectral remote sensing applications. Remote Sens. 16, 334. doi:10.3390/rs16020334

CrossRef Full Text | Google Scholar

Kokhanovsky, A., Di Mauro, B., Garzonio, R., and Colombo, R. (2021). Retrieval of dust properties from spectral snow reflectance measurements. Front. Environ. Sci. 9. doi:10.3389/fenvs.2021.644551

CrossRef Full Text | Google Scholar

Kokhanovsky, A., Lamare, M., Danne, O., Brockmann, C., Dumont, M., Picard, G., et al. (2019). Retrieval of snow properties from the sentinel-3 ocean and land colour instrument. Remote Sens. 11, 2280. doi:10.3390/rs11192280

CrossRef Full Text | Google Scholar

Kokhanovsky, A., and Rozanov, V. V. (2012). Droplet vertical sizing in warm clouds using passive optical measurements from a satellite. Atmos. Meas. Tech. 5, 517–528. doi:10.5194/amt-5-517-2012

CrossRef Full Text | Google Scholar

Kokhanovsky, A. A. (2006). Cloud Optics. Berlin, Germany: Springer.

Google Scholar

Kokhanovsky, A. A. (2021). Snow optics. Cham, Switzerland: Springer Nature.

Google Scholar

Kokhanovsky, A. A., Brell, M., Segl, K., Bianchini, G., Lanconelli, C., Lupi, A., et al. (2023). First retrievals of surface and atmospheric properties using EnMAP measurements over Antarctica. Remote Sens. 15, 3042. doi:10.3390/rs15123042

CrossRef Full Text | Google Scholar

Li, W., Stamnes, K., Chen, B., and Xiong, X. (2011). Snow grain size retrieved from near-infrared radiances at multiple wavelengths. Geophys. Res. Let. 28 (9), 1699–1702. doi:10.1029/2000gl011641

CrossRef Full Text | Google Scholar

Liou, K.-N. (2002). An Introduction to Atmopsheric Radiation. New York, NY, USA: Academic Press.

Google Scholar

Mishchenko, M. I., Dlugach, J. M., Yanovitskij, E. G., and Zakharova, N. T. (1999). Bidirectional reflectance of flat, optically thick particulate layers: An efficient radiative transfer solution and applications to snow and soil surfaces. J. Quantitative Spectrosc. Radiat. Transf. 63, 409–432. doi:10.1016/s0022-4073(99)00028-x

CrossRef Full Text | Google Scholar

Nakajima, T., and Tanaka, M. (1988). Algorithms for radiative intensity calculations in moderately thick atmospheres using a truncation approximation. J. Quantitative Spectrosc. Radiat. Transf. 40 (1), 51–69. doi:10.1016/0022-4073(88)90031-3

CrossRef Full Text | Google Scholar

Sato, N., Kikuchi, K., Barnard, S. C., and Hogan, A. W. (1981). Some characteristic properties of ice crystal precipitation in the summer season at South Pole station, Antarctica. J. Meteorol. Soc. Jpn. 59, 7772–7780. doi:10.2151/jmsj1965.59.5_772

CrossRef Full Text | Google Scholar

Six, D., Fily, M., Blarel, L., and Goloub, P. (2005). First aerosol optical thickness measurements at Dome C (East Antarctica), summer season 2003–2004. Atmos. Environ. 39, 5041–5050. doi:10.1016/j.atmosenv.2005.05.010

CrossRef Full Text | Google Scholar

Sobolev, V. V. (1975). Light Scattering in Planetary Atmospheres. Moscow: Nauka.

Google Scholar

Storch, T., Honold, H.-P., Chabrillat, S., Habermeyer, M., Tucker, P., Brell, M., et al. (2023). The EnMAP imaging spectroscopy mission towards operations. Rem. Sens. Env. 294, 113632. doi:10.1016/j.rse.2023.113632

CrossRef Full Text | Google Scholar

Turner, J., Phillips, T., Thamban, M., Rahaman, W., Marshall, G. J., Wille, J. D., et al. (2019). The dominant role of extreme precipitation events in Antarctic snow fall variability. Geophys. Res. Let. 46, 3502–3511. doi:10.1029/2018gl081517

CrossRef Full Text | Google Scholar

van de Hulst, H. C. (1974). The spherical albedo of a planet covered with a homogeneous cloud layer. Astron. Astrophys. 35, 209–214.

Google Scholar

Warren, S., and Brand, R. E. (2008). Optical constants of ice from the ultraviolet to the microwave: a revised compilation. J. Geophys. Res. 113, D14. doi:10.1029/2007jd009744

CrossRef Full Text | Google Scholar

Weinhart, A. H., Freitag, J., Hörhold, M., Kipfstuhl, S., and Eisen, O. (2020). Representative surface snow density on the East Antarctic plateau. Cryosphere 14, 3663–3685. doi:10.5194/tc-14-3663-2020

CrossRef Full Text | Google Scholar

Wiscombe, W. J. (1977). The Delta–M Method: rapid yet accurate radiative flux calculations for strongly asymmetric phase functions. J. Atmos. Sci. 34 (9), 1408–1422. doi:10.1175/1520-0469(1977)034<1408:tdmrya>2.0.co;2

CrossRef Full Text | Google Scholar

Yanovitskij, E. G. (1997). Light Scattering in Inhomogeneous Atmospheres. Berlin: Springer.

Google Scholar

Zege, E. P., Ivanov, A. P., and Katsev, I. L. (1991). Image Transfer Through Light Scattering Media. Berlin, Germany: Springer.

Google Scholar