The Approximate Analytical Solution for the Top-of-Atmosphere Spectral Reflectance of Atmosphere—Underlying Snow System over Antarctica

Abstract

The analytical solutions of the radiative transfer equation are needed for the solution of various applied atmospheric and snow optics problems. In this paper, we propose a simple analytical equation for the top-of-atmosphere (TOA) spectral reflectance. To simplify the problem under study we consider the case of Antarctica, where both snow and atmosphere are almost free of pollutants. This work is focused on the simulation of the moderate spectral resolution TOA measurements (1 nm or so) and the spectral range 400–1000 nm. The values of the coefficient of variance (CV) between the measured by the Ocean and Land Colour Instrument (OLCI) on board Sentinel-3A and modelled spectra are smaller than 10% for most cases in Antarctica. There are regions in Eastern Antarctica, where the values of CV are smaller than 5%. The areas with larger deviations between measured and retrieved spectra could be due to the presence of clouds or structures on the snow surface not captured by the proposed model.

1. Introduction

The optical signals detected on multiple satellite platforms over snow surfaces are determined by the optical properties of the snow surface and atmosphere. The solution of both direct and inverse problems of an atmosphere—underlying snow system requires simple relationships between top-of-atmosphere (TOA) reflectance R and the microphysical/optical characteristics of both snow and atmosphere. The task of this paper is to present a simple analytical relationship between the value of R as detected on a satellite with atmosphere/snow properties. Such a relationship can be established using a numerical solution of an integro-differential radiative transfer equation (RTE) [1]. However, this path is quite complicated and time-consuming. The analytical solutions of RTE are needed for the solution of various applied atmospheric and snow optics problems [2,3,4,5]. This is the main driver of this work. To simplify the problem under study we consider the case of Antarctica, where both snow and atmosphere are almost free of pollutants. This work is focused on the simulation of the moderate spectral resolution TOA measurements (1 nm or so) and the spectral range 400–1000 nm.

In the next section, we present the TOA reflectance model. Section 3 is devoted to the application of the model to the measurements performed by the Ocean and Land Colour Instrument (OLCI) on board Sentinel-3.

2. The Top-of-Atmosphere Reflectance Model

2.1. The Reflection of Light from Atmosphere—Underlying Surface System

The top-of atmosphere spectral reflectance over a Lambertian surface with albedo $r_s$ can be presented in the following way [1,6]:

$$R=\left(R_a+\frac{T_a{r}_s{}_{}}{1-r_a{r}_s}\right)T_g.$$

(1)

Here, the values of $R_a, r_a, \mathrm{and} T_a$represent the reflectance, albedo, and total transmittance, respectively, for the artificial atmosphere without gaseous absorbers and black underlying surface, $T_g$ is the gaseous transmittance. The snow is close to the Lambertian surface in the visible. However, there are differences. Therefore, to account for the non-Lambertian effects we substitute the snow albedo in the nominator of Equation (1) by the snow spectral reflectance $R_s$. Then it follows [7]:

$$R=\left(R_a+\gamma T_a{R}_s \right)T_g,$$

(2)

where $\gamma ={\left(1-r_a{r}_s\right)}^{-1}$. It follows from Equation (2) that $R\to R_s$at atmospheric optical thickness $\tau \to 0$ outside gaseous absorption bands as it should be. This is because $R_a\to 0,r_a\to 0$ and $T_a\to 1, T_g=1$ in this case by the definition. The errors of Equation (2) increase with the deviation of the underlying surface from the Lambertian one and also with the increase in optical thickness $\tau$. Then, yet another approach for the account of underlying surface effects on the value of TOA reflectance must be used (see, e.g., [8]).

However, the focus of this paper is the simulation of satellite optical signals over Antarctica, where $\tau$ is small and the error of this approximation is acceptable for the solution of most practical problems related, e.g., to snow remote sensing from space [7].

One can see that the desired analytical solution to the problem at hand can be obtained if analytical equations for all functions in Equation (2) are known. The analytical solutions for the corresponding functions are considered in the next section.

2.2. The Atmospheric Path Reflectance

The atmospheric path reflectance $R_a$can be presented in the following way using the Sobolev approximation [7,9,10]:

$$R_a=R_{ss}+R_{ms},$$

(3)

where the single scattering contribution can be presented as

$$R_{ss}=M\left(\tau \right)p\left(\theta \right),$$

(4)

$$M\left(\tau \right)=\frac{1-\exp {\left(-m\tau \right)}^{}}{4\left({\mu }_0+\mu \right)},$$

(5)

$$m={\mu }_0^{-1}+{\mu }^{-1}.$$

(6)

Here, ${\mu }_0$ is the cosine of the solar zenith angle (SZA), $\mu$is the cosine of the viewing zenith angle (VZA), $\theta$ is the scattering angle defined as

$$\cos \theta =-{\mu }_0\mu +s_{0}s\cos \phi ,$$

(7)

$\phi$ is the relative azimuthal angle (RAA),$s_0$ is the sine of the SZA, $s$is the sine of the VZA, $\tau$is the atmospheric optical thickness, $p\left(\theta \right)$ is the phase function.

The multiple light scattering contribution$R_{ms}$ can be approximated as [9]:

$$R_{ms}=1+M\left(\tau \right)q\left({\mu }_0,\mu \right)-\frac{N\left(\tau \right)}{1+0.75\left(1-g\right)\tau },$$

(8)

where

$$N\left(\tau \right)=f\left({\mu }_0\right)f\left(\mu \right),$$

(9)

$$f\left({\mu }_0\right)=\frac{1}{2}\left\{1+\frac{3}{2}{\mu }_0+\left(1-\frac{3}{2}{\mu }_0\right)\exp {\left(-\frac{\tau }{{\mu }_0}\right)}^{}\right\},$$

