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

Abstract

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

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 , the cosine of the viewing zenith angle , and the relative azimuthal angle are modeled as (Sobolev, 1975; Liou, 2002; Kokhanovsky, 2006)where is the reflectivity of the upper layer with the snow optical thickness (SOT) for the case of a black underlying surface, is the snow transmittance for the upper layer for the light propagation from the snow surface down, is the same quantity for the light propagation from the second snow layer up to the snow surface, is the spherical albedo of the upper layer, and 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, , where 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 is given by the very accurate van de Hulst approximation (van de Hulst, 1974):where a = 0.139, b = 1.17, andis the similarity parameter, is the single scattering albedo (SSA), and 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 and in the following form (Sobolev, 1975; Zege et al., 1991; Kokhanovsky, 2006; Kokhanovsky, 2021):where is the escape function and 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)where 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 and coefficients and 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 and its accurate calculation are not required. In the case of strong absorption, . 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 , and (Zege et al., 1991):wherewhere is the single scattering albedo and is the average cosine of the single light scattering angle in the upper snow layer. The parameter 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 (-)) at large values of . The value of is proportional to the similarity parameter for the upper layer as (). From Eq. 8, and as as it should be (Zege et al., 1991). The escape function is derived as follows (Kokhanovsky et al., 2021):

It follows that as as it should be. In addition, the same asymptotic solution can be derived as , andwhere 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) , which is the case for the weak absorption approximation considered here. In particular it follows, from Eq. 11 as that , where . It follows from Eq. 2, in the same case, that , where . It should be pointed out that

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 with the parameters of the upper layer. The function 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 in the nadir direction and the value of spherical albedo :

This simple equation has been derived by Kokhanovsky et al. (2024) under the assumption that , which is characteristic for snow media in the optical range of the electromagnetic spectrum, where EnMAP performs the measurements. The parameters 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 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

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 as functions of the snow grain diameter d and the wavelength . 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 and (Kokhanovsky et al., 2024):where is the attenuation parameter, is the bulk ice absorption coefficient, is the imaginary part of the ice refractive index at the wavelength λ, and is the effective grain diameter defined as , 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 =, where is the surface area of particles. The parameters and depend on the shape of ice grains and not on their size (Kokhanovsky, 2021). As

One can see that the parameters and 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 and can be found assuming the spherical shape of particles. This is not the case for the parameterwhich strongly depends on the shape of particles and, therefore, on the type of snow. The parameters and can be found using geometrical optics calculations of the parameters:as . Both parameters are proportional to the value of c, and at small values of c,

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 and ρ in terms of the real part of the ice refractive index n are (Kokhanovsky et al., 2024)

