First Retrievals of Surface and Atmospheric Properties Using EnMAP Measurements over Antarctica

Abstract

The paper presents the first retrievals of clean snow properties using spaceborne hyperspectral observations via the Environmental Mapping and Analysis Program (EnMAP). The location close to the Concordia station at the Dome C Plateau (Antarctica) was selected. At this location, the atmospheric effects (except molecular light scattering and absorption) are weak, and the simplified atmospheric correction scheme could be applied. The ice grain size, snow specific surface area, and snow spectral and broadband albedos were retrieved using single-view EnMAP measurements. In addition, we propose a technique to retrieve trace gas concentrations (e.g., water vapor and ozone) from EnMAP observations over the snow surfaces. A close correspondence of satellite and ground-measured parameters was found.

1. Introduction

The cryosphere is an integral part of the terrestrial ecosystem with important linkages and feedbacks generated through its influence on moisture fluxes, hydrology, hydrosphere, clouds, radiation, the water cycle, and climate change due to temporal changes in snow/ice cover and albedo. Therefore, snow and ice properties are monitored using ground and satellite instrumentation [1,2,3]. The measurements are performed at various temporal and spatial scales using in situ measurements and also passive and active remote sensing instrumentation in a broad spectral range. The high spatial resolution is highly relevant for studies of the cryosphere due to the horizontal variability of snow properties and impurity load (dust outbreaks, algae, and soot). The instruments with a coarser spatial resolution (say, 0.3–1 km) are not capable of resolving fine scales of snow and ice variability. In particular, the accurate information on specific snow features (e.g., spatial distribution of algal blooms) can be hardly detected using measurements performed on the scale of 0.3–1.0 km. The Environmental Mapping and Analysis Program (EnMAP) is a German hyperspectral satellite mission [4] that provides information on the evolution of aquatic and terrestrial ecosystems including the cryosphere on the spatial scale of 30 m with 224 bands in the wavelength range between 400 and 2500 nm, which means a better spectral and spatial resolution compared to many other satellite missions (e.g., the Ocean and Land Colour Instrument (OLCI) on board Sentinel-3 (S-3)). The measurements of EnMAP can be also used to assess the sub-pixel snow variability for coarse spatial resolution satellite missions, estimate the accuracy of satellite products derived on coarser spatial grids (e.g., snow fraction), and also for the calibration and validation of spatially and spectrally coarser products as optimized inputs for numerical climate products. This paper is aimed at the adaptation of the previously proposed optical snow remote sensing technique [5,6,7,8] to EnMAP measurements over areas with clean dry snow. We used the fact that the intensity of the shortwave infrared (SWIR) radiation reflected from snow depends on the snow grain size, which is the main parameter governing snow optical properties. This is also the main basis of other available snow remote sensing algorithms that rely on optical (including hyperspectral) measurements [9,10,11,12,13,14]. In particular, we derived the snow grain size, spectral and broadband snow albedo, and snow bottom-of-atmosphere reflectance (BOA) from single-view spectral top-of-atmosphere (TOA) reflectance measurements over Antarctica. The snow specific surface area was derived using the retrieved grain size. It was assumed that snow did not contain impurities and could be modeled as a horizontally homogeneous plane-parallel semi-infinite light-scattering layer. Warren et al. [15] showed that the soot content in snow at Dome C is smaller than 3.5 ppb and generally decreases with the snow depth. Such small values of soot content cannot be detected with optical instrumentation. There is some variation in the snow grain size with depth due to the layered nature of the snow surface. Solar light with different wavelengths has different penetration depths in snow; therefore, different snow volumes are sensed using multiple wavelengths. The vertical snow grain size variability was ignored because we retrieved the snow grain size using EnMAP measurements [4,16] located in the relatively narrow spectral range (1026–1235 nm) where the vertical variability of the snow grain size at Dome C, which has relatively stable environmental conditions, could be ignored.

This technique has been successfully applied to measurements by the Ocean and Land Colour Instrument on board Sentinel-3 as reported in [8]. The validation results are given in [6,7,8,17]. The respective software is available in both Python and Fortran implementations at https://github.com/GEUS-SICE/pySICE (accessed on 20 April 2023) [18]. In addition, we derived the total ozone column (TOC) and precipitable water vapor (PWV) over the selected area in Antarctica from EnMAP measurements. The retrieval was based on the analysis of the depth of the gaseous absorption bands as seen in the reflected light, which is a common technique in spaceborne trace gas remote sensing [19]. The observation of the TOC was of importance due to the frequent occurrence of ozone holes (TOC < 220 DU) over Antarctica with harmful effects on the biosphere. Since water vapor is the most important greenhouse gas in the atmosphere, accurate PWV information is essential. Heat radiated from the Earth’s surface is absorbed by water vapor molecules in the lower atmosphere. The water vapor molecules in turn radiate heat in all directions. Some of the heat returns to the underlying surface. Thus, water vapor is a second source of warmth (in addition to sunlight) at the Antarctic surface. The atmosphere in inner Antarctica is very dry with values of PWV around 0.4–0.5 mm on average [20,21,22].

2. EnMAP

EnMAP was launched on 1 April 2022 and is scheduled to operate until 2027. The swath is 30 km and the spatial resolution is 30 m. The nadir re-visit time is 27 days. An off-nadir (30 degrees) pointing feature provides an opportunity for more frequent revisit times (once every four days). EnMAP contains two detectors that have a wavelength overlap:

• Visible–near infrared (VNIR)—operating in the range of 418 nm to 993 nm;
• Short-wave infrared (SWIR)—operating in the range of 902 nm to 2445 nm.

There exist three levels of EnMAP data products that can be ordered by the users:

• L1B: top-of-atmosphere radiance; the Level 1B processor converts digital numbers into calibrated at-sensor radiances and applies dark signal, non-linearity, gain-matching, response non-uniformity and straylight corrections as well as radiometric calibration.
• L1C: orthorectified L1B data; the Level 1C processor carries out a direct georeferencing of the L1B product, which accounts for sensor-, satellite-motion-, and terrain-related geometric distortions. The L1C product is an orthorectified single data cube that is resampled and transformed to a map projection system (e.g., the UTM Universal Transverse Mercator projection with WGS84 datum).
• L2A: bottom-of-atmosphere reflectance (atmospherically corrected L1C data) and in case of water surfaces optionally normalized water-leaving reflectance or underwater reflectance. It also includes the correction for thin cirrus, haze, terrain, and adjacency effects. The land product is fully compliant with the CEOS CARD4L guidelines [23].

The spectral ranges of 1400–1450 nm and 1770–1930 nm are strongly affected by water vapor and are therefore masked out. The full width at half maximum (FWHM) of the instrument response function depends on the channel. It is in the range of 5.8–10.8 nm for the VNIR and 7.1–11.5 nm for the SWIR detectors.

The features of the instrument, which are of importance for snow research, are the capability to make the measurements over the polar regions and also the coverage of the spectral regions with the maximum sensitivity to the snow parameters (e.g., snow grain size). Further details on the instrument are given in [4,16].

This study was based on TOA reflectance data in sensor and map geometry that is not provided by the ENMAP L1B and L1C products directly. However, it can be converted using the following equation: $R=\pi I/sF{\mu }_0$, where R is the TOA reflectance, I is the TOA radiance, ${\mu }_0$ is the cosine of the solar zenith angle (SZA), F is the extraterrestrial solar flux at the surface perpendicular to the solar beam, and s is the parameter close to 1, which accounts for the Earth–Sun distance variation. The value of F was not measured by the EnMAP directly and was convolved from the solar irradiance model presented in [24].

3. Theory

3.1. Retrieval of Snow Properties

The snow grain size is an important parameter of snow surfaces. A change in grain size at the ground may serve as an indication of recent snowfalls, melting events, or blowing snow drifts. This is related to the fact that snow covers with larger snow grains are characterized by larger absorption (smaller reflection) of electromagnetic radiation in the SWIR spectral range and thus control the amount of solar radiation that is reflected back to outer space. As a matter of fact, the following relationship between the snow reflectance $R_s$ and the effective absorption length (EAL) L in the VNIR and SWIR regions of the electromagnetic spectrum exists [3,25]:

$$R_s=R_0{r}_s^f,$$

(1)

where:

$$r_s=exp\left(-\sqrt{\alpha L}\right)$$

(2)

is the spherical albedo, $\alpha =4\pi \chi /\lambda$ is the bulk ice absorption coefficient at the wavelength λ, $\chi$ is the imaginary part of the ice refractive index [26,27], and $R_0$ is the value of the snow reflectance at $\alpha =0$. It is assumed that snow is a horizontally and vertically inhomogeneous semi-infinite turbid medium, although this is not always the case for the scale of the EnMAP ground scene. The parameter f is given by the following approximate equation [25]:

$$f=\frac{u\left({\mu }_0\right)u\left(\nu \right)}{R_0},$$

(3)

where $\nu$ is the cosine of the viewing zenith angle (VZA) and ${\mu }_0$ is the cosine of the solar zenith angle (SZA). The escape function $u\left({\mu }_0\right)$ can be parameterized as follows [3]:

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

(4)

The accuracy of Equations (1) and (2) increases in the weak absorption approximation range (the single scattering albedo ${\omega }_0\to 1$) [3,6].

The values of L and $R_0$ can be derived from Equations (1) and (2) analytically using EnMAP measurements at two wavelengths $\left({\lambda }_1,{\lambda }_2\right)$ of the SWIR detector $\left(R_{meas1},R_{meas}{}_2\right)$, where atmospheric light scattering and absorption processes can be neglected. Namely, it follows from Equations (1) and (2) [8] that:

$$R_0=R_{meas,1}^{ϵ}{R}_{meas,2}^{1-ϵ}, L=\frac{1}{{\alpha }_2{f}^2}l{n}^2\left(\frac{R_{meas,2}}{R_0}\right),$$

(5)

where $ϵ=\frac{1}{1-b}, b=\sqrt{\frac{{\alpha }_1}{{\alpha }_2}}$, and ${\alpha }_i=\frac{4\pi {\chi }_i}{{\lambda }_i}.$ Here, $\left({\chi }_1,{\chi }_2\right)$ are imaginary parts of ice absorption coefficients at the wavelengths $\left({\lambda }_1,{\lambda }_2\right)$. We used the following wavelengths of EnMAP free of gaseous absorption: ${\lambda }_1=1026 \mathrm{nm}, {\lambda }_2=1235 \mathrm{nm}$. Then it follows that:

$${\chi }_1=2.298\times {10}^{-6},{\chi }_2=1.178\times {10}^{-5},$$

(6)

where we used the linear interpolation of data given in [26]. The effective absorption length (EAL) L can be used to calculate multiple snow properties including the broadband albedo [28] and BRDF (see Table 1). In particular, the effective grain diameter (EGD) d can be found using the following linear relationship [8]:

$$L=\kappa d,$$

(7)

where the constant $\kappa$ depends on the snow type and the shape of the grains. The effective grain diameter is defined as [29]:

$$d=\frac{3V}{2\Sigma },$$

(8)

where V is the average volume of grains and Σ is their average projected area [29]. We shall use the following value for the constant in Equation (7): $\kappa =16$ [8,30,31,32].

One can also derive the snow specific surface area (SSA) σ [33] from the value of EAL. Namely, it follows that:

$$\sigma =\frac{S}{M},$$

(9)

where S is the average surface area and M is the average mass of snow grains. In addition, one can derive the following from Equations (8) and (9):