(10)

$$q\left({\mu }_0,\mu \right)=3\left(1+g\right){\mu }_0\mu -2\left({\mu }_0+\mu \right).$$

(11)

The value of the average cosine of scattering angle g is determined by the following expression:

$$g=\frac{1}{2}{{\displaystyle \int }}_0^{\pi }p\left(\theta \right)\mathit{sin}\theta \mathit{cos}\theta \mathit{d}\theta .$$

(12)

Equation (8) can be substituted by other approximations or parameterisations [3,11] if needed.

Equations (3), (4) and (8) make it possible to calculate the path reflectance analytically for a given phase function $p\left(\theta \right)$and atmospheric optical thickness $\tau$. The approximate account for aerosol absorption effects can be performed multiplying $R_{ss}$ by the single scattering albedo ${\omega }_0.$It follows from Equation (8) that $R_{ms}\to 0$ as $\tau \to 0$ as it should be. The accuracy of the Sobolev approximation has been thoroughly studied [10,12,13]. It has been found that the error of this approximation is smaller than 10% for most of the observation geometries at SZA and VZA smaller than 75 degrees, and an atmospheric optical thickness smaller than 0.5, which is the case for an atmospheric layer above Antarctica.

The accuracy of the Sobolev approximation presented here increases as $g\to 0,$ which is the case for the clean atmosphere of Antarctica. Indeed, the atmospheric optical thickness outside gaseous absorption bands can be presented as

$$\tau \left(\lambda \right)={\tau }_{mol}\left(\lambda \right)+{\tau }_{aer}\left(\lambda \right),$$

(13)

where $\lambda$ is the wavelength, ${\tau }_{mol}$ is the molecular optical thickness, ${\tau }_{aer}$ is the aerosol optical thickness (AOT). Let us take the wavelength of 400 nm. In this case, one can assume that ${\tau }_{mol}=0.23$ [14] and ${\tau }_{aer}=0.02$ [15] at Dome C in Antarctica. One can see that atmospheric molecular scattering and extinction processes are far more important at 400 nm as compared to light scattering by atmospheric aerosol.

The atmospheric phase function for a vertically homogeneous nonabsorbing atmospheric layer can be presented as:

$$p\left(\theta \right)=\frac{{\tau }_{mol}{p}_{mol}\left(\theta \right)+{\tau }_{aer}{p}_{aer}\left(\theta \right)}{{\tau }_{mol}+{\tau }_{aer}},$$

(14)

where

$$p_{mol}\left(\theta \right)=\frac{3}{4}\left(1+{cos}^2\theta \right)$$

(15)

is the molecular scattering phase function and $p_{aer}\left(\theta \right)$ is the aerosol phase function. It follows from Equations (12) and (14):

$$g=\frac{{\tau }_{aer}{g}_{aer}}{{\tau }_{mol}+{\tau }_{aer}},$$

(16)

where $g_{aer}$ is the average cosine of scattering angle for atmospheric aerosol. Taking into account that $g_{aer}\approx 0.7$ at 400 nm [7], we have the following estimation: $g\approx 0.06$ at Dome C and the proposed approximation for $R_{ms}$can be used with high accuracy at 400 nm at elevated sites in Antarctica. The value of g somewhat increases at coastal sites being smaller than 0.1–0.2 at 400 nm in most cases. At longer wavelengths, the contribution of multiple scattering decreases and $R_a\to R_{ss}$. Therefore, the errors in the approximation for the value of $R_{ms} \mathrm{become}\mathrm{irrelevant}.$This explains the applicability of the Sobolev approximation derived at small values of g for the studied case of the Antarctic atmosphere. The approximation can be further improved using the truncation approximation as discussed in [10], where also the accuracy of the Sobolev approximation is studied.

We use the following approximation for the spectral aerosol optical thickness [16]:

$${\tau }_{aer}\left(\lambda \right)={\tau }_{aer}\left({\lambda }_0\right){\left(\frac{\lambda }{{\lambda }_0}\right)}^{-\alpha },$$

(17)

where ${\lambda }_0=1 \mathsf{\mu }\mathrm{m}$, and the pair $\alpha$and β = ${\tau }_{aer}\left({\lambda }_0\right)$ represent the Angström parameters. The molecular optical thickness is represented as [17]

$${\tau }_{mol}\left(\lambda \right)=\widehat{P}{\tau }_m\left(\lambda \right),$$

(18)

where $\widehat{P}=\frac{P}{P_0}$, P is the site pressure, $P_0=1013.25 \mathrm{mb}$ and ${\tau }_m\left(\lambda \right)$ isgiven by the approximate formula similar to Equation (17):

$${\tau }_m\left(\lambda \right)={\tau }_m\left({\lambda }_0\right){\left(\frac{\lambda }{{\lambda }_0}\right)}^{-v}.$$

(19)