The value of has been calculated using the following approximation (Kokhanovsky et al., 2024):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 ( (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 of turbid media with larger particles does not depend on the wavelength and can be found using the well-known equation (Kokhanovsky, 2021)where is the volumetric concentration of scatterers. The value of for snow can be derived if the snow density is known:where = 0.917 is the bulk ice density at the standard atmospheric pressure and a temperature of zero degrees. The average value of on the East Antarctica plateau is 0.355 (Weinhart et al., 2020). Therefore,where . Assuming that the upper snow layer is homogeneous, we can derivewhere 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: .

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

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 (). 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 . 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 . One can see that the reflectance is sensitive to the value of in a broad spectral range, which is not the case for the diameter .

FIGURE 3

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

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, , , and ) 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 = 1,026 nm, = 1,235 nm, and 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 in Eq. 12. In particular, from Eq. 12, we have , where 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 :

Now, instead of the assumed three coefficients, we have just two (), and the third coefficient () is calculated using Eq. 26, where can be derived from the EnMAP data, as proposed by Kokhanovsky et al. (2019):wherewhere is the bulk ice absorption coefficient at the wavelength , and is the EnMAP reflectance at the wavelength . An important point is that Eq. 26 does not contain any unknown parameters and can be readily used for the estimation of the parameter needed for the calculation of the coefficient which is used in the retrieval procedure described below.

Furthermore, we note that the solar light at the wavelength 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 (see Figure 4A). This makes it possible to derive the value of 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:

The answer is

On the other hand, the following is derived using Eq. 2 for the wavelength :where . 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 under the assumption that the parameter ε = σ:where

In the next step, we find the values of and , which minimize the functionwhere is the simulated reflectance (see Eq. 1), is the measured reflectance at the wavelength specified above, and 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 between , 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 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 and cannot be derived, and therefore, and (a single snow layer case).

After finding the pair , we use the second iteration to improve the previously derived value of not using the approximation of a semi-infinite layer. Instead, we use the following transcendent equation at the wavelength = to derive :

This equation follows directly from Eqs 1, 7. The values of derived from the minimization of the function given by Eq. 35 are used in Eq. 36. The sought parameter 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 are 1.5 times lower and larger than the value of derived in the first iteration. Afterward, we derive the improved values of and solving Eq. 36 at wavelengths = and = , respectively. In the case of the 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 :where

Eq 37 has the following physical solution:

Therefore, the optical thickness of the upper snow layer is

In the final step, we make the last iteration to find the updated value of from Eq. 36 at = and use the values of and derived at previous steps of the retrieval procedure.

The determination of the parameters , , and , 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 (Kokhanovsky, 2021) can be estimated using the integration of Eq. 1. After integration, (Kokhanovsky, 2021)

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 and direct 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 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

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

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 has been derived using Eq. 12 at a wavelength of 1,026 nm (Kokhanovsky et al., 2024). It has been found that 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 and . It is derived that for the case studied. One can see that the one-layer semi-infinite snow model with just two free parameters () 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 . 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 . However, taking into account that 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 , 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

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 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 9

FIGURE 10

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 related to the diameter of ice grains using the following simple relationship (Kokhanovsky, 2021):

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

FIGURE 11

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 on the traverse distance of 500–1,066 km from the coast. We have found that the snow-specific surface area is 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.

References

• 1

  Afanas’evV. P.BasovA. Y.BudakV. P.EfremenkoD. S.KokhanovskyA. 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. 10.3390/jmse8030202

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

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

  • Google Scholar
  • View reference in article

• 3

  ChuprovI.KonstantinovD.EfremenkoD.ZemlyakovV.GaoJ. (2022). Solution of the radiative transfer equation for vertically inhomogeneous media by numerical integration solvers: comparative analysis. Light and Eng., 21–30. 10.33383/2022-030

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  EfremenkoD.DoicuA.LoyolaD.TrautmannT. (2013). Acceleration techniques for the discrete ordinate method. J. Quantitative Spectrosc. Radiat. Transf.114, 73–81. 10.1016/j.jqsrt.2012.08.014

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  EfremenkoD.Molina GarcíaV.Gimeno GarcíaS.DoicuA. (2017). A review of the matrix-exponential formalism in radiative transfer. J. Quantitative Spectrosc. Radiat. Transf.196, 17–45. 10.1016/j.jqsrt.2017.02.015

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  GalletJ.-C.DomineF.ArnaudL.PicardG.SavarinoJ. (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. Cryosphere5, 631–649. 10.5194/tc-5-631-2011

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  GayM.FilyM.GenthonC.FrezzottiM.OerterH.WintherJ.-G. (2002). Snow grain size measurements in Antarctica. J. Glaciol.48, 527–535. 10.3189/172756502781831016

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  GrenfellT. C.WarrenS. G.MullenP. 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. 10.1029/94jd01484

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  GuanterL.KaufmannH.SeglK.FoersterS.RogassC.ChabrillatS.et al (2015). The EnMAP spaceborne imaging spectroscopy mission for Earth Observation. Remote Sens.7, 8830–8857. 10.3390/rs70708830

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  HenyeyL. G.GreensteinJ. L. (1941). Diffuse radiation in the Galaxy. Astrophysical J.93, 70–83. 10.1086/144246

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  HuY.-X.WielickiB.LinB.GibsonB.TsayS.-C.StamnesK.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. 10.1016/s0022-4073(99)00147-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  InoueR.AokiT.FujitaS.TsutakiS.MotoyamaH.NakazawaF.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. 10.5194/egusphere-2024-769

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  KokhanovskyA.BrellM.SeglK.ChabrillatS. (2024). SNOWTRAN: a fast radiative transfer model for polar hyperspectral remote sensing applications. Remote Sens.16, 334. 10.3390/rs16020334

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  KokhanovskyA.Di MauroB.GarzonioR.ColomboR. (2021). Retrieval of dust properties from spectral snow reflectance measurements. Front. Environ. Sci.9. 10.3389/fenvs.2021.644551

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  KokhanovskyA.LamareM.DanneO.BrockmannC.DumontM.PicardG.et al (2019). Retrieval of snow properties from the sentinel-3 ocean and land colour instrument. Remote Sens.11, 2280. 10.3390/rs11192280

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  KokhanovskyA.RozanovV. V. (2012). Droplet vertical sizing in warm clouds using passive optical measurements from a satellite. Atmos. Meas. Tech.5, 517–528. 10.5194/amt-5-517-2012

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

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

  • Google Scholar
  • View reference in article

• 18

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

  • Google Scholar
  • View reference in article

• 19

  KokhanovskyA. A.BrellM.SeglK.BianchiniG.LanconelliC.LupiA.et al (2023). First retrievals of surface and atmospheric properties using EnMAP measurements over Antarctica. Remote Sens.15, 3042. 10.3390/rs15123042

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  LiW.StamnesK.ChenB.XiongX. (2011). Snow grain size retrieved from near-infrared radiances at multiple wavelengths. Geophys. Res. Let.28 (9), 1699–1702. 10.1029/2000gl011641

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

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

  • Google Scholar
  • View reference in article

• 22

  MishchenkoM. I.DlugachJ. M.YanovitskijE. G.ZakharovaN. 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. 10.1016/s0022-4073(99)00028-x

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

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

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  SatoN.KikuchiK.BarnardS. C.HoganA. 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. 10.2151/jmsj1965.59.5_772

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  SixD.FilyM.BlarelL.GoloubP. (2005). First aerosol optical thickness measurements at Dome C (East Antarctica), summer season 2003–2004. Atmos. Environ.39, 5041–5050. 10.1016/j.atmosenv.2005.05.010

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

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

  • Google Scholar
  • View reference in article

• 27

  StorchT.HonoldH.-P.ChabrillatS.HabermeyerM.TuckerP.BrellM.et al (2023). The EnMAP imaging spectroscopy mission towards operations. Rem. Sens. Env.294, 113632. 10.1016/j.rse.2023.113632

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  TurnerJ.PhillipsT.ThambanM.RahamanW.MarshallG. J.WilleJ. D.et al (2019). The dominant role of extreme precipitation events in Antarctic snow fall variability. Geophys. Res. Let.46, 3502–3511. 10.1029/2018gl081517

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

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

  • Google Scholar
  • View reference in article

• 30

  WarrenS.BrandR. E. (2008). Optical constants of ice from the ultraviolet to the microwave: a revised compilation. J. Geophys. Res.113, D14. 10.1029/2007jd009744

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  WeinhartA. H.FreitagJ.HörholdM.KipfstuhlS.EisenO. (2020). Representative surface snow density on the East Antarctic plateau. Cryosphere14, 3663–3685. 10.5194/tc-14-3663-2020

  • CrossRef
  • Google Scholar

• 32

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

  • CrossRef
  • Google Scholar

• 33

  YanovitskijE. G. (1997). Light Scattering in Inhomogeneous Atmospheres. Berlin: Springer.

  • Google Scholar

• 34

  ZegeE. P.IvanovA. P.KatsevI. L. (1991). Image Transfer Through Light Scattering Media. Berlin, Germany: Springer.

  • Google Scholar