$$\sigma =\frac{\xi }{\rho d},$$

(10)

where ξ = 3S/2Σ. Therefore, it follows that:

$$\sigma =\frac{C}{\rho L},$$

(11)

where:

$$C=\frac{3\kappa S}{2\Sigma }.$$

(12)

The calibration constant C depends on the shape of grains and snow type. It does not depend on the size of snow grains and can be derived from independent measurements of the SSA (say, using microtomography or adsorption techniques [34] and EAL using, e.g., snow albedo measurements. In particular, it follows from Equation (2) that:

$$L=\frac{l{n}^2{r}_s\left({\lambda }_{swir}\right)}{\alpha \left({\lambda }_{swir}\right)},$$

(13)

where ${\lambda }_{swir}$ is the wavelength in the SWIR region of the electromagnetic spectrum (e.g., 1235 nm). It follows for the surface area of randomly oriented convex particles [35] that: S =$\gamma \Sigma$, where $\gamma =4$ and therefore C = 6$\kappa .$ We used this assumption in the retrieval of the SSA based on spaceborne observations (see Table 1). It should be pointed out that concave particles exist in snow as well [3]. In particular, it follows that $\gamma =7.35$ for the so-called concave Koch snowflakes [29]. In this case, the constant C is almost two times larger as compared to the assumption used in this paper. Therefore, the surface area is substantially increased for concave particles having the same parameter Σ as for convex particles. The expected range of the parameter C is between 90 and 180, and therefore more accurate estimations of the calibration constant C based on experimental measurements are needed.

The bottom-of-atmosphere snow spectral reflectance and snow spherical/plane (spectral and broadband) albedo can be easily derived if the value of L is known [3] (see Table 1 and

Table 2. The parameters of the spherical broadband albedo (sBBA) parameterization for different spectral regions. The parameters for the plane broadband albedo (pBBA) parameterization are the same except p is substituted by $u_{}^2\left({\mu }_0\right)p$[36].

BBA	a	b	$$p,{\mathrm{mm}}^{-1}$$
Spectral range: 0.3–0.7 μm	0	1	7.86 × 10−5
Spectral range: 0.7–2.5 μm	0.2335	0.6600	3.27 × 10−2
Spectral range: 0.3–0.25 μm	0.5721	0.3612	2.35 × 10−2

). The snow bidirectional reflectance distribution function for given observation and illumination conditions can be also derived. It is related to the reflectance R by a simple equation: BRDF = R/π. We concluded that the EAL is the major parameter of snow covers and must be reported in current and future snow remote sensing algorithms based on the measurements in the SWIR range of the electromagnetic spectrum.

Table 1. The EnMAP clean snow products [3,15,36]. The plane albedo is defined as an integral of the azimuthally averaged reflection function with respect to the viewing zenith angle. The spherical albedo is the integral of the plane albedo with respect to the solar zenith angle.

Satellite Snow Product	Abbreviation/Units	Equation	Comments
Effective absorption length	EAL, mm	L
Effective grain diameter	EGD, mm	d = L/κ	κ = 16
Specific surface area	SSA, ${\mathrm{m}}^2/\mathrm{kg}$	σ = 6κ/ρL	ρ = 0.917 g${\mathrm{cm}}^{-3}$ is the density of bulk ice
Broadband albedo Plane BBA Spherical BBA	BBA pBBA sBBA	$r=a+b{e}^{-\sqrt{pL}}$	The values of a, b, and p depend on the spectral ranges used to compute the BBA (see Table 2)
Spherical albedo	SA	$r_s=e^{-\sqrt{\alpha L}}$	$\alpha$ is the bulk ice absorption coefficient
Plane albedo	PA	$r_p=e^{-u\left({\mu }_0\right)\sqrt{\alpha L}}$	$u\left({\mu }_0\right)$ is given by Equation (4)
BOA reflectance	BOAR	$R_{BOA}=R_0{e}^{-f\sqrt{\alpha L}}$	$f$ is given by Equation (3)

Table 1. The EnMAP clean snow products [3,15,36]. The plane albedo is defined as an integral of the azimuthally averaged reflection function with respect to the viewing zenith angle. The spherical albedo is the integral of the plane albedo with respect to the solar zenith angle.

Satellite Snow Product	Abbreviation/Units	Equation	Comments
Effective absorption length	EAL, mm	L
Effective grain diameter	EGD, mm	d = L/κ	κ = 16
Specific surface area	SSA, ${\mathrm{m}}^2/\mathrm{kg}$	σ = 6κ/ρL	ρ = 0.917 g${\mathrm{cm}}^{-3}$ is the density of bulk ice
Broadband albedo Plane BBA Spherical BBA	BBA pBBA sBBA	$r=a+b{e}^{-\sqrt{pL}}$	The values of a, b, and p depend on the spectral ranges used to compute the BBA (see Table 2)
Spherical albedo	SA	$r_s=e^{-\sqrt{\alpha L}}$	$\alpha$ is the bulk ice absorption coefficient
Plane albedo	PA	$r_p=e^{-u\left({\mu }_0\right)\sqrt{\alpha L}}$	$u\left({\mu }_0\right)$ is given by Equation (4)
BOA reflectance	BOAR	$R_{BOA}=R_0{e}^{-f\sqrt{\alpha L}}$	$f$ is given by Equation (3)

3.2. Retrieval of Precipitable Water Vapor and Total Ozone Column

The differences between the retrieved bottom-of-atmosphere reflectance (see Table 1) and the EnMAP-measured top-of-atmosphere reflectance can be attributed to atmospheric scattering and absorption processes. Therefore, there is a way to obtain atmospheric properties from EnMAP measurements.

To demonstrate this statement, we applied the retrieval procedure described above to the EnMAP measurements at Dome C in the vicinity of the Concordia station (75.1°S, 123.35°E) in Antarctica (EnMAP data take ID: 4946; acquisition datetime: 2022-10-29 00:11:38 (UTC+0); tile: 002; processing version: 010111; processing date: 13 December 2022; off-nadir pointing (across-track): 13.84 degrees; off-nadir pointing (along-track): 0.16 degrees; solar zenith angle: 67.26 degrees; solar azimuth angle: 53.31 degrees).

The EnMAP L1B top-of-atmosphere radiance product in sensor geometry data was used. The data cubes were scaled band-wise based on the gain and offset information provided in the metadata to mW/${\mathrm{m}}^2$/sr/nm. Slight EnMAP inherent across-track radiometric misregistrations, which result in striping artifacts in the L1B data product and the derivatives based on them, were compensated using an in-house de-striping algorithm (which will be implemented in the official data product in the future). The de-striped TOA radiance data was then converted into TOA reflectance as discussed above. All derived snow parameters in this paper are based on the TOA reflectance in sensor geometry generated based on the de-striped L1B data set.

A comparison of the derived spectral BOA reflectance and EnMAP-measured TOA reflectance is shown in Figure 1. The calculation of BOA spectral reflectance was based on two retrieved snow parameters (L = 2.3163 mm and $R_0=0.9534$) and the bulk ice spectral absorption coefficient provided in [26] (in the near infrared and SWIR) and in [27] (in the visible). As underlined in the previous section, the SWIR channels located at $1026 \mathrm{nm} \mathrm{and} 1235 \mathrm{nm}$ were used. As can be seen from the analysis of the EnMAP TOA reflectance given in Figure 1, the measurements at these channels were not affected by gaseous absorption and atmospheric light scattering effects in an optically thin atmosphere at Dome C [37]. Actually, BOA and TOA reflectances at these channels almost coincided, which was used in this work to generate snow properties from EnMAP measurements.

The difference in BOA and TOA reflectances around 450 nm (where atmospheric light absorption processes could be neglected) was very small, which highlighted the fact that this difference could be attributed to molecular scattering processes with negligible contribution from aerosol particles present at very small concentrations at this highly elevated Antarctic site.

The difference in curves shown in Figure 1 in the spectral range of 400–1300 nm was mostly due to gaseous absorption. The depths of gaseous absorption bands (see, e.g., the ozone absorption band at 600 nm in Figure 1) could be used to retrieve their total columns from EnMAP reflectance spectra [19]. The differences in the range of reflectance below 0.3—where BOA and TOA reflectance must almost coincide for the optically thin Antarctic atmosphere (outside gaseous absorption bands)—were mostly due to the limited-range applicability of the theory described in the previous section. The presented theory is valid only for weakly absorbing, strongly light-scattering media [6,25]. In addition, the layered nature of snowpack may play a role [38,39,40].

We proposed to use the measurements at the wavelength ${\lambda }_3=1128.47 \mathrm{nm}$ to estimate the precipitable water vapor from EnMAP data. As one can see in Figure 1, the TOA atmosphere reflectance just outside of the water absorption band located at 1130 nm was well described by the derived BOA reflectance. Therefore, we can propose the following simple analytical model for the TOA reflectance in the studied water absorption band:

$$R_{TOA}\left(\lambda \right)=R_{TOA,nogas}\left(\lambda \right)T_{H_{2}O}\left(\lambda \right),$$

(14)

where $R_{TOA,nogas}\left(\lambda \right)$ can be approximated by the snow reflection function $R_s\left(\lambda \right)$ given by Equations (1) and (2) and the pair ($L,R_0$) is derived as explained above. Equation (14) makes it possible to determine the water vapor spectral transmission function:

$$T_{H_{2}O}\left(\lambda \right)=\frac{R_{TOA}\left(\lambda \right)}{R_s\left(\lambda \right)}.$$

(15)

We can also introduce the slant water vapor optical thickness ${\tau }_{H_{2}O}\left(\lambda \right)$ such that:

$$T_{H_{2}O}\left(\lambda \right)=exp\left(-{\tau }_{H_{2}O}\left(\lambda \right)\right)$$

(16)

or:

$${\tau }_{H_{2}O}\left(\lambda \right)=ln\left(T_{H_{2}O}^{-1}\left(\lambda \right)\right)$$

(17)

One can see that the task of the precipitable water vapor (PWV) determination from TOA EnMAP spectra was reduced to the determination of the precipitable water vapor from the spectrum ${\tau }_{H_{2}O}\left(\lambda \right)$. The spectrum ${\tau }_{H_{2}O}\left(\lambda \right)$ can be simulated for various illumination and observation geometries as well for air pressure, temperature, and water vapor profiles using suitable radiative transfer models. Corresponding look-up-tables can be used to derive the value of PWV from the analysis of the spectrum ${\tau }_{H_{2}O}\left(\lambda \right)$.

In this work, we used the approximation for the slant water vapor optical thickness proposed in [41,42,43]:

$${\tau }_{H_{2}O}\left(\lambda \right)=x^n\left(\lambda \right),$$

(18)

where n = 0.646 for the water vapor absorption band located at 1130 nm [41] and:

$$x\left(\lambda \right)=BM{N}_{H_{2}O}k\left(\lambda \right).$$

(19)

Here, $N_{H_{2}O}$ is the precipitable water vapor column to be determined, $k\left(\lambda \right)$ is the absorption coefficient of water vapor, $M$ is the airmass factor, and B is the correction coefficient, which accounts for the air pressure and temperature profiles. We have used the following expression for the correction coefficient [41]:

$$B={\left(P/P_0\right)}^a{\left(T_0/T\right)}^b,$$

(20)

where a = 0.781, b = 0.439 for the considered absorption band $T_0=273.16 \mathrm{K}, \mathrm{and} P_0=1013.25 \mathrm{hPa}.\mathrm{The} \mathrm{pair}$ (P, T) is the average pressure and temperature for a given location calculated using respective vertical profiles. It was derived in this work using the average profiles for these parameters at Dome C for October (491 hPa and 229 K, respectively). The set was obtained through an average of all the radiosoundings (up to the height H = 7 km; see Figure 3 in [44]) performed by the Antarctic Meteo-Climatological Observatory (AMCO) (https://www.climantartide.it/ (accessed on 6 June 2023)) between 2012 and 2017 while binning the soundings monthly. It should be pointed out that the average ground values of these parameters for October at Dome C at the ground level were equal to 639 hPa and 220 K, respectively. The AMCO reported ground pressure and temperature were equal to 651 hPa and 225 K for the moment of satellite measurements, which were close to the average values. The recorded wind speed was 4.6 m/s, and the relative humidity was 51%. The sky was clean, and clouds were absent.

Equation (19) makes it possible to avoid the use of look-up tables and derive the PWV from the following analytical equation:

$$N_{H_{2}O}=\frac{x\left({\lambda }_3\right)}{BMk\left({\lambda }_3\right)},$$

(21)

where we used the following EnMAP channel: ${\lambda }_3=1128.45 \mathrm{nm}$. It follows from Equations (18) and (21) that:

$$N_{H_{2}O}=\frac{{\tau }_{H_{2}O}^{1/n}\left({\lambda }_3\right)}{BMk\left({\lambda }_3\right)},$$

(22)

where ${\tau }_{H_{2}O}^{}\left({\lambda }_3\right)$ can be obtained from Equations (1), (15), and (17) and $R_{TOA}\left({\lambda }_3\right)=R_{meas}\left({\lambda }_3\right)$. This equation makes it possible to determine the precipitable water vapor for a given air mass factor M and the correction coefficient $B$ using the EnMAP measurements at the wavelength ${\lambda }_3$. We used the following value for the absorption coefficient of water vapor [41]: $k\left({\lambda }_3\right)=1.793 {\mathrm{cm}}^{-1},$ where we accounted for the relatively broad (11 nm) width of the used EnMAP channel.

The same arguments as explained above can be applied to other trace gases, including ozone. In particular, we could derive the total ozone column (TOC) using arguments similar to those used to obtain Equation (22) [31,45]:

$$N_{O_3}=\frac{{\tau }_{O_3}^{}\left({\lambda }_0\right)}{M{C}_{abs, O_3}^{}\left({\lambda }_0\right)},$$

(23)

where:

$${\tau }_{O_3}\left(\lambda \right)=ln\left[\frac{R_{TOA,nogas}\left({\lambda }_0\right)}{R_{TOA}\left({\lambda }_0\right)}\right]$$

(24)

and we assumed that ${\lambda }_0=599.267 \mathrm{nm}$ (close to the center of the ozone Chappuis absorption band). We also used the fact that n = B = 1 for the ozone Chappuis absorption band [43,45]. The spectra of ozone absorption cross-section $C_{abs, O_3}^{}$($\lambda$) measured by Serdyuchenko et al. [46] were used. The value of $C_{abs, O_3}^{}$($\lambda$) was almost not influenced by air temperature and pressure for typical atmospheric conditions in the spectral range studied, which is an important feature as far as practical applications of the described algorithm are concerned. In particular, we used the following cross-section value at T = 213 K [46]:

$$C_{H_{2}O}^{}\left(599.27nm\right)=5.06707\times {10}^{-21}\frac{c{m}^2}{molec}$$

(25)

Then, the value of the total ozone column was expressed in molec/${\mathrm{cm}}^2$. Usually the value of TOC is expressed in Dobson units (DUs). Using the fact that 1 DU = 2.689 × ${10}^{16}\frac{\mathrm{molec}}{{\mathrm{cm}}^2}$, we could derive from Equations (23) and (24) the total ozone column expressed in DU:

$$N_{O_3}=K\frac{{\tau }_{O_3}\left({\lambda }_0\right)}{M},$$

(26)

where K = 7339.26 DU. It should be pointed out that the retrieved ${\tau }_{O_3}^{}\left({\lambda }_0\right)$ and thus the TOC derived using Equation (26) under the assumption that $R_{TOA,nogas}\left({\lambda }_0\right)=R_s\left({\lambda }_0\right)$ was biased because the atmospheric light scattering could not be neglected in this spectral range. This conclusion could be drawn based on the analysis of Figure 1 in the spectral range of 400–700 nm as well. In particular, one can see that the derived BOA reflectance and TOA reflectance differed around 400 nm due to mostly molecular scattering effects. Therefore, we estimated the value of $R_{TOA,nogas}\left({\lambda }_0\right)$ using the third-order polynomial derived from EnMAP measurements at the channels 429.29, 486.94, 706.40, and 839.73 nm. In addition, we used the geometrical approximation for the airmass factor:

$$M=1/\nu +1/{\mu }_0.$$

(27)

This completes the description of our retrieval procedure based entirely on analytical equations. The proposed theory makes it possible to perform fast and accurate retrievals of snow properties and also PWV and TOC over snow using EnMAP measurements.

4. Application of the Retrieval Algorithm to L1B Top-of-Atmosphere EnMAP Radiance Data

We applied the described retrieval algorithm for the whole EnMAP scene (about ${10}^6 \mathrm{pixels}$) in the vicinity of Concordia station at Dome C (75°06′00″S, 123°19′58″E, 3233 m elevation) in Antarctica (29 October 2022, 00:14:38 UTC, SZA = 67.26 degrees, VZA = 13.84 degrees). The position of the EnMAP scene relative to the position of the Concordia station is shown in Figure 2. The station itself was not included in the scene. However, as shown below, the variability in the retrieved characteristics was low at the site studied. Therefore, the results derived at some distance from the station were representative at the station location. In addition, it should be pointed out that the interpretation of the results as derived exactly at the position of the station was complicated by the nearby buildings and the snow disturbance due to various activities at the station.

The results for the snow specific surface area, the effective grain diameter, BBA, TOC, and PWV are shown in Figure 3 and Figure 4.

The histograms of the retrieved parameters closely followed Gaussian distributions $f\left(x\right)=A{e}^{-\frac{1}{2 } {\left(\frac{x-\mu }{\sigma }\right)}^2 }$, where $A=\frac{1}{\sqrt{2\pi }\sigma },$ σ is the standard deviation and $\mu$ is mean shown in

Table 3. The statistical characteristics of the retrieved parameters.

Parameter	Average Value	Standard Deviation	CV (%)
Effective grain diameter, mm	0.1429	0.0078	5.5
pBBA (0.3–2.5 μm)	0.8291	0.0015	0.2
SSA, ${\mathrm{m}}^2$/kg	45.93	2.51	5.5
TOC, DU	193.67	13.98	7.2
PWV, mm	0.172	0.0058	3.4

, where we present average values for the retrieved parameters for the scene studied. The coefficient of variance (CV; average/STDV) was rather small for all retrieved parameters on the scale of the EnMAP scene (30 km × 30 km) at the location studied (see Table 3).

The variation in surface parameters shown in Figure 3 was rather small and could be expected due to variability in the local snow cover characteristics (existence and abundance of sastrugi, their orientation, etc.). One can see that the spatial distributions of BBA and SSA were similar. This was due to the fact that they were linked to the spatial distribution of the effective grain diameter shown in Figure 3. Both the SSA and BBA decreased with the size of ice grains as one might expect.

The spatial distributions of total ozone and water vapor column derived on the scale of 30 m shown in Figure 4 could be partially attributed to the uncertainties of the algorithm and violation of assumptions regarding the underlying surface properties. Therefore, additional validation of the data shown in Figure 4 is needed. It is common for satellite optical instrumentation to deliver total ozone data on the scale of several kilometers and even tens of kilometers, which corresponded just one number as reported in Table 3 (194$\pm 14$ DU). The retrievals as shown in Figure 4 are of additional value at the borders of ozone holes and also in the cases where sharp variations in gaseous emissions at the surface (plumes) were present.

We also derived the spectral plane/spherical albedo and snow reflection function using the algorithm described above. Data for a randomly chosen pixel is presented in Figure 5. One can see that the plane and spherical albedos behaved in expected way by approaching 1.0 in the visible and having several oscillations in the SWIR due to the corresponding bulk ice absorption bands. The results at the wavelengths above 1300 nm were less accurate due to the decreased accuracy of the weak absorption approximation [25] as discussed before.

5. The Comparison with Ground Measurements

We performed a comparison of the data shown in Table 3 with ground measurements performed at Concordia station in Antarctica on 29 October 2022. The ground and satellite measurements were not precisely aligned in space and time. However, as it follows from Table 3, the variations in various parameters at the Dome C Plateau were quite low, and the average parameters presented in Table 3 well represent the underlying snow–atmosphere system. The list of instruments used for validation is given in Table 4.

The spectral plane albedo (PA) was measured at Dome C using the Autosolex instrument [27], an automated visible–near infrared (VNIR) spectrometer that has been acquiring surface albedo measurement every 12 min between October and March since 2012. The Autosolex is composed of two measurement heads located approximately 1.5 m above the surface, each containing two upward- and downward-facing cosine receptors. At such a height, the downward-facing collectors were estimated to receive most of the reflected signal (75%) from the surface from an area of approximately 3.5 m in radius. The receptors were connected to an Ocean Optics Maya2000 Pro spectrometer with fiber optics passing through an optical switch. The Autosolex measured the downward and upward (reflected) solar radiation. Using the calibration and processing steps described in [27], the spectral albedo was calculated as the ratio of the upwelling to downwelling spectral irradiances, providing usable data across the 400–1050 nm range. Owing to the negligible amount of diffuse illumination, the bi-hemispherical reflectance measured by the instrument (here called albedo) was directly compared to the directional–hemispherical reflectance (plane albedo) retrieved from EnMAP images. The spectral albedo measured at the ground by head 1 on 29 October 2022, 03:00 UTC+0 (63 degrees solar zenith angle) is given in Figure 6. The measurements of head 2 were strongly affected by a steep surface just under the sensor during the period and were discarded. The slope under head 1 was much less, but nevertheless due to the low Sun in this season, this slope was sufficient to disturb our measured albedo in the visible range with an unexpected trend and values above 1. These artefacts are documented in [47].

The EnMAP measurements were recorded at a solar zenith angle of 67 degrees; the retrieved plane albedo at Dome C is also shown in the same figure. The selected acquisition of the ground measurements was 3 h later when the Sun was at the highest in the sky in order to minimize the slope effect. EnMAP is less affected by such effects due to the larger spatial resolution (30 m). The coincidence of the satellite and ground-based albedo in the spectral range of 700–1050 nm was excellent. The slight difference in the visible was mostly related to the problems with ground measurements. The small differences in the NIR (say, around 975 nm) could be due to the errors related to the spectral ice refractive index database.

Table 4. The ground optical instrumentations used for the validation of satellite products.

Instrument	Target	Spectral Range	Reference
VNIR spectrometer	Spectral albedo, specific surface area	0.4–1.1 μm	[27]
Pyranometer	Broadband albedo	0.35–2.5 μm	[48]
UV radiometer	Total ozone	0.300, 0.306, 0.310, 0.314, 0.325, 0.338, and 0.364 μm	[49]
FTIR spectrometer	Precipitable water vapor	7.1–100 μm	[50,51,52]

Table 4. The ground optical instrumentations used for the validation of satellite products.

Instrument	Target	Spectral Range	Reference
VNIR spectrometer	Spectral albedo, specific surface area	0.4–1.1 μm	[27]
Pyranometer	Broadband albedo	0.35–2.5 μm	[48]
UV radiometer	Total ozone	0.300, 0.306, 0.310, 0.314, 0.325, 0.338, and 0.364 μm	[49]
FTIR spectrometer	Precipitable water vapor	7.1–100 μm	[50,51,52]

The broadband albedo (BBA) we derived was around 0.83, which was consistent with ground observations of the BBA at Concordia station on 29 October 2022 (see Figure 7a). The slight overestimation of the EnMAP-derived BBA could be due to imperfect temporal and spatial collocation as well as difficulties related to the BBA automatic ground measurements. The albedometer worked in the automatic regime, therefore, the measurements could be affected by sunglint, 3D effects (e.g., sastrugi), and diamond dust. Our results were consistent with other ground and also OLCI/S-3 observations [15,48,53,54,55,56]. In particular, Lanconelli et al. [48] reported that the austral summer albedo of 2007/2008 measured at the ground at Dome C and also at the South Pole were essentially equal (see Figure 7b) and had a value of 0.83$\pm 0.03$, which was similar to our finding (see Table 3). It should be mentioned that the austral summer albedo is quite robust and does not vary considerably from year to year. Because the BBA is determined by the effective grain diameter for clean snow, this indirectly confirmed the derived EGD using EnMAP SWIR measurements.

We also derived the NIR and visible-plane BBA (0.69 and 0.99). These values were close to the earlier reported values (0.66 and 0.96) but at another location and time period in Antarctica [57].

The snow specific area (SSA) derived from the ground spectral albedo measurements reported above (see Figure 6) was 52.7${ \mathrm{m}}^2/\mathrm{kg}$ [58], which was consistent with EnMAP satellite observations at Dome C on the same day (46 ${\mathrm{m}}^2/\mathrm{kg}$, see Figure 8). The difference between the satellite and ground observations of the SSA could be due to different assumptions on the snow grain shape used in the retrieval process as well as a temporal and spatial mismatch of satellite and ground observations. It should be pointed out that the 13% difference between the ground and satellite observations was close to the error of the ground measurements themselves [59]. The validation of the technique for the determination of the SSA at other sites was performed using OLCI/S-3 data [6,7].

The effective ice grain diameter (EGD) we retrieved was 0.14 mm. Therefore, the respective effective radius r = d/2 we retrieved was 0.07 mm, which was consistent with the ground observations at Dome C and also at other places in Antarctica [38,60]. In particular, Gay et al. [60] measured the value of the grain radius in the range of 0.05–0.18 mm at the top layer of the snow at Dome C. The derived value of the grain size was closer to the lower border of this interval, although one should account for the differences in the definition of the effective grain size (see Equation (8)) derived from EnMAP observations and the grain size measured on the ground using digital image processing, where the mean radius of all the convex parts of the ice/air boundaries or the mean convex radius was measured. Only for spherical particles did both sizes coincide. Grenfell et al. [38] measured the values of the snow grain radius in the range of 0.05–0.15 mm at the Vostok station in Antarctica. These measurements were consistent with those performed at Dome C and confirmed small spatial variations in the grain size in central Antarctica [60].

The total ozone column (TOC) was measured on the ground as well. The UV-RAD ultraviolet radiometer was used for this purpose. The instrument was designed to monitor solar UV radiance and the ozone column in the polar regions [61]. The instrument measures UV irradiance with seven channels equipped with narrow-band interference filters with a width at half maximum of about 1 nm. In addition, two optical passbands elements attenuate the wings of the filters, and the seven channels are centered at 300, 306, 310, 314, 325, 338, and 364 nm. The filters are mounted on a rotating wheel that is contained inside an isolated and thermostatted box together with the photomultiplier used as a sensor. The estimation of the ozone column from the UV-RAD measurements was performed by using the Stamnes method developed for the case of global irradiance measurements [62]. This approach compares the ratio of irradiances measured at two wavelengths with the corresponding ratio calculated for different solar zenith angles and ozone columns using radiative transfer modeling. One of the wavelengths is chosen to be among the spectral band of stronger ozone absorption, while the other has to be negligibly absorbed. For the operation at Concordia, 314 and 338 nm were selected as the first and second wavelengths, and the lookup table with the model estimations of the corresponding irradiance ratios was composed using Tropospheric Ultraviolet–Visible (TUV) computer code [63]. The scanning of all UV-RAD channels usually took 100–120 s depending on the weather conditions, and the radiometer performed such scanning at 5 min intervals. The irradiance at 338 nm was additionally measured after the 310 nm channel to reduce the time between wavelengths selected for the ozone calculation. The ozone column was determined from the measurements carried out for SZAs below 75°, and the values obtained within the day were smoothed using an averaging procedure. The root mean square (RMS) of deviations from the smoothed curve was considered as an estimate of the measurement error during the day. In clear sky conditions, such an error was found to be between 2 and 5 DU, while in the presence of broken clouds it could increase up to 10 DU. The last comparison of the UV-RAD with a similar radiometer was performed in summer 2019, and the stability of the instrument was continuously controlled by comparing the measured model irradiances and those computed using TUV. Figure 9 illustrates the temporal ozone evolution for 28–30 October 2022. The EnMAP-retrieved TOC (193.7 DU) is shown by a red point that was close to the ground measurements. The average values of the TOC derived from ground measurements were 187 DU for 28 October and 183 DU for 29 October (a 3.5 and 5.8% difference, respectively, from the EnMAP-derived value). A comprehensive validation of the satellite TOC estimations using OLCI/S-3 over snow in the Chappuis absorption band was conducted in [45]. It was found that the average bias of total ozone OLCI retrievals was below 3%. We expected the same accuracy for the EnMAP 30 m measurements and plan to validate our EnMAP TOC retrieval technique using a larger dataset in our next publication.

The total precipitable water vapor (PWV) in the Dome C region, East Antarctic Plateau, was obtained through the analysis of the infrared atmospheric radiance spectra measured by the REFIR-PAD instrument. The results of measurements are given in Figure 10. The Radiation Explorer in the Far-Infrared–Prototype for Applications And Development (REFIR-PAD) [50,51,52] is a Fourier-transform spectroradiometer designed to measure the emitted atmospheric radiance in the 100–1500 cm⁻¹ (6.67–100 µm) spectral range with a spectral resolution up to 0.25 cm⁻¹, and it was successfully deployed in the 2006–2011 period in several campaigns (both in ground-based zenith-looking and stratospheric balloon-borne nadir-looking configurations). In December 2011, the REFIR-PAD instrument was permanently installed at Concordia Station, Antarctica, where it performs a continuous monitoring of the downwelling infrared radiation with a nominal resolution of 0.4 cm⁻¹ and a measurement repetition rate of about 12 min. To retrieve the vertical profiles of temperature and water vapor, the measured spectra were processed through an inversion code based on the LBLRTM atmospheric forward model described in [64]. The retrieval was performed on a 21-level logarithmically spaced vertical grid covering the range between the ground and 23 km. The grid provided a layer thickness that increased with altitude to compensate for the decrease in vertical resolution, which is an inherent property of the retrieval in the zenith-looking observation geometry. The number of fitted levels was four for water vapor (3 m, 170 m, 1 km, and 6 km) and six for temperature (ground, 3 m, 30 m, 170 m, 1 km, 6 km). The PWV values corresponding to the retrieved vertical profiles were provided as a direct output by the LBLRTM forward model; the uncertainty on the PWV was estimated through the fitting errors provided by the MINUIT routine and was about 5%. The PWV retrieval process was successfully validated through comparison with a microwave radiometer, a Raman lidar, and radiosondes [51]. It followed that the retrieved PWV was 0.20 mm for the case studied in this work (0.03 mm or 15% difference from the value we derived (0.17 mm at 00:14 UTC on 29 October 2022).

These results were very encouraging and point to the fact that reliable algorithms for the determination of the TOC and PWV derived from EnMAP data can be developed in the future.

6. Conclusions

We performed for the first time retrievals of the atmosphere–underlying snow surface system parameters using the EnMAP observations over an area close to the Concordia Station in Antarctica. The retrieval of the snow grain size was based on measurements at two EnMAP SWIR channels (1026 nm and 1235 nm) for which the influence of atmospheric parameters on the TOA reflectance could be neglected. This made it possible to determine the spectral bottom-of-atmosphere snow reflectance as well as the snow specific surface area and broadband albedo in the framework of the weak absorption approximation of the radiative transfer theory [25,31,65] under the assumption that the snow was a semi-infinite layer that contained no impurities and was homogeneous in the vertical and horizontal directions. The manual cloud screening procedure was used.

Clearly, the deviations in the derived BOA reflectance from the EnMAP TOA reflectance in the SWIR were mostly due to light absorption by atmospheric trace gases. This assumption made it possible to estimate the total precipitable water vapor (PWV) from the EnMAP TOA reflectance at the wavelength of 1128.5 nm over snow. The total ozone column (TOC) was derived using the EnMAP measurements at a wavelength of 600 nm following the procedure described in [45] for OLCI/S-3. We did not perform the estimation of aerosol characteristics from EnMAP measurements because the BOA and TOA reflectances in the visible outside ozone absorption band were very close (see Figure 1) and the differences originated mostly due to molecular light scattering, pointing to the fact that aerosol optical thickness in the visible was very low during the EnMAP measurements (0.02 or even lower at 550 nm). The preliminary validation of the derived results was performed with independent ground measurements, and satisfactory close correspondence was found between the derived and retrieved parameters. A detailed validation of the proposed retrieval algorithm will be presented in a future work.

Overall, it was shown that EnMAP measurements can be used to assess the environmental conditions around Concordia station in Antarctica. The parameters of snow cover such as the snow grain size and snow albedo were derived from EnMAP measurements. It was found that the diameter of ice grains was quite small (0.14 mm), which was consistent with independent ground measurements [60]. The derived broadband snow albedo was equal to 0.83, which was also consistent with local and long-term ground measurements [48]. The retrieval technique was very fast and based on the analytical solution of the inverse problem described in this paper and elsewhere [15]. The technique can be applied to multiple current and future hyperspectral instrumentation positioned at the bottom of the atmosphere, the top of the atmosphere, and also within the terrestrial atmosphere.

The EnMAP measurements confirmed that the environment around Concordia is very dry with low aerosol and total ozone loads (as known, the ozone hole occurred in October 2022 over Concordia). The average value of the total ozone derived using EnMAP measurements (194 $\pm$14 DU) was close to that measured on the ground (183 DU). The derived PWV we found (0.17$\pm$ 0.01 mm) was close to the ground measurements (0.2 mm) as well. The coefficients of variance of both the TOC and PWV were very small for the region studied. Therefore, the derived TOC and PWV spatial distributions maps were influenced by the noise related to the assumptions used in the retrieval process and the underlying surface variability and require further investigations.

Although the derived parameters were known with a high precision for the Concordia station due to efforts related to the maintenance of ground-based instrumentation, this is not the case for other places in Antarctica, where ground measurements are absent and EnMAP measurements can be used for the assessment of key environmental parameters as shown in this paper. A similar spaceborne instrumentation can be also used for surveys of the atmosphere and underlying surface of other planets.