One can assume that ${\tau }_m\left({\lambda }_0\right)=0.0084$and $\upsilon =4.0932$. Both parameters can vary depending on temperature, humidity and pressure profiles [14,17,18,19]. The site pressure at a given altitude H can be calculated using the following equation: $P=P_0\exp \left(-\frac{H}{H_0}\right)$ assuming the scale height $H_0$ (e.g., 6 km). Alternatively, the value of $\widehat{P}=P/P_0$ can be derived from the information on the atmospheric pressure $P$ at a given location (e.g., from the European Centre for Medium-Range Weather Forecasts (ECMWF) re-analysis (https://www.ecmwf.int/en/forecasts/datasets (accessed on 23 March 2022))). We shall assume that the spectral aerosol asymmetry parameter is described by the following formula [7]:

$$g_{aer}\left(\lambda \right)=g_0+g_1{e}^{-\frac{\lambda }{{\lambda }_0}},$$

(20)

where the wavelength $\lambda$ is in microns and

$$g_0=0.5263,g_1=0.4627,{\lambda }_0=0.4685\mathsf{\mu }\mathrm{m}.$$

(21)

Taking into account that the spectral TOA reflectance is determined mostly by the underlying snow reflectance for the case of a cloudless sky, we can use a simplified approximation for the aerosol phase function $p_{aer}\left(\theta \right)$ in Equation (14). This function can be calculated using Mie theory under the assumption of aerosol particle size distribution and the complex refractive index of aerosol particles or their mixture. In this work we assume that the aerosol phase function can be presented in the following analytical form [20]:

$$p_{aer}\left(\theta \right)=c{p}_{aer,1}\left(\theta \right)+\left(1-c\right)p_{aer,2}\left(\theta \right),$$

(22)

$$\begin{gathered} p_{aer,1}\left(\theta \right)=\frac{1-g_{aer,1}^2}{{\left(1-2g_{aer,1}\cos \theta +g_{aer,1}^2\right)}^{\frac{3}{2}}}, \\ p_{aer,2}\left(\theta \right)=\frac{1-g_{aer,2}^2}{{\left(1-2g_{aer,2}\cos \theta +g_{aer,2}^2\right)}^{\frac{3}{2}}}. \end{gathered}$$

(23)

We shall assume that

$$g_{aer,1}^{}=0.8, g_{aer,2}^{}=-0.45.$$

(24)

It follows for the asymmetry parameter for the case under study:

$$g_{aer}^{}=c{g}_{aer,1}^{}+\left(1-c\right)g_{aer,2}^{}$$

(25)

and, therefore,

$$c=\frac{g_{aer}-g_{aer,2}}{g_{aer,1}-g_{aer,2}},$$

(26)

where we assume that $g_{aer}$ is given by Equation (20). This concludes the description of the approach to derive the atmospheric path reflectance.

2.3. The Atmospheric Spherical Albedo

Namely, the following expression for the atmospheric plane albedo derived using the same technique as used for the determination of the atmospheric path reflectance outlined above can be used [9].

$$r_p\left({\mu }_0\right)=1-\frac{f\left({\mu }_0\right)}{1+0.75(1-g)\tau },$$

(27)

where $f\left({\mu }_0\right)$ is given by Equation (10) and light absorption processes are neglected. It follows from Equation (27) as $\tau$$\to$0:

$$r_p\left({\mu }_0\right)=w\tau ,$$

(28)

where w = ${\left(2{\mu }_0\right)}^{-1}$ − 0.75 g, which is a valid approximation at small values of $g.$ In particular, Equation (28) coincides with a similar formula given in [21,22] for the case of g = 0.

The spherical albedo is defined as

$$r_a=2{\displaystyle \underset{0}{\overset{1}{\int }}{r}_p}\left({\mu }_0\right){\mu }_{0}d{\mu }_0.$$

(29)

The substitution of Equation (27) to Equation (29) gives after integration:

$$r_a=1-\frac{1+\psi \left(\tau \right)}{1+0.75(1-g)\tau },$$

(30)

where

$$\psi (\tau )=\left[1+\frac{\tau }{2}\right]\frac{{\tau }^2}{2}{E}_1\left(\tau \right)-\left[1+\tau \right]\frac{\tau }{4}\exp (-\tau )$$

(31)

and

$$E_1(\tau )={\displaystyle \underset{1}{\overset{\infty }{\int }}\frac{\exp (-x\tau )}{x}dx}$$

(32)

is the exponential integral [23]. The exponential integral has the following properties [20,23]:

$$E_1(\tau \to \infty )=\exp (-\tau )/\tau \to 0,$$

(33)

$$\tau E_1(\tau \to 0)={\tau }^2-\frac{{\tau }^3}{4}+\frac{{\tau }^4}{18}-\frac{{\tau }^5}{96}-\tau \ln \tau -\gamma \tau ,$$

(34)

where $\gamma =0.5772157$ is the Euler’s constant. Therefore, one can see that $r_a$$\to 1$ as τ$\to \infty$ as it should be for nonabsorbing case considered here.

The atmospheric optical thickness is usually smaller than 0.5 for polar atmospheres and one can use the approximation given by Equation (34) for the product $\tau E_1(\tau )$ in Equation (31). Then it follows:

$$\psi (\tau )=\frac{1}{2}\left(\tau +\frac{{\tau }^2}{2}\right)\left({\tau }^2-\frac{{\tau }^3}{4}+\frac{{\tau }^4}{18}-\frac{{\tau }^5}{96}-\tau \ln \tau -\gamma \tau \right)-\frac{1}{4}\left(1+\tau \right)\tau \exp (-\tau ).$$

(35)

Therefore, one derives:

$$\begin{gathered} r_a=1- \\ \left[\frac{1-\frac{\tau }{4}\left(1+\tau \right)\exp (-\tau )+\frac{{\tau }^2}{2}\left(1+\frac{{\tau }^{}}{2}\right)\left({\tau }^{}-\frac{{\tau }^2}{4}+\frac{{\tau }^3}{18}-\frac{{\tau }^4}{96}-\ln \tau -\gamma \right)}{1+0.75(1-g)\tau }\right]. \end{gathered}$$

(36)

It follows from Equation (36) as $\tau$$\to 0:$$r_a\to b\tau ,$ where $b=1-3g/4$. One can see that $r_a\approx \tau$ at small values of atmospheric optical thickness and g = 0. The value of $r_a<\tau$ at g > 0.

The intercomparison of the dependence of spherical albedo on the value of atmospheric optical thickness derived using Equation (36) and the exact radiative transfer calculations is shown in Figure 1, where we also present the results derived using the following simple formula (approximation—2) valid for the case of pure molecular scattering [24]:

$$r_a=\tau (a\exp (-\tau /\alpha )+b\exp (-\tau /\beta )+c)$$

(37)

where a = 0.180, b = 0.583, c = 0.215, $\alpha$ = 0.168, $\beta$ = 1.092. The maximal error of Equation (37) is smaller than 1%. It is based on the parameterisation of exact radiative transfer (RT) calculations for the case of molecular scattering [24]. The error of approximation (36) is smaller than 2% at $\tau \le 1$ for the molecular scattering case. We also show the calculations for the case of aerosol and molecular scattering mixture. It is assumed that aerosol optical thickness is 0.05 and molecular optical thickness changes from 0.01 to 1.0. The aerosol parameters were taken from [25]. One can see that the errors of the approximation do not increase substantially in the presence of a thin aerosol layer. It should be pointed out that one is interested not in the value of $r_a$ but rather in the quantity $\gamma ={(1-r_a{r}_s)}^{-1}$ in the problem under study (see Equation (2)). The error of an approximation for the value of $\gamma$ is even lower.

2.4. The Atmospheric Transmittance Outside Gaseous Absorption Bands

The transmission function $T_a$ is approximated as follows [10]:

$$T_a=t_a^m,$$

(38)

where m is given by Equation (6). The value of $t$is calculated using the following approximation [10]:

$$t_a^{}=e^{-B\tau }.$$

(39)

Here,

$$B=\frac{1}{2}{{\displaystyle \int }}_{\frac{\pi }{2}}^{\pi }p\left(\theta \right)\mathit{sin}\theta d\theta$$

(40)

is the so-called backscattering fraction. The substitution of Equation (14) to Equation (40) gives:

$$B=\frac{{\tau }_{mol}{B}_{mol}+{\tau }_{aer}{B}_{aer}}{{\tau }_{mol}+{\tau }_{aer}}$$

(41)

where $B_{mol}=0.5,$

$$B_{aer}=cB\left(g_{aer,1}\right)+\left(1-c\right)B\left(g_{aer,2}\right)$$

(42)

and

$$B(G)=\frac{1-G}{2G}\left[\frac{1+G}{\sqrt{1+G^2}}-1\right]$$

(43)

It follows at small values of G:

$$B\left(G\right)\approx \frac{1}{2}+\frac{G\left(G^2-3\right)}{2(1+G^2+\left(1-G^2\right)\sqrt{1+G^2)}}.$$

(44)

The accuracy of Equations (38)–(40) has been studied in [10], where it has been found this approximation can be used with an accuracy better than 5% at solar zenith angles smaller than 70 degrees for the continental and water soluble aerosol layers with an aerosol optical thickness smaller than 0.5, which is the case for polar atmospheres.

2.5. The Atmospheric Gaseous Transmittance

Atmospheric gaseous transmittance in Equation (2) can be presented as a product of transmittances due to various gaseous components in the atmosphere. In this work, we are interested primarily in the calculation of radiative characteristics in the range 400–1000 nm, where the main atmospheric gases which contribute significantly to the atmospheric absorption processes, are molecular oxygen, ozone, and water vapor. Therefore, to simplify the calculations, the contribution of other gases is ignored. Such an approach has been also used in several previous papers on the subject [2,11]. In particular, we shall assume that the gaseous transmittance is given by the following product:

$$T_g=T_{O_2}^{}{T}_{O_3}^{}{T}_{H_{2}O}^{},$$

(45)

where $T_{O_2}^{},T_{O_3}^{},\mathrm{and}{T}_{H_{2}O}^{}$ are the gaseous transmittances due to molecular oxygen, ozone and water vapor.

2.5.1. Ozone

We shall assume that the ozone transmittance $T_{O_3}^{}$ can be derived using the following expression [17,27]:

$$T_{O_3}^{}=\exp (-m_{O_3}{\tau }_{O_3}(\lambda )),$$

(46)

where $m_{O_3}$ is the ozone airmass factor (AMF) and ${\tau }_{O_3}(\lambda )$ is the vertical optical density (VOD) of ozone:

$${\tau }_{O_3}(\lambda )={\displaystyle \underset{0}{\overset{z_{top}}{\int }}{C}_{abs,O_3}(\lambda ,z)n_{O_3}(z)}dz,$$

(47)

where $C_{abs,O_3}(\lambda ,z)$ is the ozone absorption cross-section, z is the altitude, $z_{top}$ is the top of atmosphere altitude, $n_{O_3}(z)$ is the number concentration of ozone molecules at the altitude z.

We use the geometrical approximation for the AMF (see Equation (6)). In the spectral region of interest, the dependence of the ozone absorption cross section on both temperature and pressure (and, therefore, the altitude) can be ignored. Then one derives from Equation (47):

$${\tau }_{O_3}(\lambda )=N{C}_{abs,O_3}(\lambda ),$$

(48)

where

$$N_{O_3}={\displaystyle \underset{0}{\overset{z_{top}}{\int }}{n}_j(z)}dz$$

(49)

is the total ozone column (TOC) measured in molecules per ${\mathrm{cm}}^2$. The value of $C_{abs,O_3}$ is measured in ${\mathrm{cm}}^2$ (per molecule) and can be parameterised as follows in the spectral range under study [28]:

$$C_{abs,O3}(\lambda )=AF(\lambda ),$$

(50)

where $A=18.48\times {10}^{-21}{\mathrm{cm}}^2$/molecule

$$F(\lambda )=\frac{\zeta \left(\lambda \right)}{{(1+\zeta \left(\lambda \right))}^2},$$

(51)

$$\zeta =\exp (\frac{w\left(\lambda \right)-w_p}{\Delta }),$$

(52)

$w={10}^7/\lambda$, where $\lambda$ is given in nm and $w$ is the wavenumber with the dimension ${\mathrm{cm}}^{-1}.$ The authors of [28] suggest (after correction of a misprint) the following values for the parameters in Equation (52): $w_p=\mathrm{16,811}{\mathrm{cm}}^{-1},\mathsf{\Delta }=877{\mathrm{cm}}^{-1}(\mathrm{at}\mathrm{wavenumbers}w<w_p)$ and $\mathsf{\Delta }=1210 {\mathrm{cm}}^{-1}$ at $w\ge w_p$. We show the results of calculations of $C_{abs,O3}(\lambda )$ according to this equation and also measurements of $C_{abs,O3}$ [29] in Figure 2a as the function of the wavelength $\lambda ={10}^7/w$, where $w$ is given in 1/cm and the dimension of the wavelength is nm. One can see that theoretical and experimental results differ in the regions, where the spectral curve $C_{abs,O3}(\lambda )$ has oscillations. The relative accuracy is better than 25% in the spectral range 500–700 nm (see Figure 2b), which suits our purposes (especially because we are interested in the value of exp(−$m{\tau }_{O_3}$), which is derived with much better accuracy as compared to the ozone VOD). Summing up, we assume that the value of $T_{O_3}^{}$ can be approximated as follows:

$$T_{O_3}=\exp (-\kappa F\left(\lambda \right)),$$

(53)

where $\kappa =A{m}_{O_3}{N}_{O_3}.$ The error of the ozone transmittance calculation using simple Equation (53) as compared to the case, where experimental values of $C_{abs,O3}$ are used, is shown in Figure 3 at $N_{O3}$ = 300 DU (1 DU = 2.69$\times {10}^{16}$ ozone molecules/${\mathrm{cm}}^2$) and the ozone airmass factor $m_{O_3}$ = 1. One can see that the respective error is below 1% for the case studied, which is smaller than the error of the respective spectral reflectance measurements.

2.5.2. Water Vapor

The transmittance of water vapor in the spectral range 900–1000 nm can be calculated using the following simple equation [28]:

$$T_{H_{2}O}=\exp \left\{-s^{x_1}\right\},$$

(54)

where

$$s=N_{H_{2}O}^{ef}{c}_{H_{2}O}(\lambda ),$$

(55)

$$N_{H_{2}O}^{ef}=Q_{H_{2}O}{m}_{H_{2}O}{N}_{H_{2}O},$$

(56)

$$Q_{H_{2}O}={(\overline{P}/P_0)}^n{(t_0/\overline{t})}^k,$$

(57)

$$c_{H_{2}O}(\lambda )={\displaystyle \sum _{j=1}^2{B}_j{F}_j\left(\lambda \right)},$$

(58)

$$F_j(\lambda )=\frac{{\zeta }_j(\lambda )}{{(1+{\zeta }_j\left(\lambda \right))}^2},{\zeta }_j(\lambda )=\exp (\frac{w(\lambda )-w_j}{{\mathsf{\Delta }}_j}),$$

(59)

$w={10}^7/\lambda$, where $\lambda$ is given in nm and $w$ is the wavenumber with the dimension ${\mathrm{cm}}^{-1},$ n = 0.775, k = 0.721, $t_0$ = 273.16 K, $P_0=1013.25\mathrm{hPa}$, $(\overline{P},\overline{t})$ are average values of pressure and temperature for a given location calculated using respective vertical profiles, $x_1=0.649$; $m_{H_{2}O}$ is the water vapor airmass factor, which can be found using geometrical approximation (see Equation (6)). The constants $B_j$, $w_j$ and ${\mathsf{\Delta }}_j$ are given in

Table 1. The parameters of water vapor absorption model [28].

λ, nm	${\mathit{B}}_{\mathit{j}}$$,{\mathbf{cm}}^{-1}$	${\mathit{w}}_{\mathit{j}}$$,{\mathbf{cm}}^{-1}$	${\Delta }_{\mathit{j}}$$,{\mathbf{cm}}^{-1}$	Comment
910 (j = 1)	0.744	11,099	23.4 $(w<w_j)$ 73.8 $(w\ge w_j)$	Weak absorption band
940 (j = 2)	7.560	10,697	23.1 $(w<w_j)$ 110.2 $(w\ge w_j)$	Strong absorption band

below. The value of the Precipitable Water vapor (PWV) content $N_{H_{2}O}$, which is the input parameter for the model, is given in cm. When measured in linear units (cm, mm), the water vapor column is represented by the height (or depth) $N_{H_{2}O}$ that the water column would occupy if the vapor were condensed into liquid and spread evenly across the column. Using the density of water, we can also report water vapor in ${\mathrm{gcm}}^{-2}$ = 1 cm or ${\mathrm{kgm}}^{-2}$ = 1 mm. The example of water vapor gaseous transmittance calculated in the framework of this model is presented in Figure 4.

Equation (54) is similar to Equation (53) except $x_1$ and Q differ from the unity as assumed in Equation (54). This difference is due to the fact that most of the ozone is located in a single layer positioned in the stratosphere, which is not the case for water vapor with the largest concentrations closer to the surface. The value of the effective PWV content $N_{H_{2}O}^{ef}$ can be estimated from Equations (56) and (57) under the assumption that the value of $Q_{H_{2}O}$ is known. The determination of $N_{H_{2}O}^{}$ requires information on water vapor airmass factor and also temperature and pressure profiles for a given location.

2.5.3. Oxygen

The oxygen transmittance is found using the approximation similar to Equation (54) [30]:

$$T_{O_2}=\exp (-s^{x_2}),$$

(60)

where

$$s=N_{O_2}^{ef}{c}_{O_2}(\lambda ),$$

(61)

$$N_{H_2{O}_2}^{ef}=Q_{O_2}{m}_{O_2}{N}_{O_2}.$$

(62)

The parameter $Q_{O_2}$ is given by formula similar to Equation (57) except n = 0.9353, k = 0.1936, $m_{O_2}$ is the oxygen AMF, which can be estimated using geometrical approximation (see Equation (6)), $N_{O_2}$ is the oxygen amount in atm-cm, $x_2$ = 0.5641.

The spectral function $c_{O_2}(\lambda )$ has been tabulated in [30]. We have parameterised this function using the double Gauss approximation (at the wavelengths smaller or equal to 764 nm):

$$c_{O_2}(\lambda )=\mathsf{\Lambda }\exp (-1.7{(\lambda -{\lambda }_1)}^2)+0.32\mathsf{\Lambda }\exp (-0.7{(\lambda -{\lambda }_2)}^2),$$

(63)

where ${\lambda }_1$ = 760.75 nm, ${\lambda }_2=763.36\mathrm{nm}$, $\mathsf{\Lambda }=1.8\times {10}^{-5}{(\mathrm{cm}-\mathrm{atm})}^{-1}$ and $\lambda$ is expressed in nm. At the wavelengths $\lambda$ larger than 764 nm, the following parameterisation has been used:

$$c_{O_2}(\lambda )=\frac{\Upsilon }{1+\exp \left[(\lambda -{\lambda }_s)/{\Delta }_s\right]},$$

(64)

where $\Upsilon =$ 8.419 $\times {10}^{-6}{(\mathrm{cm}-\mathrm{atm})}^{-1}$, ${\lambda }_s$ = 764.11 nm, ${\mathsf{\Delta }}_s$ = 0.85036 nm. The intercomparison of $c_{O_2}(\lambda )$ derived using Equations (63), (64) and data presented in [30] is shown in Figure 5. It follows that the proposed parameterisations can indeed be used to model gaseous transmittance in the molecular oxygen A-band. The transmittance $T_{O_2}$(λ) at $N_{O_2}$ = 8.706853 $\times {10}^5$ cm-atm and $m_{O_2}=1$ is given in Figure 6 both using Equations (63) and (64), and tabular data for the spectral function given in [30]. Good correspondence between the two approaches for the calculation of $T_{O_2}$(λ) is found.

2.6. Reflectance of Solar Light from a Snow Surface

The only function in Equation (2), which we did not study so far is the spectral snow reflectance $R_s\left(\lambda \right).$ It is modelled using the approximation proposed in [31,32,33,34,35] for the case of a clean vertically and horizontally homogeneous plane—parallel layers of turbid media:

$$R_s(\mu ,{\mu }_0,\phi )=R_0(\mu ,{\mu }_0,\phi )r_s^{\xi },$$

(65)

where

$$\xi =u\left({\mu }_0\right)u\left(\mu \right)/R_0(\mu ,{\mu }_0,\phi ),$$

(66)

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

(67)

$$r_s=\exp \left\{-\sqrt{y}\right\},$$

(68)

$$y=L{\alpha }_i,$$

(69)

$$R_0=\frac{a+b\left({\mu }_0+{\mu }_0\right)+c{\mu }_0\mu +p_s\left(\theta \right)}{4\left({\mu }_0+\mu \right)},$$

(70)

where a = 1.247, b = 1.186, c = 5.157, ${\alpha }_i=\frac{4\pi \chi }{\lambda }\mathrm{is}\mathrm{the}\mathrm{bulk}\mathrm{ice}\mathrm{absorption}\mathrm{coefficient},$ χ is the part of the complex ice refractiveve index [5], and L is the effective absorption length (proportional to the average ice grain size). The phase function of a snow layer$p_s\left(\theta \right)$ is approximated as [5,6]:

$$p_s\left(\theta \right)=11.1\exp \left(-0.087\theta \right)+1.1\exp \left(-0.014\theta \right),$$

(71)

where θ is the scattering angle in degrees. Equation (65) follows from the analysis of the radiative transfer equation valid at the snow single scattering albedo $\omega$close to one, which is a valid approximation in the spectral range under study [5]. Additionally, it is assumed that 1 − ω is proportional to the product of the snow grain size and the bulk ice absorption coefficient and the snow phase function does not change considerably in the visible and near infrared regions of the electromagnetic spectrum [5].

In the case of polluted snow, the equations given below remain valid except the parameter y is changed to the following expression [33,35]:

$$y=L\left[{\alpha }_i+f{(\lambda /{\lambda }_0)}^{-\upsilon }\right],$$

(72)

where $\upsilon$ is the absorption Angström parameter of impurities in snow, ${\lambda }_0$ = 1 μm and f is the parameter, which is proportional to the relative impurity concentration $C=C_{imp}/C_{ice}$, where $C_{imp}$ is the volumetric concentration of impurities, $C_{ice}$ is the volumetric concentration of ice grains. Namely, it follows [35]:

$$f=\frac{C}{K}{\kappa }_0$$

(73)

where $K$ is the absorption enhancement parameter ($\approx 1.8$) and ${\kappa }_0$ is the volumetric absorption coefficient of impurities at the wavelength 1 μm. The value of ${\kappa }_0$ is equal to the ratio of the average absorption cross-section of impurities (say, dust particles) to their average volume [35]. One can see that the description of optical characteristics of polluted snow requires two additional parameters $(f,\upsilon ).$

3. The Intercomparison of Derived Parameterisation with Spaceborne Measurements of Spectral Top-of-Atmosphere Reflectance

The TOA spectral reflectance function calculated using Equation (2) is given in Figure 7. We also show the simulation and measurement results for the top-of-atmosphere OLCI (https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-3-olci (accessed on 18 September 2022)) reflectance at Dome C in Antarctica. The specification of both OLCI channels and their widths is used in the averaging procedure (assuming box function). The derived bottom of atmosphere (BOA) snow reflectance is also shown for comparison. The TOA reflectance measured by OLCI is given by crosses in Figure 7. The assumed parameters of the atmosphere—underlying snow system are given in

Table 2. The parameters of the model used in Figure 7 for the case of clean snow underlying surface.

N	Parameter	Value	Comment
1	α	0.008	atmospheric aerosol
2	β	1.3	atmospheric aerosol
3	$N_{O_2}$, cm-atm	$8.706853\times {10}^4$	oxygen
4	$N_{O_3}$, DU	250	ozone
5	$N_{H{O}_2}$, mm	0.33	water vapor
6	$P$, hPa	650	surface pressure
7	$\overline{P}$, hPa	325	average pressure
8	$\overline{t}$, K	233	average temperature
9	L, mm	2.24	snow effective absorption length

. The molecular optical thickness has been modelled using the following equation:

$${\tau }_{mol}=0.015+20.64\exp (-\lambda /86.16),$$

(74)

where λ is the wavelength in nm. Equation (74) has been derived by fitting the spectral molecular optical thickness tabular data provided in [14] for the case of the Antarctica Dome C conditions.

The spectral imaginary part of the ice refractive index has been found as proposed in [5] using the dataset presented in [36] at wavelengths below 600 nm and that given in [37] at larger wavelengths. The linear interpolation has been applied to derive the values between the spectral grid points presented in respective tables.

Figure 7 shows that OLCI spectra are significantly influenced by the gaseous absorption (see, e.g., the oxygen absorption band located around 761 nm). The decrease in OLCI spectra towards 400 nm is due to atmospheric light scattering. This decrease can be used to estimate the atmospheric aerosol load over snow for elevated atmospheric aerosol load conditions. As a matter of fact, the aerosol load in Antarctica is in general too low [15] to be measured from a satellite taking into account instrument calibration errors and inaccuracy of bottom-of-atmosphere albedo estimations. The influence of atmospheric scattering can be neglected at wavelengths above 500 nm for the case studied. Measurements at 1020 nm are affected by the light absorption effects in snow cover with larger values of reflectance for smaller and less absorbing ice grains. Measurements in the vicinity of ozone (600 nm) and water vapor (940 nm) absorption bands can be used to estimate the total ozone and the atmospheric precipitable water vapor content. They can be derived from the depth of corresponding absorption bands as seen in the reflected light [5,38]. The atmospheric oxygen column is known in advance for any height of the underlying surface. Therefore, the measurements of reflectance in the oxygen A-band can be used to estimate the height of the underlying surface, and ground atmospheric pressure and develop various cloud screening procedures.

We have found that the difference between theoretical results derived using Equation (2) and observations for all channels is smaller than 5% except for the channel corresponding to the maximal absorption by oxygen in the OLCI spectrum, where one needs to use the instrument spectral response function (instead of box function) in the averaging procedure. It should be pointed out that some differences seen in Figure 7 are due to errors in OLCI measurements. In particular, theoretical results suggest that the spectral reflectance must decrease in the range of 400–450 nm (towards a shorter wavelength). This is mostly due to molecular light scattering effects for the case studied. Such a decrease is not seen in the OLCI spectrum given in Figure 7 by red crosses, which makes it possible to suspect the problem with OLCI measurements in this spectral range. As a matter of fact, such an anomaly in OLCI spectra has already been detected. In particular, the use of OLCI calibration coefficients (gains), which leads to a decrease in OLCI reflectance in the range of 400–450 nm and a change in OLCI spectral reflectance at some other channels, have been ellaborated in [39]. As it follows from Figure 7, the difference between the modelled and measured spectra decreases after the application of gains (except for the last channel). The difference at the last OLCI channel is due to the fact that the value of L in Table 2 has been derived from OLCI measurements without the application of gains. This difference is reduced in case the value of L is retrieved from the measurements at the last OLCI channel after the application of the gains.

The total ozone column and precipitable water vapor content are also important atmospheric parameters, which can be found. It follows that the value of TOC = 250 DU for the case shown in Figure 7. ECMWF data for DOME C provided in OLCI files at this date and time give TOC = 291 DU. Therefore, the difference is around 14% for the TOC for the case studied. The study of errors related to the determination of TOC from OLCI measurements has been performed in [27]. The derived value of PWV for the case shown in Figure 7 is 0.33 mm. This value is close to the average value of PWV at the site for autumn (PWV = 0.4 mm), where averaging has been performed for 4 years (2010–2014, see Table 4 presented in [40]). Further discussion of the PWV content as retrieved using OLCI over land surfaces and related errors including the influence of instrumental effects on the retrieval is given in [38].

As a matter of fact, the theory presented above makes it possible to estimate the product X = mQN accurately. The correct value of N can be derived from X: N = X/Qm using improved estimation of the AMF m and also the parameter Q. Further discussion of this topic is out of the scope of this paper.

The effective grain size can be derived from the value of the effective absorption length using the following approximate equation [41]: d = L/16. Therefore, the effective grain diameter d = 0.14 mm for the case shown in Table 2, which is consistent with grain size measurements in Antarctica as reported in [42].

We also applied the theory described above to a large OLCI dataset in Eastern Antarctica. Figure 8 shows the distribution of the coefficient of variance (normalised root-mean-square difference between model and spaceborne measurements of TOA reflectance outside the oxygen absorption band):

$$CV=\frac{RMSD}{\overline{R}},$$

(75)

where

$$RMSD=\sqrt{\frac{{\displaystyle \sum _{j=1}^J{(R_{appr}({\lambda }_j)-R_{meas}({\lambda }_j))}^2}}{J}},$$

(76)

$$\overline{R}=\frac{{\displaystyle \sum _{j=1}^J{R}_{meas}^{}{\left({\lambda }_j\right)}_{}}}{J}$$

(77)

and J = 18 (oxygen absorption A-band measurements are excluded). As it follows from the comparison of Figure 8 and Figure 9, the largest values of CVs are correlated with the largest values of OLCI reflectance in the oxygen A-band signifying the presence of clouds [43]. The parameters L, $N_{O_3}$,$N_{H_{2}O}^{ef}{}_{}$ for the OLCI scene studied have been estimated using the following equations, which follow from Equations (2), (46), (54) and (65) under the assumption that the atmospheric scattering contribution is neglected [34,37]:

$$L=\frac{{\ln }^2(R_0/R_{21}^{})}{{\alpha }_{21}{\xi }^2},$$

(78)

$$N_{O_3}=\mathrm{K}\frac{\ln (R_{s,7}/R_7)}{m_{O_3}},$$

(79)

$$N_{H_{2}O}^{ef}{}_{}={\ln }^{\epsilon }(R_{s,20}/R_{20}),$$

(80)

where $\mathrm{K}$ = 9349.3 DU, $\epsilon$ = 1/$x_1$, ${\alpha }_{21}=\frac{4\pi {\chi }_{21}}{{\lambda }_{21}}$ is the bulk ice absorption coefficient 1020 nm, ${\chi }_{21}$ is the imaginary part of the ice refractive index at the OLCI channel ${\lambda }_{21}=1020 \mathrm{nm},$ $\xi$ is given by Equation (66), $R_7$, $R_{20},$$R_{21}$ are the OLCI TOA measurements at the wavelengths 620, 940 and 1020 nm, respectively, $R_{s,7}$$R_{s,20}$$R_{s,21}$ are the estimated OLCI TOA reflectances at the OLCI channels 7 (620 nm), 20 (940 nm) and 21 (1020 nm). Air mass factors in Equations (79) and (80) are estimated using the geometrical approximation (see Equation (6)). The value of molecular optical thickness has been modelled using Equation (74). The imaginary part of the ice refractive index has been found as specified above. The parameter $N_{O_2}^{ef}$ has been estimated using the formula similar to Equation (80) except OLCI measurements and snow reflectance at the channel 764.38 nm have been used. Other parameters have been selected as specified in Table 2.

To increase the accuracy of the derived spectral reflectance one can substitute the approximate value of $R_0$ by that estimated from OLCI measurements [34]:

$$R_0=R_{meas}^z(865 \mathrm{nm})R_{meas}^{1-z}(1020\mathrm{nm}),$$

(81)

where $z={(1-b)}^{-1}$, $b=\sqrt{\frac{{\alpha }_{17}}{{\alpha }_{21}}}$, ${\alpha }_{17}$ is the bulk ice absorption coefficient at channel 17 of OLCI (865 nm), $R_{meas}^{}$ is the measured OLCI reflectance at a given OLCI channel.

The histogram of the coefficient of variance corresponding to the data shown in Figure 8 is given in Figure 10. We see that the values of CV are smaller than 10% for most cases. There are regions in Eastern Antarctica, where the deviations are below 5%. The areas with larger deviations could be due to the presence of clouds or structures on the snow surface not captured by the proposed model. The differences can be further decreased if OLCI gains are used as discussed above.

4. Conclusions

In this work, we have proposed an analytical approximation for the clear sky top-of-atmosphere spectral reflectance as measured by satellite instrumentation at a moderate spectral resolution in the visible and near-infrared regions of the electromagnetic spectrum (400–1000 nm) over extended snow surfaces, such as those existing in Antarctica. Both the atmosphere and surface are relatively clean in Antarctica, which makes the solution to the studied problem possible. One can also apply the technique to other similar areas with underlying snow and ice surfaces (e.g., Greenland, Tibetan Plateau, etc.). Although the numerical solution of the radiative transfer equation leads to more accurate results as compared to those presented in this work, it does not mean that this results in better accuracy as compared to the experimental measurements of spectral reflectance. This is related to the fact that exact modelling of light reflectance from natural snow and ice surfaces as detected on a satellite is hardly possible because of the various structures on the snow surfaces (e.g., sastrugi), close-packed media effects, layered nature and complex microstructure of snow fields containing both external and internal impurities. The technique can be extended to larger wavelengths taking into account other absorption bands of atmospheric gases. Additionally, Mie calculations [44] or look-up tables can be used for the spectral aerosol phase function and optical thickness in case the atmospheric aerosol load is substantial and the aerosol parameters are known in advance.

The spectral range 400–1000 nm is of special importance because the TOA reflectance is more sensitive to both snow and atmospheric aerosol properties as compared to longer wavelengths, where the contribution of atmospheric aerosol and molecular light scattering can be neglected and the reflectance of light from snow decreases. The fast radiative transfer FORTRAN code for the calculation of the reflectance over snow (SNOWTRAN) can be freely downloaded from the site https://sites.google.com/site/kokhanovsky2016/codes (accessed on 18 September 2022).