Simulation of snow albedo and solar irradiance profile
with the Two-streAm Radiative TransfEr in Snow
(TARTES) v2.0 model

Ghislain Picard, Quentin Libois

<figure><img src="image_1.png" /></figure>

To cite this version:

Ghislain Picard, Quentin Libois. Simulation of snow albedo and solar irradiance profile with the Two-streAm Radiative TransfEr in Snow (TARTES) v2.0 model. Geoscientific Model Development, 2024, 17 (24), pp.8927-8953. 10.5194/gmd-17-8927-2024. hal-04849062

HAL Id: hal-04849062

https://hal.science/hal-04849062v1

Submitted on 20 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Simulation of snow albedo and solar irradiance profile with the Two-streAm Radiative TransfEr in Snow (TARTES) v2.0 model

Ghislain Picard<fnref n="1" /> and Quentin Libois<fnref n="2" />

**Correspondence:** Ghislain Picard (ghislain.picard@univ-grenoble-alpes.fr)

Received: 19 April 2024 – Discussion started: 25 June 2024

Revised: 23 August 2024 – Accepted: 10 October 2024 – Published: 19 December 2024

**Abstract.** The Two-streAm Radiative TransfEr in Snow (TARTES) model computes the spectral albedo and the profiles of spectral absorption, irradiance, and actinic fluxes for a multi-layer plane-parallel snowpack. Each snow layer is characterized by its specific surface area, density, and impurity content, in addition to shape parameters. In the landscape of snow optical numerical models, TARTES distinguishes itself by taking into account different shapes of the particles through two shape parameters, namely the absorption enhancement parameter *B* and the asymmetry factor *g*. This is of primary importance as recent studies working at the microstructure level have demonstrated that snow does not behave as a collection of equivalent ice spheres, a representation widely used in other models. Instead, *B* and *g* take specific values that do not correspond to any simple geometrical shape, which leads to the concept of the “optical shape of snow”. Apart from this specificity, TARTES combines well-established radiative transfer principles to compute the scattering and absorption coefficients of pure or polluted snow, as well as the δ-Eddington two-stream approximation to solve the multi-layer radiative transfer equation. The model is implemented in Python, but conducting TARTES simulations is also possible without any programming through the SnowTARTES web application, making it very accessible to non-experts and for teaching purposes. Here, after describing the theoretical and technical details of the model, we illustrate its main capabilities and present some comparisons with other common snow radiative transfer models (AART, DISORT-Mie, SNICAR-ADv3) as a validation procedure. Overall the agreement on the spectral albedo, when in compatible conditions (i.e., with spheres), is usually within 0.02

and is better in the visible and near-infrared range compared to longer wavelengths.

## 1 Introduction

Snow, a porous medium made of ice and air, is by far the most reflective material in the solar spectral range on Earth. Any fluctuation in the snow cover extent or changes in the surface snow properties have consequences for the global radiative budget and the climate (Qu and Hall, 2007; Räisänen et al., 2017). Snow albedo (also known as hemispherical reflectance, Schaepman-Strub et al., 2006) is the primary variable controlling the amount of solar energy absorbed in the snowpack (Flanner et al., 2011). Also of importance is the depth at which this absorption occurs. The deeper solar radiation is absorbed, the less likely the corresponding heat is to be transferred back to the surface by thermal conduction, where it can eventually be evacuated to the atmosphere through longwave emission or turbulent mixing. Hence, the warming and potential melt of the snowpack depend on the vertical profile of absorbed sunlight (Dombrovsky et al., 2019). This profile is often approximated by an exponential function decreasing with depth, with the decay length called the *e*-folding depth or penetration depth (Kokhanovsky, 2022). This quantity is also of great importance for photochemical processes (King and Simpson, 2001; Domine et al., 2008).

Snow optical properties are driven by the physical properties of the snow and the impurities it contains, along with the illumination conditions. The snow microstructure (i.e., the arrangement of ice and air at the micrometer scale) controls the absorption of radiation in the near-infrared range

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, IGE, 38000 Grenoble, France</footnote>

<footnote n="2">²CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France</footnote>
(Wiscombe and Warren, 1980) and is key to understanding some of the snow–albedo feedback loops that amplify climate change in snow-covered regions (Hall, 2004; Qu and Hall, 2007; Picard et al., 2012; Box et al., 2022). It is common to represent snow microstructure as a collection of grains with prescribed shape and size (Warren and Wiscombe, 1980; Grenfell and Warren, 1999). However, it is known that snow on the ground is generally not granular, so equivalent concepts have emerged to quantify more general microstructures. For instance, the specific surface area (SSA) (the ratio between the total air-interface surface area and the mass of ice, Domine et al., 2006) advantageously replaces the grain size as it can be rigorously defined and calculated for any porous medium made of two phases (ice and air here) even when distinct grains are not apparent. Regarding the shape, the situation is less advanced, but geometrical metrics related to the chord length distribution can provide useful information (Malinka, 2014; Krol and Löwe, 2016; Dumont et al., 2021; Robledano et al., 2023).

The thickness of the snow cover is another major driving variable in the case of shallow covers and a dark underlying ground, especially in the visible range, where light penetrates deepest (Perovich, 2007). Another important variable is the roughness of the surface, which tends to decrease the albedo (Warren et al., 1998; Leroux and Fily, 1998; Larue et al., 2020). The presence of liquid water also slightly changes the absorption coefficient in a few spectral bands (e.g., 980–1000 nm), which can be detected in spectral albedo measurements (Dumont et al., 2017; Donahue et al., 2022). More importantly it induces fast structural changes of the microstructure (Colbeck, 1982; Brun et al., 1989), usually leading to a rapid decrease in SSA and thus in albedo. Impurities, whether they are of mineral, organic, or biological origin, can also greatly affect the absorption in the visible range and drive both the albedo and the penetration (Chevrollier et al., 2022; Réveillet et al., 2022; Di Mauro et al., 2024). At last, the illumination characteristics (angular and spectral distributions) also play a role because the snow reflectance depends on the wavelength and the incidence angle. As such, neither broadband nor spectral snow albedo is strict surface snow properties, because of this dependency on the illumination characteristics. It means that even without any change in the snow properties, the surface albedo may change with changing environmental conditions (presence of clouds, sun elevation, etc.).

To account for these multiple factors, numerous snow optical models have been developed (e.g., He and Flanner, 2020). Even though most rely on the radiative transfer (RT) principles, they greatly differ in the method used to solve the RT equation (RTE, e.g., Chandrasekhar, 1960), in the representation of the medium (snow microstructure, the 1D or 3D geometry of the snowpack, surface roughness, heterogeneity of the snowpack, impurities, etc.), in the fundamental constants used (e.g., real and imaginary parts of the ice refractive index), and in the output optical quantities (e.g., albedo,

absorption profile, actinic flux, transmittance). Each model has its niche of applications, from efficient but approximate code suitable to large-scale climate models to very precise solvers to investigate detailed optical behaviors (e.g., bidirectional reflectance distribution function, BRDF). To cite a few, pioneering works used phenomenological (Dunkle and Bevans, 1956; Bohren, 1987) or more rigorous (Wiscombe and Warren, 1980; Warren and Wiscombe, 1980) two-stream approximations to solve the RTE and Mie theory (Mie, 1908) or geometrical optics (Bohren and Barkstrom, 1974) to represent spherical ice particles suspended in the air. DISORT (Stamnes et al., 1988a, b) is a general robust and popular solver frequently applied to snow, usually in combination with the Mie theory, hence assuming spherical grains (Glendinning and Morris, 1999; Green et al., 2002; Gallet et al., 2011; Carmagnola et al., 2013; Dang et al., 2019). Here, we refer to this combination as “DISORT-Mie”. An advantage of DISORT is its ability to account for the radiation propagation in many directions (multi-stream). To account for other particle shapes and calculate the BRDF, an efficient code was proposed by Mishchenko et al. (1999) (available as a Python package here: https://github.com/ghislainp/mishchenko_brf, last access: 24 March 2024), but it cannot handle layering. TUV-snow is a DISORT-based coupled snow–atmosphere model specifically designed for UV radiation and photochemistry applications (Lee-Taylor and Madronich, 2002; France et al., 2011). PBSAM (Aoki et al., 2011) and SNICAR (Flanner and Zender, 2005) are fast two-stream solvers suitable for surface albedo calculation in climate simulations (Onuma et al., 2020; Usha et al., 2020). SNICAR is currently one of the most actively developed codes with a large panel of state-of-the-art parameterizations, for instance to account for snow algae (Cook et al., 2017) or to account for ice layers with the recent replacement of the two-stream solver by the adding–doubling solver (Flanner et al., 2021). Importantly, all these models are unidimensional, meaning that they describe the snowpack as a stack of homogeneous, horizontally infinite, and flat layers. This configuration is known as plane parallel. Conversely, three-dimensional RT models are necessary to account for snowpack with 3D structures, embedded objects, or light sources. Some examples are models for rough surfaces (Warren et al., 1998; Larue et al., 2020; Robledano et al., 2022), for explicit photon trajectory calculations in the snow microstructure (Kaempfer et al., 2005; Picard et al., 2009; Xiong and Shi, 2014; Letcher et al., 2022; Robledano et al., 2023), and for interactions with embedded objects or instruments (Gallet et al., 2009; Picard et al., 2016).

In this rich landscape, the Two-streAm Radiative TransfEr in Snow (TARTES) spectral model differs by two main aspects. First it relies on a simple yet state-of-the-art representation of the snow microstructure, allowing us to represent it with four parameters only: the SSA, the density, the asymmetry factor *g* that quantifies the forward scattering of snow, and the absorption enhancement parameter *B* that quantifies
the lengthening of photon paths inside the ice phase due to multiple internal reflections. As a consequence, TARTES is not restricted to spherical particles, and any pair of g and B values can be used, possibly not corresponding to any particular idealized geometrical shape. It also offers the possibility of conforming with the asymptotic approximation radiative transfer (AART, Kokhanovsky and Zege, 2004) when the snowpack is semi-infinite, while still being able to simulate a multi-layered snowpack. The second main difference from other models is the use of the Python language. This facilitates rapid tests (e.g., in notebooks) and implementation of new features (e.g., impurities, Tuzet et al., 2019). While TARTES has been used over a decade (e.g., Libois et al., 2013; Shao et al., 2018; van Dalum et al., 2019; Tuzet et al., 2020; Manninen et al., 2021; Veillon et al., 2021), this paper aims at a first comprehensive and formal description of the model. Some minor but significant adjustments using the latest results on microstructure (Robledano et al., 2023) are also included, as well as a brief presentation of the ecosystem of tools relevant to snow optical computations built around TARTES.

Section 2 provides the detailed derivation of the model; Sect. 3 presents the TARTES software and the associated ecosystem; and Sect. 4 presents the results including self-consistency checks, a comparison with other models, and a presentation of its specific capabilities. Section 5 discusses how TARTES fits with the concept of the “optical shape of snow” and presents the model limitations. Section 6 concludes this study.

## 2 The physics behind TARTES

TARTES relies on the $\delta$-Eddington approximation to solve the plane-parallel RTE and compute the spectral upwelling and downwelling fluxes within a multi-layered snowpack. To this end the single-scattering properties of each layer are computed from the SSA, density, snow grain shape, and amount of light-absorbing impurities (Fig. 1). This section provides all the theoretical details on which TARTES is built, as well as new formulations added in version 2.0.

### 2.1 The plane-parallel radiative transfer equation

The steady-state unpolarized RTE describes the intensity (or radiance) field $I$ propagating in an absorbing and scattering slab at depth $z$ in a direction $(\theta, \phi)$, where $\phi$ is the azimuth angle, $\theta$ is the zenith angle defined as the angle between the inward normal to the surface and the direction of light propagation, and $z$ is defined positive from the surface downward. Such a medium is characterized by its extinction coefficient $\sigma_e$ (m$^{-1}$), scattering coefficient $\sigma_s$ (m$^{-1}$), and scattering phase function $p(\theta, \phi, \theta', \phi')$ (unitless). The absorption coefficient is obtained by energy conservation as $\sigma_a = \sigma_e - \sigma_s$. The phase function describes the probability

for light to be scattered into the direction $(\theta, \phi)$ when coming from the direction $(\theta', \phi')$. Here we consider a horizontal multi-layered snowpack. Each layer of the snowpack is assumed to be isotropic and to have homogeneous optical properties. The snowpack is illuminated at the surface by solar radiation that can be a combination of diffuse (i.e., isotropic) and direct light incident at zenith angle $\theta_0$ and azimuth angle $\phi_0$. There are no internal light sources in the snowpack, and no thermal emission is considered since TARTES focuses on the shortwave range (200–4000 nm). Along the direction $s$ defined by $(\theta, \phi)$, $I$ decreases due to extinction (absorption and scattering) and increases due to scattering from all other directions $(\theta', \phi')$, so the RTE reads

$$ \frac{dI(z, \theta, \phi)}{ds} = -\sigma_e I(z, \theta, \phi) + \frac{\sigma_s}{4\pi} \int_{0}^{\pi} \int_{0}^{2\pi} p(\theta, \phi, \theta', \phi') \times I(z, \theta', \phi') d\phi' \sin\theta' d\theta', \quad (1) $$

where the phase function is normalized so that $\frac{1}{4\pi} \int_{0}^{\pi} \int_{0}^{2\pi} p(\theta, \phi, \theta', \phi') d\phi' \sin\theta' d\theta' = 1$. Defining $\mu = \cos\theta$ and $\mu' = \cos\theta'$, noting that $dz = \mu ds$, and further defining the optical thickness such that $d\tau = \sigma_e dz$, Eq. (1) becomes

$$ \mu \frac{dI(\tau, \mu, \phi)}{d\tau} = -I(\tau, \mu, \phi) + \frac{\omega}{4\pi} \int_{-1}^{1} \int_{0}^{2\pi} p(\mu, \phi, \mu', \phi') \times I(\tau, \mu', \phi') d\phi' d\mu', \quad (2) $$

where the single-scattering albedo $\omega = \sigma_s/\sigma_e$. Since TARTES focuses on radiative fluxes through horizontal surfaces, Eq. (2) can be azimuthally integrated, which reads

$$ \mu \frac{dI(\tau, \mu)}{d\tau} = -I(\tau, \mu) + \frac{\omega}{2} \int_{-1}^{1} p(\mu, \mu') I(\tau, \mu') d\mu', \quad (3) $$

where we have defined the azimuthally averaged intensity $I(\tau, \mu) = \frac{1}{2\pi} \int_{0}^{2\pi} I(\tau, \mu, \phi) d\phi$ and the azimuth-independent

phase function $p(\mu, \mu') = \frac{1}{2\pi} \int_{0}^{2\pi} p(\mu, \phi, \mu', \phi') d\phi$.

When the snowpack is at least partly illuminated by a beam source (e.g., direct solar radiation), it is useful to write $I(\tau, \mu) = I_{\text{dir}}(\tau, \mu) + I_{\text{diff}}(\tau, \mu)$, where the direct intensity $I_{\text{dir}}$ corresponds to light that has not been scattered and $I_{\text{diff}}$ is the diffuse intensity. At the surface the direct intensity is $F_0 \delta(\mu - \mu_0) \delta(\phi - \phi_0)$, where $\mu_0 = \cos\theta_0$, and at depth

$$ I_{\text{dir}}(\tau, \mu) = \frac{F_0}{2\pi} \delta(\mu - \mu_0) e^{-\tau/\mu_0}, \quad (4) $$

where $F_0$ is the intensity of the solar beam at the surface and $\delta$ the Dirac function.
<figure><img src="image_4.png" /><figcaption>Figure 1. Main inputs, options, and flow of computations in TARTES v2.0.</figcaption></figure>

Reporting Eq. (4) in Eq. (3), we obtain the RTE for the diffuse intensity:

$$ \mu \frac{dI_{\text{diff}}(\tau, \mu)}{d\tau} = -I_{\text{diff}}(\tau, \mu) + \frac{\omega}{2} \int_{-1}^{1} p(\mu, \mu') I_{\text{diff}}(\tau, \mu') d\mu' + \frac{\omega}{4\pi} p(\mu, \mu_0) F_0 e^{-\tau/\mu_0}. \quad (5) $$

From now on $I_{\text{diff}}$ will simply be referred to as $I$.

## 2.2 The δ-Eddington approximation of the phase function

Within each layer, snow is assumed to be isotropic, so the phase function depends only on the scattering angle $\Theta$ between the incident and scattered light, and we can write $p(\mu, \phi, \mu', \phi') = p(\cos \Theta)$. This angle is such that

$$ \cos \Theta = \mu\mu' + \sqrt{(1 - \mu^2)(1 - \mu'^2)} \cos(\phi - \phi'). \quad (6) $$

$p(\cos \Theta)$ can be expanded in Legendre polynomials $P_l$:

$$ p(\cos \Theta) = \sum_{l=0}^{\infty} \omega_l P_l(\cos \Theta), \\ \text{where } \omega_l = \frac{2l + 1}{2} \int_{-1}^{1} p(\cos \Theta) P_l(\cos \Theta) d\cos \Theta. \quad (7) $$

By virtue of normalization $\omega_0 = 1$, and the mean cosine of the scattering angle of the phase function, called the asymmetry factor $g$, is such that $\omega_1 = 3g$. Using the addition theorem of spherical harmonics (Chandrasekhar, 1960), it can finally be shown that

$$ p(\mu, \mu') = \sum_{l=0}^{\infty} \omega_l P_l(\mu) P_l(\mu'). \quad (8) $$

Hence, the two-term truncation of the phase function reads

$$ p(\mu, \mu') = 1 + 3g\mu\mu'. \quad (9) $$

To handle the strong forward scattering of snow particles, TARTES relies on the δ-Eddington approximation, which consists in writing the phase function as the sum of a strictly forward scattering component (a Dirac) and a two-term phase function (Eq. 9). Joseph et al. (1976) proposed weighting both contributions so that the asymmetry factor is conserved and the second moment of the phase function equals $g^2$ (i.e., the second moment of the Henyey-Greenstein phase function Henyey and Greenstein, 1941, with asymmetry factor $g$). This reads

$$ p(\mu, \mu') = 2g^2\delta(\mu - \mu') + (1 - g^2)(1 + 3g^*\mu\mu'), \quad (10) $$

with $g^* = \frac{g}{1+g}$.

Combining Eqs. (5) and (10) we obtain

$$ \mu \frac{dI(\tau^*, \mu)}{d\tau} = -I(\tau^*, \mu) + \frac{\omega^*}{2} \int_{-1}^{1} (1 + 3g^* \mu \mu') I(\tau^*, \mu') d\mu' + \frac{\omega^*}{4\pi} (1 + 3g^* \mu \mu_0) F_0 e^{-\tau^*/\mu_0}, \quad (11) $$

where the following variable changes have been made:

$$ \tau^* = \tau(1 - \omega g^2), \quad (12) $$

$$ \omega^* = \frac{(1 - g^2)\omega}{(1 - \omega g^2)}. \quad (13) $$

Hence, the δ-Eddington approximation of the phase function consists in solving Eq. (5), with $\tau$, $\omega$, and $g$ replaced by $\tau^*$, $\omega^*$, and $g^*$. $g^*$ is less than $g$, so the scaled phase function is less forward peaking than the original phase function, which reduces the errors in the following two-stream resolution of the RTE. Note that in this approximation the solution for direct radiation is scaled accordingly, meaning that direct radiation can propagate deeper in the snowpack, because light scattered in the forward direction is treated as unscattered light.
## 2.3 Equations for fluxes and Eddington approximation

In TARTES, we are interested in the vertical downward and upward fluxes in the snowpack, $F^-$ and $F^+$ respectively. These quantities are defined as

$$F^{-}(\tau^{*}) = 2\pi \int_{0}^{1} I(\tau^{*}, \mu) d\mu, \qquad (14)$$

$$F^{+}(\tau^{*}) = 2\pi \int_{0}^{1} I(\tau^{*}, -\mu) d\mu. \qquad (15)$$

Integrating Eq. (11) over both positive and negative values of $\mu$ results in two differential equations:

$$\begin{aligned} \frac{dF^{-}(\tau^{*})}{d\tau^{*}} = & -2\pi \int_{0}^{1} I(\tau^{*}, \mu) d\mu \\ & + \pi\omega^{*} \int_{0}^{1} \int_{-1}^{1} (1 + 3g^{*}\mu\mu') I(\tau^{*}, \mu') d\mu' d\mu \\ & + \omega^{*}\gamma_{4} F_{0} e^{-\tau^{*}/\mu_{0}}, \end{aligned} \qquad (16)$$

$$\begin{aligned} \frac{dF^{+}(\tau^{*})}{d\tau^{*}} = & 2\pi \int_{0}^{1} I(\tau^{*}, -\mu) d\mu \\ & - \pi\omega^{*} \int_{0}^{1} \int_{-1}^{1} (1 - 3g^{*}\mu\mu') I(\tau^{*}, \mu') d\mu' d\mu \\ & - \omega^{*}\gamma_{3} F_{0} e^{-\tau^{*}/\mu_{0}}, \end{aligned} \qquad (17)$$

with

$$\gamma_{4} = \frac{1}{4}(2 + 3g^{*}\mu_{0}) \quad \text{and} \quad \gamma_{3} = \frac{1}{4}(2 - 3g^{*}\mu_{0}). \qquad (18)$$

Next the Eddington approximation is used, which consists in expanding the intensity $I(\tau^*, \mu)$ as

$$I(\tau^*, \mu) = I_0(\tau^*) + \mu I_1(\tau^*), \qquad (19)$$

so that

$$F^{-}(\tau^{*}) = 2\pi \left[ \frac{I_{0}(\tau^{*})}{2} + \frac{I_{1}(\tau^{*})}{3} \right], \qquad (20)$$

$$F^{+}(\tau^{*}) = 2\pi \left[ \frac{I_{0}(\tau^{*})}{2} - \frac{I_{1}(\tau^{*})}{3} \right]. \qquad (21)$$

This reads

$$2\pi I(\tau^*, \pm\mu) = \frac{1}{2}[(2 \pm 3\mu)F^{-}(\tau^*) + (2 \pm 3\mu)F^{+}(\tau^*)], \qquad (22)$$

and therefore

$$2\pi \int_{0}^{1} I(\tau^{*}, \pm\mu) d\mu = \frac{1}{4}[(4 \pm 3)F^{-}(\tau^{*}) + (4 \mp 3)F^{+}(\tau^{*})]. \qquad (23)$$

Eventually,

$$\begin{aligned} & \pi\omega^{*} \int_{0}^{1} \int_{-1}^{1} (1 \pm 3g^{*}\mu\mu') I(\tau^{*}, \mu') d\mu' d\mu \\ = & \frac{\omega^{*}}{4} [(4 \pm 3g^{*})F^{-}(\tau^{*}) + (4 \mp 3g^{*})F^{+}(\tau^{*})]. \end{aligned} \qquad (24)$$

Substituting Eqs. (23) and (24) into Eqs. (16) and (17) we obtain

$$\begin{aligned} \frac{dF^{-}(\tau^{*})}{d\tau^{*}} = & -\frac{1}{4}[7F^{-}(\tau^{*}) + F^{+}(\tau^{*})] \\ & + \frac{\omega^{*}}{4}[(4 + 3g^{*})F^{-}(\tau^{*}) + (4 - 3g^{*})F^{+}(\tau^{*})] \\ & + \omega^{*}\gamma_{4} F_{0} e^{-\tau^{*}/\mu_{0}}, \end{aligned} \qquad (25)$$

$$\begin{aligned} \frac{dF^{+}(\tau^{*})}{d\tau^{*}} = & \frac{1}{4}[F^{-}(\tau^{*}) + 7F^{+}(\tau^{*})] \\ & - \frac{\omega^{*}}{4}[(4 - 3g^{*})F^{-}(\tau^{*}) + (4 + 3g^{*})F^{+}(\tau^{*})] \\ & - \omega^{*}\gamma_{3} F_{0} e^{-\tau^{*}/\mu_{0}}, \end{aligned} \qquad (26)$$

which can be factorized as

$$\frac{dF^{-}(\tau^{*})}{d\tau^{*}} = \gamma_{2}F^{+}(\tau^{*}) - \gamma_{1}F^{-}(\tau^{*}) + \omega^{*}\gamma_{4}F_{0}e^{-\tau^{*}/\mu_{0}}, \qquad (27)$$

$$\frac{dF^{+}(\tau^{*})}{d\tau^{*}} = \gamma_{1}F^{+}(\tau^{*}) - \gamma_{2}F^{-}(\tau^{*}) - \omega^{*}\gamma_{3}F_{0}e^{-\tau^{*}/\mu_{0}}, \qquad (28)$$

where

$$\begin{aligned} \gamma_{1} = & \frac{1}{4}[7 - \omega^{*}(4 + 3g^{*})] \\ \text{and} \quad \gamma_{2} = & -\frac{1}{4}[1 - \omega^{*}(4 - 3g^{*})]. \end{aligned} \qquad (29)$$

This corresponds to two coupled first-order differential equations, with matrix A such that

$$A = \begin{pmatrix} -\gamma_1 & \gamma_2 \\ -\gamma_2 & \gamma_1 \end{pmatrix}, \qquad (30)$$

which has two eigenvalues $k_e$ and $-k_e$, with $k_e = \sqrt{\gamma_1^2 - \gamma_2^2}$, and corresponding eigenvectors

$$v_1 = \begin{pmatrix} 1 \\ 1/\Gamma \end{pmatrix}, \quad \text{and} \quad v_2 = \begin{pmatrix} 1 \\ \Gamma \end{pmatrix}, \qquad (31)$$
where $\Gamma = \frac{\gamma_1 - k_e}{\gamma_2}$. A particular solution of this system is sought as

$$F_p^-(\tau^*) = G^- e^{-\tau^*/\mu_0}, \quad (32)$$

$$F_p^+(\tau^*) = G^+ e^{-\tau^*/\mu_0}. \quad (33)$$

Inserting these expressions into Eqs. (27) and (28) results in two equations with two unknowns, which give $G^-$ and $G^+$:

$$G^- = \frac{\mu_0^2 \omega^* F_0}{(k_e \mu_0)^2 - 1} [(\gamma_1 + 1/\mu_0)\gamma_4 + \gamma_2 \gamma_3], \quad (34)$$

$$G^+ = \frac{\mu_0^2 \omega^* F_0}{(k_e \mu_0)^2 - 1} [(\gamma_1 - 1/\mu_0)\gamma_3 + \gamma_2 \gamma_4]. \quad (35)$$

This overall gives the following solutions for the total (i.e., sum of direct and diffuse) downward and upward fluxes:

$$F_{\text{tot}}^-(\tau^*) = Ae^{-k_e \tau^*} + Be^{k_e \tau^*} + (G^- + \mu_0 F_0) e^{-\tau^*/\mu_0}, \quad (36)$$

$$F_{\text{tot}}^+(\tau^*) = \Gamma A e^{-k_e \tau^*} + \frac{B}{\Gamma} e^{k_e \tau^*} + G^+ e^{-\tau^*/\mu_0}, \quad (37)$$

where A and B are unknowns that will be determined according to the boundary conditions. Note that these above solutions are consistent with those used by Toon et al. (1989). In addition to the fluxes, the actinic flux can be derived:

$$F_{\text{act}}(\tau^*) = 2\pi \int_{-1}^{1} I(\tau^*, \mu) d\mu. \quad (38)$$

Given Eq. (19), and considering the contribution of the direct radiation to the actinic flux, this overall reads

$$F_{\text{act}}(\tau^*) = 2(F^-(\tau^*) + F^+(\tau^*)) + F_0 e^{-\tau^*/\mu_0}. \quad (39)$$

## 2.4 Alternative formulation to match the AART theory

The form of Eqs. (36) and (37) is common to all two-stream methods (e.g., Meador and Weaver, 1980). Considering a semi-infinite snowpack, it is clear that $\Gamma$ corresponds to the asymptotic diffuse albedo and $k_e$ to the asymptotic flux extinction coefficient. To allow a perfect match of the two-stream solution with the AART theory in the case of a single layer, we propose testing a new variant denoted by TARTES$\infty$AART in the following that uses alternative expressions for $\Gamma$, $k_e$, $G^-$, and $G^+$. More specifically, following Kokhanovsky and Zege (2004) we let

$$\Gamma = \exp \left( -4 \sqrt{\frac{1-\omega}{3(1-g)}} \right), \quad (40)$$

$$k_e = \sqrt{3(1-\omega)(1-g)}. \quad (41)$$

The parameters $G^-$ and $G^+$ are chosen so that the direct albedo of the AART theory is obtained in the case of a semi-infinite snowpack, which is given by

$$\alpha_{\text{dir}}(\mu_0) = \exp \left( -\frac{12}{7}(1+2\mu_0) \sqrt{\frac{1-\omega}{3(1-g)}} \right). \quad (42)$$

It implies that

$$\Gamma A + G^+ = \alpha_{\text{dir}} \mu_0 F_0. \quad (43)$$

We also have $A + G^- = 0$ because the incident diffuse radiation is zero, but we need another constraint on $G^-$ and $G^+$. We set their sum equal to that of the $\delta$-Eddington approximation, which is

$$G^- + G^+ = \frac{3}{2} G_0 (1 + g(1 - \omega)), \quad (44)$$

where

$$G_0 = \frac{\mu_0^2 \omega F_0}{(k_e \mu_0)^2 - 1}. \quad (45)$$

This finally reads

$$G^- = \frac{\frac{3}{2} G_0 (1 + g(1 - \omega)) - \alpha_{\text{dir}} \mu_0 F_0}{\Gamma + 1}, \quad (46)$$

$$G^+ = \frac{3}{2} G_0 (1 + g(1 - \omega)) - G^- . \quad (47)$$

When these new formulas are used, $\tau^*$ in the previous Eqs. (11)–(39) must be changed to $\tau$, without $\delta$ scaling. Note that we also tested the $\delta$ scaling with these AART formulas and found similar performance to the present TARTES$\infty$AART formulas without scaling. These results are not reported.

## 2.5 Extension to a multi-layered snowpack

The equations derived so far all considered a unique homogeneous layer. For a multi-layered snowpack, the fluxes within each layer of the snowpack have the general form given by Eqs. (36) and (37), but to determine the actual fluxes, the constants A and B should be determined for each layer, which amounts to $2N$ unknowns ($A_i, B_i, i \in \{1, n\}$) for a snowpack with $N$ layers. These unknowns are deduced from the continuity of $F_{\text{tot}}^\pm(\tau^*)$ at the layer interfaces $(2(N-1)$ equations) and the top and bottom boundary conditions (2 equations). Continuity of the diffuse fluxes at $\tau_i^*$ between layers $i$ and $i+1$ reads

$$A_i e^{-k_{e,i} \tau_i^*} + B_i e^{k_{e,i} \tau_i^*} + G_i^- e^{-\tau_i^*/\mu_0} = A_{i+1} e^{-k_{e,i+1} \tau_i^*} + B_{i+1} e^{k_{e,i+1} \tau_i^*} + G_{i+1}^- e^{-\tau_i^*/\mu_0}, \quad (48)$$

$$\Gamma_i A_i e^{-k_{e,i} \tau_i^*} + \frac{B_i}{\Gamma_i} e^{k_{e,i} \tau_i^*} + G_i^+ e^{-\tau_i^*/\mu_0} = \Gamma_{i+1} A_{i+1} e^{-k_{e,i+1} \tau_i^*} + \frac{B_{i+1}}{\Gamma_{i+1}} e^{k_{e,i+1} \tau_i^*} + G_{i+1}^- e^{-\tau_i^*/\mu_0}. \quad (49)$$

From now on we use the notation $k_i = k_{e,i}$ and define $A'_i = A_i e^{-k_i \tau_{i-1}^*}$ and $B'_i = B_i e^{k_i \tau_{i-1}^*}$. Note also that $\tau_0 = 0$. The boundary conditions at the top of the snowpack, where the
diffuse flux is $F_0^{\text{diff}}$, and at the bottom, where the underlying surface is assumed to be Lambertian and characterized by its albedo $\alpha_b$, read

$$A_1 + B_1 + G_1^- = F_0^{\text{diff}}, \quad (50)$$

$$ \begin{aligned} & \Gamma_N A_N e^{-k_N \tau_N^*} + \frac{B_N}{\Gamma_N} e^{k_N \tau_N^*} + G_N^+ e^{-\tau_N^* / \mu_0} \\ & = \alpha_b \left( A_N e^{-k_N \tau_N^*} + B_N e^{k_N \tau_N^*} + (G_N^- + \mu_0 F_0) e^{-\tau_N^* / \mu_0} \right). \end{aligned} \quad (51) $$

The linear system formed by these 2N independent equations can be written as

$$MX = V, \quad (52)$$

with

$$X = {}^t(A'_1, B'_1, \dots, A'_i, B'_i, \dots, A'_N, B'_N). \quad (53)$$

The matrix M reads

$$ \begin{pmatrix} \frac{1}{e_1^-} & \frac{1}{\tau_1^+} & 0 & 0 & 0 & \dots & 0 & 0 \\ \frac{1}{\Gamma_1 e_1^-} & \frac{1}{\Gamma_1 \tau_1^+} & -1 & -1 & 0 & \dots & 0 & 0 \\ 0 & 0 & -\frac{e_2^-}{\Gamma_2} & \frac{e_2^+}{-1/\Gamma_2} & 0 & \dots & 0 & 0 \\ 0 & 0 & \frac{\Gamma_2 e_2^-}{\Gamma_2 \tau_2^+} & \frac{\Gamma_2 e_2^+}{\tau_2^-} & \dots & \dots & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \ddots \\ 0 & 0 & 0 & 0 & \dots & \dots & \frac{-\Gamma_N}{(\Gamma_N - \alpha_b) e_N^-} & \frac{-\Gamma_N}{-1/(\Gamma_N - \alpha_b) e_N^+} \\ 0 & 0 & 0 & 0 & \dots & \dots & \frac{-\Gamma_N}{(\Gamma_N - \alpha_b) e_N^-} & \frac{-\Gamma_N}{(1/\Gamma_N - \alpha_b) e_N^+} \end{pmatrix}. \quad (54) $$

and

$$ V = {}^t(F_0^{\text{diff}} - G_1^-, \dots, dG_i^- e^{-\tau_i^* / \mu_0}, dG_i^+ e^{-\tau_i^* / \mu_0}, \dots, [\alpha_b(G_N^- + \mu_0 F_0) - G_N^+] e^{-\tau_N^* / \mu_0}), \quad (55) $$

where we used the notation $e_i^{\pm} = e^{\pm k_i d\tau_i^*}$, $d\tau_i^* = \tau_i^* - \tau_{i-1}^*$ is the optical depth of layer i, and $dG_i^{\pm} = G_i^{\pm} - G_i^{\pm}$. The matrix M can be tridiagonalized by consecutively performing the following replacement operations on the lines $M_j$ of M for an even $2 \le j < N$:

$$1. \quad M_j - \Gamma_{i/2+1} L_{j+1} \rightarrow M_j,$$

$$2. \quad (1 - \Gamma_{i/2} \Gamma_{i/2+1}) M_{j+1} - \Gamma_{i/2} M_j \rightarrow M_{j+1}.$$

The new matrix M reads

$$ \begin{pmatrix} \frac{1}{(1-\Gamma_1 \Gamma_2) e_1^-} & \frac{1}{(1-\Gamma_2 \Gamma_1) e_1^+} & \frac{0}{(\Gamma_1^2 - \Gamma_2^2)} & \frac{0}{(\Gamma_1 - \Gamma_2) e_1^-} & \frac{0}{(\Gamma_1 - \Gamma_2) e_1^+} & \dots & 0 & 0 \\ 0 & 0 & \frac{0}{(1-\Gamma_1 \Gamma_2) e_2^-} & \frac{0}{(1-\Gamma_1^2 - \Gamma_2^2) e_2^-} & \frac{(\alpha_1 - 1/\alpha_2)}{\Gamma_1 \Gamma_2} & \dots & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots & \dots & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots & \dots & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \ddots \\ 0 & 0 & 0 & 0 & \dots & \dots & \frac{(\Gamma_1^2 - 1)}{(\Gamma_1^2 - \alpha_b) e_N^-} & \frac{(\Gamma_1^2 - \alpha_b) e_N^-}{(1/\Gamma_1 - \alpha_b) e_N^+} \end{pmatrix} $$

Accordingly, the new vector V reads

$$ \begin{aligned} V = {}^t(F_0^{\text{diff}} - G_1^-, \dots, (dG_i^- - \Gamma_{i+1} dG_i^+) e^{-\tau_i^* / \mu_0}, \\ (dG_i^+ - \Gamma_i dG_i^-) e^{-\tau_i^* / \mu_0}, \dots, [\alpha_b(G_N^- + \mu_0 F_0) - G_N^+] e^{-\tau_N^* / \mu_0}). \end{aligned} \quad (56) $$

The 2N unknowns are efficiently retrieved by inversion of the tridiagonal system. Then the fluxes at each interface are calculated as follows:

$$F_{\text{tot}}(\tau_i^*) = A'_i e^{-k_i^* d\tau_i^*} + B'_i e^{k_i^* d\tau_i^*} + (G_i^- + \mu_0 F_0) e^{-\tau_i^* / \mu_0}, \quad (57)$$

$$F_{\text{tot}}^{+}(\tau_i^*) = \Gamma_i A'_i e^{k_i^* d\tau_i^*} + \frac{B'_i}{\Gamma} e^{k_i^* d\tau_i^*} + G_i^+ e^{-\tau_i^* / \mu_0}. \quad (58)$$

Note that in practice the fluxes can be computed at any requested depth. To this end, the layer i corresponding to this depth is first identified. The ratio of the distance between the above interface and the requested depth, as well as the thickness of the layer, is used to scale $d\tau_i^*$ in the above solutions.

## 2.6 Computed quantities

Beyond the fluxes (and actinic fluxes) that are the native variables returned by the above equations, it is possible to compute the energy absorbed by layer i as the sum of the absorbed upwelling $E_u$ and downwelling $E_d$ energy using

$$E_i = \frac{F_{\text{tot}}^{+}(\tau_i^*) - F_{\text{tot}}^{+}(\tau_{i-1}^*)}{E_u} - \frac{(F_{\text{tot}}(\tau_i^*) - F_{\text{tot}}(\tau_{i-1}^*))}{E_d}, \quad (59)$$

$$ \begin{aligned} E_u = \Gamma_i A'_i (e^{-k_i^* d\tau_i^*} - 1) &+ \frac{B'_i}{\Gamma_i} (e^{k_i^* d\tau_i^*} - 1) \\ &+ G_i^+ (e^{-\tau_i^* / \mu_0} - e^{-\tau_{i-1}^* / \mu_0}), \end{aligned} \quad (60) $$

$$ \begin{aligned} E_d = A'_i (e^{-k_i^* d\tau_i^*} - 1) &+ B'_i (e^{k_i^* d\tau_i^*} - 1) \\ &+ (G_i^- + \mu_0 F_0) (e^{-\tau_i^* / \mu_0} - e^{-\tau_{i-1}^* / \mu_0}), \end{aligned} \quad (61) $$

while the energy absorbed by the ground is given by

$$ E_{\text{bottom}} = (1 - \alpha_b)(A'_N e_N^- + B'_N e_N^+ + (G_N^- + \mu_0 F_0) e^{-\tau_N^* / \mu_0}). \quad (62) $$

The albedo of the snowpack is also calculated as the ratio of the upward to the downward flux at the surface:

$$ \alpha = \frac{1}{\mu_0 F_0 + F_0^{\text{diff}}} \left( \Gamma_1 A_1 + \frac{B_1}{\Gamma_1} + G_1^+ \right). \quad (63) $$

So far we have not specified anything about the spectral dimension of incident light. Implicitly all above derivations are valid for monochromatic radiation, so TARTES is in essence a monochromatic model. Since the single-scattering properties of the snowpack are wavelength-dependent, the matrix M and the vector V are computed at each relevant wavelength. Broadband quantities are thus obtained by summing the contribution of all wavelengths. For instance, the broadband albedo $\bar{\alpha}$ is obtained through spectral integration:

$$ \bar{\alpha} = \frac{\sum_{i=1}^{N} \alpha(\lambda_i)(F_0^{\text{diff}}(\lambda_i) + \mu_0 F_0(\lambda_i))}{\sum_{i=1}^{N} F_0^{\text{diff}}(\lambda_i) + \mu_0 F_0(\lambda_i)}. \quad (64) $$
## 2.7 Treatment of diffuse incident radiation

As seen previously, the incident radiation in TARTES can be direct or diffuse. Wiscombe (1977) has shown that in the case of diffuse radiation the performance of the δ-Eddington approximation was limited, sometimes leading to negative values of albedo. This is why Warren and Wiscombe (1980) and Dang et al. (2019) computed the diffuse albedo as an angular average of direct albedos. In TARTES, the most accurate strategy to handle diffuse radiation is to compute the integrated sum of direct radiation coming from all directions, following an angular distribution such that $p(\theta_0) = \cos\theta_0$. Hence, it requires integrating the solutions for direct incident light at various angles. As only the vector V depends on incident light characteristics, to compute the optical properties of a snowpack at various angles of incidence, M has to be calculated only once, which is computationally relatively efficient.

An alternative strategy proposed in TARTES is to consider that diffuse radiation can be approximated by direct radiation at an effective zenith angle $\theta_{diff}$ such that (see Eq. 42)

$$ \frac{3}{7}(1 + 2\cos\theta_{\text{diff}}) = 1, \qquad (65) $$

which corresponds approximately to an angle of 48.2°. This alternative is the default option in TARTES (hereinafter referred to as 48.2°). Note that in the initial version of TARTES (Libois et al., 2013) the diffuse albedo was calculated using a direct component at an angle of 53° based on equivalence tests using DISORT-Mie. This angle was changed to the approximate value of 48° on 22 June 2022 and is now obtained with the exact calculation (Eq. 65).

Note that despite the problem identified by Wiscombe (1977), the pure diffuse boundary condition of the two-stream method is implemented in TARTES (hereinafter denoted by “2S”) and should be selected for testing only as done in Sect. 4.1.1. To avoid negative albedo, we set $\alpha = 0$ when a negative value is obtained, which occurs when the ice absorption is very large (1400–1600 nm) (Sect. 4.1.1). The other quantities calculated by TARTES are not corrected, and we discourage using the 2S option in general.

## 2.8 Treatment of optically deep layers and snowpacks

When a layer is too thick, the terms $e_i^\pm$ become either extremely large or small, and in both cases they cannot be handled numerically. To avoid this, when a layer is too thick (practically when $k_i d\tau_i^* > 200$), its optical depth is modified so that $k_i d\tau_i^* = 200$. In addition, when a snowpack is very deep, energy does not penetrate through the whole snowpack; it is essentially absorbed in the topmost layers. To save computation time, the snowpack used for the calculations is reduced to the top n layers, where n is the smallest integer such

that

$$ \sum_{i=1}^{n} k_{i}^{*} d\tau_{i} > 30. \qquad (66) $$

At the same time, the optical thickness of the last layer is set to 30/$k_n$ and the underlying albedo is set to 1 to ensure the underlying surface does not absorb energy.

## 2.9 Single-scattering properties of snow

The previous sections detailed how the fluxes and vertical profiles of absorbed energy within a multi-layered snowpack are computed. This section details how the single-scattering properties of snow, namely $\sigma_e$, $\omega$, and g, are determined from the snow physical properties, namely SSA, density, grain shape, and impurity contents. It is worth having in mind that the RTE applies to a continuous medium. As snow is a porous medium, it is common to define an optically equivalent continuous medium to represent it. In practice, its extinction coefficient is determined from the number concentration $N$ ($m^{-3}$) of snow grains and the average extinction cross section of snow grains (Kokhanovsky and Zege, 2004):

$$ \sigma_e = N \bar{C}_{\text{ext}}. \qquad (67) $$

This strategy, which was originally developed to compute for instance the optical properties of clouds (Stephens, 1978), implicitly assumes that scatterers are independent. Although this is unlikely to be the case in snow which is a dense medium, this formalism remains widely used since it has proved its efficiency in simulating snow optical properties in the solar spectrum where ice absorption is relatively low. Actually it remains efficient as long as the asymmetry factor is also computed assuming independent scatterers (Kokhanovsky, 2004). Using this representation, snow density is related to the average volume of snow grains $\bar{V}$: $\rho = N\rho_{\text{ice}}\bar{V}$. We further assume that snow grains are large compared to the wavelength of solar radiation, so that $C_{\text{ext}} = 2\Sigma$ (where $\Sigma$ is the projected area of an individual grain), and that the grains are convex so that $\Sigma = S/4$, where S is the total surface area of a grain. Hence, $\sigma_e$ finally reads

$$ \sigma_e = \frac{\rho \text{SSA}}{2}, \qquad (68) $$

where SSA equals by definition $\frac{S}{\rho_{\text{ice}}\bar{V}}$. Note that this expression was originally known for convex particles only (e.g., Libois et al., 2013) but was then applied to a more general porous medium (Malinka, 2014).

The single-scattering albedo is computed after Kokhanovsky and Macke (1997), who propose an analytical expression depending on the refractive index $m = n - i\chi$ and grain shape S, based on Monte Carlo computations relying on the geometrical optics approximation:

$$ (1 - \omega) = \frac{1}{2}(1 - W(n))(1 - e^{-\psi(n, S)c}), \qquad (69) $$
where

$$c = \frac{24\pi\chi}{\rho_{\text{ice}}\lambda\text{SSA}}, \qquad (70)$$

and $\rho_{\text{ice}}$ is the bulk density of ice. W does not depend on grain shape (for randomly oriented, convex, particles, Kokhanovsky and Macke, 1997) and is assumed to depend linearly on n based on tabulated values (Kokhanovsky, 2004, p. 61), so that in TARTES

$$W(n) = 0.0611 + 0.17(n - 1.3). \qquad (71)$$

Likewise,

$$\psi(n, S) = \frac{2}{3} \frac{B(n, S)}{1 - W(n)}, \qquad (72)$$

where $B(n, S)$ is the absorption enhancement parameter. Note that at low absorption Eq. (69) collapses to Eq. (6) of Libois et al. (2013). TARTES proposes three options to compute B. The first option (the only one in the previous TARTES version) is a linear dependence on n based on Kokhanovsky and Macke (1997), so that

$$B(n, S) = B_0(S) + 0.4(n - 1.3), \qquad (73)$$

where $B_0(S)$ can be prescribed by the user to account for a particular geometrical shape. For instance, for spherical particles $B_0 = 1.25$ (Libois et al., 2013).

The second option is

$$B(n, S) = n^2, \qquad (74)$$

which stems from the recent work on random media and is now the recommended and default option (Malinka, 2014; Robledano et al., 2023). The last option gives the user the possibility to prescribe a constant value or a value per wavelength.

The asymmetry factor g also depends on the detailed snow microstructure, but in the granular representation it can be computed from the scattering phase function of individual grains. At weakly absorbing wavelengths g mainly depends on snow grain shape S, but at absorbing wavelengths it also depends on the ice imaginary part of the refractive index $\chi$ and SSA. In TARTES g is calculated in consistency with B, according to the three options. In the first option $g(n, S)$ is computed after Kokhanovsky and Macke (1997):

$$g(n, S) = g_{\infty}(n) - [g_{\infty}(n) - g_0(n, S)]e^{-y(n, S)c}. \qquad (75)$$

$g_{\infty}(n)$ is the asymmetry factor of a purely absorbing sphere, and $g_0(n, S)$ is the asymmetry factor of the non-absorbing particle of shape S. $g_0$ and $g_{\infty}$ are both assumed to depend linearly on n, so that

$$g_{\infty}(n) = 0.9751 - 0.105(n - 1.3), \qquad (76)$$

$$g_0(n, S) = g_0(S) - 0.38(n - 1.3), \qquad (77)$$

where $g_0(S)$ is prescribed by the user. Again the dependence on n corresponds to that of spheres (Kokhanovsky, 2004). In

the second option, $g = 0.82$ (Robledano et al., 2023), and in the third option, the user prescribes the value as a constant or as one value per wavelength.

Finally y is also assumed to depend linearly on n. In TARTES the expression corresponding to spheres is taken so that (Kokhanovsky, 2004)

$$y(n) = 0.728 + 0.752(n - 1.3). \qquad (78)$$

Note that for the three options, the variables W, $g_{\infty}$, and y are calculated by linear relationships (Eqs. 71, 76, and 78) corresponding to spheres. This may result in inconsistencies with respect to B and g, but given that at present these three variables have not been investigated for snow, we prefer to keep the relationship used up to now.

## 2.10 Impurities

At the wavelengths where ice is very weakly absorbing, the optical properties of snow are very sensitive to the presence of light-absorbing impurities. TARTES can account for such impurities, in practice black carbon (BC), dust, and humic-like substances (HULIS). For the sake of simplicity it is assumed that impurities are external to snow grains. When impurities are added in realistic, low quantities, $N_i$, we assume that the extinction coefficient of snow is unchanged but the absorption coefficient is altered. This supposes that impurity scattering is negligible. According to simulations at 1030 nm (Fair et al., 2022), this applies to BC in any case, as well as to dust except for fine particles (< 1 µm) in high concentration (e.g., > 500 ppm). When this approximation is valid, it follows that the single-scattering co-albedo is

$$(1 - \omega) = (1 - \omega)_{\text{snow}} + \frac{1}{\sigma_e} \sum_i N_i C_{\text{abs}}^i, \qquad (79)$$

where $C_{\text{abs}}^i$ is the average absorption cross section of impurities i. Rewritten as a function of the bulk mass concentration $c_i$ (kg kg$^{-1}$) and using Eq. (69), it reads

$$(1 - \omega) = \frac{1}{2}(1 - W(n))(1 - e^{-\psi(n, s)c}) - \frac{2}{\lambda\text{SSA}} \sum_i \text{MAE}^i c_i, \qquad (80)$$

where MAE is the mass absorption efficiency (in m$^2$ kg$^{-1}$, e.g., Caponi et al., 2017).

To calculate MAE, TARTES uses different formulations according to the particle size. For small particles compared to the wavelength (applies to BC and HULIS), the absorption cross section $C_{\text{abs}}$ of an impurity of type i, of volume $V_i$ and with complex refractive index $m_i$, is given by Kokhanovsky (2004):

$$C_{\text{abs}}^i = - \frac{6\pi V_i}{\lambda} \text{Im} \left( \frac{m_i^2 - 1}{m_i^2 + 2} \right). \qquad (81)$$
Dividing by $\rho_i V_i$, where $\rho_i$ is the impurity bulk density, yields the mass absorption efficiency:

$$MAE^i = -\frac{6\pi}{\lambda\rho_i}\mathrm{Im}\left(\frac{m_i^2 - 1}{m_i^2 + 2}\right). \qquad (82)$$

In TARTES, the characteristics of BC can be taken from Bond and Bergstrom (2006) ($\rho_{BC} = 1800 \text{ kg m}^{-3}$ and $m_{BC} = 1.95 - 0.79i$) or by default from SNICAR-ADv3 (Flanner et al., 2021). The values for the HULIS are taken from Hoffer et al. (2006).

For large particles (applies to dust), the absorption is not simply related to the imaginary part of their refractive index and the volume but depends on the shape, size, and other impurity particularities. While Mie theory applies to spherical particles – and other more complex theories to more general shapes (Mishchenko et al., 1996) – the computation is usually intensive. Instead TARTES directly uses tabulated MAE values that can be obtained from independent calculations or from in situ measurements. More precisely, the MAE is calculated from the MAE at a specific wavelength (usually $\lambda_0 = 400 \text{ nm}$ or $\lambda_0 = 550 \text{ nm}$) and the spectral dependence given by the Ångström absorption exponent (AAE) such as

$$MAE(\lambda) = MAE(\lambda_0) \left(\frac{\lambda}{\lambda_0}\right)^{-AAE}. \qquad (83)$$

In TARTES v2.0, values from Caponi et al. (2017) obtained from different regions in the world are implemented. For small particles (PM$_{2.5}$), available locations are Libya, Morocco, Algeria, Mali, Saudi Arabia, Kuwait, Namibia, China, and Australia, and for large particles (PM$_{10}$), available locations are Libya, Algeria, Bodele, Saudi Arabia, Namibia, China, Arizona, Patagonia, and Australia (Caponi et al., 2017).

## 3 Numerical implementation

TARTES was initially implemented in Python and then converted in part to Fortran to be integrated in the detailed snowpack model Crocus (Vionnet et al., 2012), where it can be used to compute the profile of absorption in the snowpack (Libois et al., 2015; Tuzet et al., 2017). Both versions are based on the same equations, though the Python version has more features beyond absorption calculation. Note also that the Fortran version is not systematically updated along with the Python version.

### 3.1 Python version

TARTES v2.0 is compatible with Python 3.7 and higher. The four main entry points for the user are the functions `albedo`, `absorption_profile`, `irradiance_profile`, and `actinic_profile`. The inputs and outputs are listed in Table 1. These functions take as inputs the wavelengths at which the computation

is performed, the properties of each layer (thickness, SSA, density, $B_0$ and $g_0$, impurities), the name of the ice refractive index database (Warren and Brandt, 2008; Picard et al., 2016), the bottom albedo, and the illumination conditions (solar zenith angle, total radiation flux $F_0 + F_0^{\text{tot}} = F_0 + F_0^{\text{diff}}$, and fraction of direct radiation $F_0/F_0^{\text{tot}}$). The outputs are the albedo, the absorption, the upwelling and downwelling irradiances, and the actinic flux for each wavelength, for each function respectively. In addition, a function allows the calculation of broadband albedo given the incident spectrum distribution.

These user functions internally call the core function `tartes`, which computes the intrinsic optical properties of each layer from their physical properties and impurity contents, and then call the function `two_stream rt`, which solves the radiative transfer equation, eventually allowing us to compute all the quantities mentioned here above (Fig. 1).

The code has a test suite (10 tests currently) that can be automatically run using the software `pytest` to check the conformity of the code.

Compared to the previous version, TARTES v2.0 proposes several options to compute diffuse radiation and to compute the semi-infinite layer albedo that appears in the two-stream equation. The default values for the grain shape parameters are changed based on recent advances (Robledano et al., 2023) and likewise for the black carbon properties to match SNICAR-ADv3. The code has also been improved with the factorization of the impurity calculations (all types are now using MAE, either tabulated or calculated from the refractive index), type hinting, automatic strict formatting, modern packaging, and continuous integration for automatic publication on PyPI. The documentation has been improved, and the conformity between the equations presented in this paper and the code has been carefully checked. In addition, the SnowTARTES web application has been enhanced to calculate the irradiance profile.

### 3.2 Fortran version

The original TARTES code was converted in a Fortran version suitable for the integration in Crocus (Vionnet et al., 2012). Crocus predicts the evolution of a multi-layered snowpack from local meteorological forcings. For this the absorption of the solar energy in each layer needs to be computed at each time step. This is critical to compute the energy budget of the snowpack and consequently the temperature profile and gradients that control snow metamorphism (Flanner and Zender, 2005). The original optical scheme in Crocus (Brun et al., 1989) estimates the albedo, hence the total absorbed energy, in three spectral bands from empirical relationships with the grain shape and size, and then it applies an ad hoc Beer-Lambert law to distribute the absorbed energy in the layers. Using a proper radiative transfer model has many advantages: the higher spectral resolution (10 nm by default) allows us to resolve the spectral features of snow without
<table><caption>**Table 1.** Symbols, input parameters, and outputs of the user functions in TARTES. Parameters marked with (l) are given for each layer. The (z) dimension is given by the user and is independent of the layers.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Description</th><th>Default value</th><th>Unit/type</th></tr></thead><tbody><tr><td colspan="5">Inputs</td></tr><tr><td>wavelength</td><td>λ</td><td>Wavelength</td><td></td><td>m</td></tr><tr><td>SSA(l)</td><td>SSA</td><td>Specific surface area</td><td></td><td>m<sup>2</sup> kg<sup>-1</sup></td></tr><tr><td>density(l)</td><td>ρ</td><td>Density for each layer</td><td></td><td>kg m<sup>-3</sup></td></tr><tr><td>thickness(l)</td><td></td><td>Thickness of each layer</td><td></td><td>m</td></tr><tr><td>z(z)</td><td></td><td>Depth at which irradiance is calculated</td><td></td><td>m</td></tr><tr><td>g0(l)</td><td>g0</td><td>Asymmetry parameter (n = 1.3, non-absorbing)</td><td>0.82</td><td>–</td></tr><tr><td>B0(l)</td><td>B0</td><td>Absorption enhancement (non-absorbing)</td><td>n\, ^2</td><td>–</td></tr><tr><td>shape_parameterization(l)</td><td></td><td>Method to compute B(λ) and g(λ)</td><td>robledano23</td><td>string</td></tr><tr><td>impurities(l)</td><td>c<sub>i</sub></td><td>Impurity concentration in each layer</td><td></td><td>kg kg<sup>-1</sup></td></tr><tr><td>impurities_type(l)</td><td></td><td>Type of impurities</td><td>SootSNICAR3</td><td>string</td></tr><tr><td>refrac_index</td><td>n</td><td>Real and imaginary parts of refractive index</td><td>p2016<br>Picard et al. (2016)</td><td>value or data set</td></tr><tr><td>soilalbedo</td><td>α<sub>b</sub></td><td>Spectral albedo of underlying layer</td><td>0</td><td>–</td></tr><tr><td>dir Frac</td><td>F<sub>0</sub>/F<sub>0</sub><sup>tot</sup></td><td>Fraction of directional flux</td><td></td><td>–</td></tr><tr><td>totflux</td><td>F<sub>0</sub><sup>tot</sup></td><td>Total spectral incident flux (direct + diffuse)</td><td></td><td>any unit</td></tr><tr><td>diff_method</td><td></td><td>Diffuse radiation compute method</td><td>aart eq<br>48.2°</td><td>string</td></tr><tr><td>infprop_method</td><td></td><td>Infinite property calculation method</td><td>delta_eddington</td><td>string</td></tr><tr><td>sza</td><td>θ<sub>0</sub></td><td>Solar zenith angle</td><td></td><td>degrees</td></tr><tr><td colspan="5">Outputs</td></tr><tr><td>albedo</td><td>α</td><td>Spectral albedo</td><td></td><td>–</td></tr><tr><td>broadband albedo</td><td>α&#773;</td><td>Broadband albedo</td><td></td><td>–</td></tr><tr><td>absorption(l)</td><td>E<sub>i</sub></td><td>Spectral absorption</td><td></td><td>as F<sub>0</sub><sup>tot</sup></td></tr><tr><td>irradiance(z)</td><td>F<sub>tot</sub><sup>+</sup>, F<sub>tot</sub><sup>-</sup></td><td>Spectral irradiance</td><td></td><td>as F<sub>0</sub><sup>tot</sup></td></tr></tbody></table>

resorting to questionable spectral averages; it computes the profile of energy in a way that is consistent with the albedo; the direct and diffuse fluxes and the solar zenith angle dependence are treated properly; and absorption by the soil, snow grain shape, and impurities are fully accounted for. However, the computation time can be significant, especially if the spectral integration is highly resolved. This is the reason why studies have developed optimized wavelength sampling strategies to reduce the number of computations needed (van Dalum et al., 2019; Veillon et al., 2021).

Although the translation of TARTES in Fortran was initially motivated by the integration in Crocus, it is a self-contained model that can be integrated into any other model. The source code is part of SURFEX (http://www.umr-cnrm.fr/surfex/, last access: 1 March 2024, Masson et al., 2013). We call this version TARTES.F hereafter.

TARTES.F was developed in 2014 exactly following the Python version. It contained the same physics and parameters. However, it has not been updated since then and now slightly differs from the most recent Python version. The main difference is the ice refractive database that is based on Warren and Brandt (2008) only (Picard et al., 2016, was not available). In 2019, impurities were added (Tuzet et al., 2017) using specific MAE values that slightly differ from the Python version (not used here).

In terms of performance, a simulation for a two-layer snowpack at 106 wavelengths repeated 10 000 times takes 3.4 s with TARTES.F on a commodity laptop, while the Python version takes 148 s, about 40 times longer.
## 3.3 Related software and libraries for snow optics

### 3.3.1 SnowTARTES web application

SnowTARTES is an interactive web application (https://snow.univ-grenoble-alpes.fr/snowtartes, last access: 24 March 2024) meant to easily compute spectral albedo and irradiance profiles using TARTES without writing code. SnowTARTES uses the Python version and thus provides exactly the same results. The snowpack is described layer by layer in a text form. In order to conduct a sensitivity analysis, any of the input parameters can be prescribed as a range (start, end, step) instead of a single value. These ranges are combined, launching multiple calculations (limited to a maximum of 10 for the sake of visibility in the plot). The calculated albedo spectra are immediately plotted (one curve for each calculation) and can be downloaded as a comma-separated-value-formatted file.

SnowTARTES offers a selection of grain shapes to chose from, each one corresponding to a pair ($B_0, g_0$) based on Li-bois et al. (2013) and Robledano et al. (2023). Similarly, it offers a selection of seven background albedo spectra (grass, ice, several soils, along with 0 and 1) extracted from the Johns Hopkins University Spectral Library (ECOSTRESS and formally ASTER library, https://speclib.jpl.nasa.gov/, last access: 24 March 2024).

### 3.3.2 Snowoptics Python package

The Python package Snowoptics (https://github.com/ghislainp/snowoptics, last access: 24 March 2024) provides a series of functions to compute albedo, extinction, and BRDF for a semi-infinite homogeneous snowpack using the AART theory (Kokhanovsky and Zege, 2004; Kokhanovsky, 2012). The arguments of these functions are similar (name and unit) to TARTES. The package also provides functions to compute the effect of the slope on measured albedo and to correct from the slope effect based on Picard et al. (2019), which is not available in TARTES. Here, Snowoptics is used in Sect. 4.2 to compare TARTES with AART.

### 3.3.3 AtmosRT Python package

Because of the impact of the illumination geometry on the spectral albedo and of the spectral distribution of incident ir-radiance on the broadband albedo, it is necessary to know the spectral solar irradiance. There are many available ra-diative transfer models for the atmosphere that can provide this information for TARTES calculations, such as SMARTS (Gueymard, 2001), MODTRAN (Berk et al., 2014), or li-bRadtran (Emde et al., 2016), but the offer in Python is lim-ited. PyRTM is an unmaintained software package (https://github.com/Queens-Applied-Sustainability/PyRTM, last access: 24 March 2024) providing a Python 2 interface to ac-cess two general atmospheric radiative transfer models writ-ten in Fortran, namely SBDART (Ricchiazzi et al., 1998)

and SMARTS. From this, we developed the AtmosRT pack-age, including support for Python 3 and a few additional minor improvements for ease of use. The function *atmo-spheric_incident_spectrum* is implemented in TARTES to perform simple calculations with SBDART through At-mosRT and directly provides the total flux and the direct fraction at each wavelength as required by the TARTES func-tions. The function takes the solar zenith angle and cloud op-tical depth as input and uses the default cloud optical proper-ties of SBDART (Ricchiazzi et al., 1998).

## 4 Results

Several simulations with TARTES and other models are com-pared in this section. By default, unless specified, we con-sider a semi-infinite homogeneous snowpack, i.e., made of a single layer, thick enough so that the bottom boundary does not influence the albedo and the presented profiles. The layer has an SSA of $20 \text{ m}^2\text{kg}^{-1}$, and the density is $350 \text{ kgm}^{-3}$, but this latter variable has no influence on thick-snow albedo in the conventional radiative transfer framework used in TARTES (Malinka, 2023). Other particular condi-tions of the simulations are indicated for each case.

### 4.1 Self-consistency checks

TARTES has a few options to select between different ap-proximations or different modes of calculations. Here, we compare these approximations and check their consistency.

#### 4.1.1 Diffuse illumination

The three methods to take into account diffuse illumination in TARTES are compared in Fig. 2 for the default TARTES version and in Fig. 3 for the TARTES∞AART version, for the spectral and broadband albedo of the default semi-infinite snowpack under typical illumination for wintertime alpine snow. The reference method uses integration (denoted by “Integr.”) of *n* simulations in direct illumination mode with a solar zenith angle varying from 0 to 90° (with *n* = 128 regular steps in cosine of the angle). The results show a close agreement with the direct calculations at a zenith an-gle of 48.2°, the equivalent angle obtained by matching the diffuse and direct expressions of AART albedo (Sect. 2.7). The latter approximation yields relatively accurate albedo, with a deviation lower than 0.003 at wavelengths < 1300 nm and never exceeding 0.007 in the investigated wavelength range (up to 2000 nm). The broadband albedo difference is virtually zero (0.0006). Similar results are obtained for the TARTES∞AART variant. The maximum error is 0.011 in this case, and the broadband albedo difference is 0.001.

In contrast, the formulation using diffuse illumination for the upper boundary condition of the two-stream ap-proximation (2S method) performs poorly in the original TARTES model, with deviations reaching nearly 0.02 at
wavelengths < 1300 nm and barely acceptable (> 0.03) be-
yond 1400 nm, although the broadband albedo difference
of 0.0045 probably remains small enough for most appli-
cations. The inherent difficulty in handling diffuse radia-
tion with the δ-Eddington approximation was discussed ex-
tensively by Wiscombe (1977) and Wiscombe and Warren
(1980). Conversely, the 2S method works well with the
TARTES∞AART as it gives nearly equivalent results to
the direct simulation at 48.2° (the dotted curves overlap in
Fig. 3b). This agreement is expected as TARTES∞AART
was designed to be equivalent to AART, and 48.2° is the
equivalent angle for the direct and diffuse illumination in this
theory.

As a conclusion, we suggest not to use the 2S method for practical applications, since it provides almost no benefit compared to the 48.2° method (at best an agreement, at worst a large error) over the spectral range of interest, and the calculation is only slightly faster than for a direct beam calculation (23 ms vs. 34 ms for the computation in Fig. 2 involving 320 wavelengths).

For the profiles of irradiance, absorption, and actinic flux,
only two methods are implemented: the integration and the
direct 48.2° calculation. Figure 4 shows the profile of ab-
sorption close to the surface of the semi-infinite snowpack at
800 nm, a wavelength maximizing absorption since snow co-
albedo (the proportion of absorbed radiation) is greater than
at shorter wavelengths, and the incoming solar radiation is
still relatively large compared to longer wavelengths. Only
the original TARTES version is used. At this wavelength,
90 % of the absorption occurs in the topmost ≈ 4 cm (Fig. 4a)
for the snowpack considered here. The profiles of absorption
obtained with both methods look similar, with maximum dif-
ferences of about 1.5 %, reaching close to the surface, as seen
in Fig. 4b.

Based on these results and calculations with a shallower
snowpack (not shown), the direct 48.2° calculation was cho-
sen as the default method to simulate diffuse radiation in
TARTES. The integration is in principle more accurate but
requires many more computations (solving the linear system
for 128 angles instead of 1) even though measuring the exe-
cution time (51 ms instead of 34 ms) does not show a differ-
ence in the same proportion because only the constant vector
of the linear system depends on the angle, not the matrix. In
practice, users who prefer the accuracy offered by the inte-
gration method can explicitly set this option.

4.1.2 Consistency between albedo, profile of irradiance, and profile of absorption calculations

As the albedo, irradiance, and absorption profiles in snow and
the absorption below the snowpack are computed by three
distinct Python functions, it is worth checking that energy
is conserved across these different quantities. To this end
we first compare the albedo with the absorption profile. The
semi-infinite snowpack is split in numerical layers of 1 cm

(top first meter), even though their properties are identical
(SSA = 20 m² kg⁻¹, density = 350 m² kg⁻¹). Comparing the
sum of all layer absorption A (divided by the incident irradi-
ance, set arbitrarily to 1 W m⁻²) with the co-albedo 1-α, we
found a residual numerical error < 2 × 10⁻¹⁶ for all wave-
lengths, which is close to the machine 64 bit floating point
limit. Likewise, we checked that the calculated albedo per-
fectly matches the ratio of upwelling and downwelling irra-
diance at the surface, calculated from the vertical profiles of
irradiance. These two tests are part of the automatic test suite
available in the TARTES code base.

4.2 Comparison of TARTES with the asymptotic
analytical radiative transfer (AART)

The comparison between TARTES and the AART to simu-
late the diffuse and direct albedos of a semi-infinite snow-
pack (a single thick homogeneous layer) is presented in
Figs. 5, 6, and 7.

Used with the same parameters (SSA, density, B₀, and g₀),
the comparison between AART and TARTES∞AART for
both the diffuse and direct illuminations demonstrates that
the formulation in Eqs. (40)–(41) and (46)–(47) allows our
code to conform to the AART analytical expression in the
case of a semi-infinite homogeneous snowpack. The maxi-
mum error is indeed 4 × 10⁻¹⁶, corresponding to numerical
rounding errors. This new formulation can be useful in cases
where the conformity with AART is essential.

AART is also indistinguishably similar to the original
TARTES model in the visible range and up to 1400 nm on
the spectrum plots in Figs. 5a and 6a. However, the resid-
uals in Figs. 5b and 6b highlight that the differences are
much larger than numerical rounding errors. Nonetheless,
they remain < 0.015 up to 1400 nm. At longer wavelengths
the differences become noticeable and increase up to around
0.024. Interestingly, the differences between the models are
mainly significant in the domain of strong ice absorption,
where none of these models is expected to be valid. Indeed,
AART is meant to be valid only for weak absorption (see
Fig. 8 of Kokhanovsky and Zege, 2004), and TARTES uses
the δ-Eddington approximation, which likewise is only rele-
vant for weak absorption (Wiscombe, 1977). Depending on
the application, the reported errors in the longer wavelengths
of the solar spectrum can be either negligible (e.g., broad-
band albedo calculations, absorption calculation) or major
(e.g., spectroscopic applications at 1550 nm as in Gallet et al.,
2009).

To further explore the difference between the models,
Fig. 7 shows the albedo at 1300 nm as a function of SSA,
B₀, and g₀. Again we observe the perfect similitude between
AART and TARTES∞AART. On the other hand, the differ-
ence with the original TARTES model is small except for
very low SSA, where it reaches a maximum absolute value of
0.015. This error is small in absolute albedo but corresponds
to a significant relative error considering that albedo is < 0.2.
<figure><img src="image_5.png" /></figure>

**Figure 2.** Comparison of albedo under diffuse illumination computed by three methods: (1) by integrating over all incident angles, (2) using a unique direct beam at the equivalent angle predicted by the AART (48.2°), and (3) using diffuse radiation for the upper boundary condition in the two-stream formulation (abbreviated 2S). The snowpack is semi-infinite with SSA = 20 m² kg⁻¹. The broadband albedo (ω) is indicated in the legend for each spectrum.

<figure><img src="image_6.png" /></figure>

**Figure 3.** Same as Fig. 2 but for the TARTES∞AART variant.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Comparison of the absorption profile under diffuse illumination computed by two methods: (1) by integrating over all incident angles and (2) using a unique direct beam at the equivalent angle predicted by the AART (48.2°). The snowpack is semi-infinite with SSA = 20 m² kg⁻¹.</figcaption></figure>

The differences when $B_0$ and $g_0$ are varied are even weaker over the investigated range. Exploring other wavelengths (results not shown) indicates that the difference is much smaller in the visible range and progressively increases in the near-infrared range. These results obtained for the SSA and $B_0$, along with the wavelength dependence, again confirm that the models mainly diverge when the single-scattering albedo is low (it is lower for large grains than for smaller ones), which is when the absorption is strong.

To provide a more general recommendation to future users, we explored a wide range of usual conditions ($\lambda$ < 1400 nm, SSA in 5–100 m² kg⁻¹, and SZA in 0–70°) and found a maximum difference of 0.025 (at SSA = 5 m² kg⁻¹, SZA = 40°, and $\lambda$ = 1400 nm). Furthermore, 90% of the simulations show a very small difference < 0.002, leading to the conclusion that, for a semi-infinite snowpack, TARTES is virtually equivalent to AART in most usual conditions.

## 4.3 Comparison of TARTES with other numerical snow radiative transfer models

TARTES is now compared to two widely used models: DISORT-Mie (with 16 streams) and SNICAR-ADv3 (online tool, http://snow.engin.umich.edu/, last access: 6 December 2024). Since the first model (and the second until a recent version) is limited to spherical particles, we consider this shape for all the simulations. In TARTES, this is achieved by letting $B_0$ = 1.25 and $g_0$ = 0.895 (Libois et al., 2013). The conditions of simulations are as similar as possible among the models. For instance, the same ice refractive index (Picard et al., 2016) is used for the three models. TARTES.F is also included in this comparison, but with a different ice refractive index (Warren and Brandt, 2008) as discussed below.

## 4.4 Clean semi-infinite snowpack

Figure 8 shows diffuse albedo simulations for the different models for a semi-infinite snowpack with SSA = 20 m² kg⁻¹. Overall the agreement is very good, with virtually unnoticeable differences in Fig. 8a, except for TARTES∞AART, which stands out for wavelengths higher than 1400 nm. The residual albedo panel (Fig. 8b) reveals small differences of around 0.01 and occasionally up to 0.03 in amplitude, which may be significant for some applications. From this comparison no outliers or particularly similar models emerge (except the new TARTES∞AART formulation). Furthermore, the presence of spikes and oscillations suggests numerical issues rather than physical differences or bugs. This is also suggested by the differences between the Python and Fortran versions of TARTES, despite a common theory and initial code. On average, TARTES appears closest to DISORT-Mie with a root mean square difference (RMSD) of 0.0035, followed by TARTES.F with a RMSD = 0.0039 and SNICAR-ADv3 with a RMSD = 0.0061. The agreement is overall much better for the shorter wavelengths < 1400 nm as also noted for the comparison with AART. The differences in broadband albedo reported in Fig. 8a are very small.

For $\lambda$ < 490 nm, DISORT-Mie yields an albedo of 1.0, which is certainly a rounding error. In this highly reflective domain, the 32 bit float arithmetic used by DISORT-Mie is certainly insufficient. The difference observed in the visible range for TARTES.F is explained by the use of the refractive index database (Warren and Brandt, 2008) (hard-coded in the Fortran code) that has been recently updated with stronger absorption values (Picard et al., 2016) used for the other model simulations.

## 4.5 Clean two-layer snowpack

Figure 9 shows the comparison for a typical two-layer snowpack made of a thin layer of fresh snow (a 1 cm thick layer with SSA = 50 m² kg⁻¹ and density 150 kg m⁻³) on top of slightly aged snow (an infinitely thick layer with SSA = 20 m² kg⁻¹). SNICAR-ADv3 is not included because the online version used in this paper does not handle multiple layers. Overall the results are similar to the one-layer case. The maximum difference is about 0.03, the average difference is much smaller, and the errors share some similar patterns with the former comparison.

As for the single-layered snowpack, TARTES∞AART stands out in the longer wavelengths. It was mainly implemented to check the conformity of TARTES with AART from a theoretical point of view but appears to provide no practical benefit over the original TARTES model. For this reason, we do not further consider it. The original TARTES model is kept as the default in the following and in the code.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of the diffuse albedo spectra computed by TARTES (with 48.2°), TARTES∞AART (with 48.2°), and standalone AART under diffuse illumination. Panel (a) shows the albedo and panel (b) the difference with respect to TARTES. The snowpack is semi-infinite with SSA = 20 m² kg⁻¹.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Same as Fig. 5 but for direct illumination (SZA = 60°).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison of direct albedo (SZA = 60°) at 1300 nm as a function of SSA, $B_0$, and $g_0$ computed by TARTES and AART. The snowpack is semi-infinite with SSA = 20 m² kg⁻¹.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of diffuse albedo spectra calculated by different numerical models for a semi-infinite snowpack with SSA = 20 m² kg⁻¹.</figcaption></figure>

### 4.5.1 Clean thin snowpack

The extreme case of a 1 cm thick snowpack with a perfectly black underlying surface (bottom albedo $\alpha_B$= 0) is presented in Fig. 10. In principle, the two-stream approximation is less adequate in these conditions (corresponding to $\tau^* = 7–17$ depending on the wavelength) compared to the DISORT-Mie model. The results indeed show a degradation in the visi-

ble range, where the underlying surface has an impact (the albedo is notably lower than in Fig. 8). Nevertheless, the difference remains under 0.01 in this range, which is comparable to the more favorable case of the thick snowpack.

### 4.6 Snowpack polluted with black carbon and dust

Light-absorbing impurities can be taken into account by all models considered here. Here we compare simulations with
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Comparison of diffuse albedo spectra calculated by different numerical models for a two-layer snowpack. A 1 cm layer with $SSA = 50 \text{ m}^2 \text{ kg}^{-1}$ and density $= 150 \text{ kg m}^{-3}$ is overlying a semi-infinite layer with $SSA = 20 \text{ m}^2 \text{ kg}^{-1}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison of diffuse albedo spectra calculated by different numerical models for a thin snowpack (1 cm thick) with $SSA = 20 \text{ m}^2 \text{ kg}^{-1}$ and density $= 350 \text{ kg m}^{-3}$ overlying a dark surface ($\alpha_b = 0$).</figcaption></figure>
particles much smaller than the wavelength (BC) first and then much larger than the wavelength (dust).

Figure 11 shows semi-infinite diffuse albedo spectra obtained by the three models for BC concentrations of 100, 500, and 2000 ng g⁻¹. For the sake of testing and validation, we use such extreme values compared to the typical amount found in snow (Bisiaux et al., 2012; Kang et al., 2020). Contrary to the clean snowpack, significant differences are obtained for concentrations of 500 and 2000 ng g⁻¹. The difference reaches 0.04 around 500 nm between DISORT-Mie and TARTES at the maximal concentration. Such a value, in the visible range where the solar irradiance is maximum, has a dramatic effect on the absorption. This is reflected in the differences in broadband albedo, which are limited to 0.001–0.006 between SNICAR-ADv3 and TARTES but reach 0.005–0.022 between DISORT-Mie and TARTES depending on the concentration.

These differences are likely explained by the different representations of the impurities in these models. For maximal comparability, all the simulations assume the same refractive index and density of BC. They are taken from SNICAR-ADv3 (Flanner et al., 2021), which is also the default in TARTES v2.0. The three models also assume that the particles are suspended in the air (a representation known as external mixing, Flanner et al., 2012). However, the size of the particles differs for reasons inherent to each model. SNICAR-ADv3 and DISORT-Mie use Mie theory (spherical particles); however, SNICAR-ADv3 simulates a lognormal distribution of particle sizes, with a mean mass-weighted radius of 67 nm for black carbon, while our implementation of DISORT-Mie is limited to a mono-disperse collection of spheres. A radius of 67 µm is taken for Fig. 11, but this is not strictly equivalent to the SNICAR-ADv3 configuration. On the other hand, TARTES uses the Rayleigh approximation that is only valid for small particles. Comparing MAE predicted by TARTES and reported for SNICAR-ADv3 in Fig. 3a in Flanner et al. (2021), we observe (not shown) a slightly higher MAE with TARTES at 400 nm (+20 %), an agreement at 490 nm, and a slightly lower MAE at 1000 nm (−23 %), which explains the lower TARTES albedo at wavelengths < 500 nm and higher at longer wavelengths observed in Fig. 11.

To illustrate the influence of the particle size, Fig. 12 shows the variations of albedo at 550 nm predicted by DISORT-Mie. TARTES compares very well with DISORT-Mie for the lowest radii (blue marker overlapping the violet curve in Fig. 12), while SNICAR-ADv3 slightly overestimates the DISORT-Mie albedo at 67 µm (green marker). This results from the fact that SNICAR simulations with BC mean radius of 67 µm include actual impurities with radii smaller and larger than 67 µm (it assumes a lognormal size distribution) for which the absorption is larger than at 67, as 67 µm appears to be very close to the minimum absorption. Nevertheless, the variations of albedo in the range 0–150 nm are relatively small (< 0.02), which seems acceptable with respect to the other uncertainties on the properties of the light-absorbing particles (concentration and density).

Figure 13 shows a comparison of different dust types (different origins and sizes) and models. The concentration is 100 µg g⁻¹, which is realistic in alpine snow (Dumont et al., 2020; Di Mauro et al., 2024). The agreement is again fairly good between TARTES and SNICAR-ADv3, which results from the similar MAE values for the Libya PM₂.₅ and Algeria PM₂.₅ dusts (400 nm = 110 and 73 m² kg⁻¹) used in TARTES, compared to the values for the Sahara used in SNICAR-ADv3 (Flanner et al., 2021, Table 2, Fig. 3).

## 4.7 Simulations of profiles of absorption, irradiance, and actinic flux.

Figure 14 illustrates TARTES ability to calculate profiles of absorption, irradiance, and actinic flux for a three-layer snowpack with SSA = 50, 20, and 20 m² kg⁻¹ from top to bottom; density = 300, 300, and 350 kg m⁻³; and thickness = 10 and 20 cm, with the last layer being infinitely thick. The wavelength is 600 nm. The profiles are presented relative to the incident flux; this is why the x-axis label of each graph is unitless.

Regarding the absorption profile, the convention in TARTES is to return the total radiation absorbed in every layer which is suitable for a direct input in thermodynamic calculations. The unit is the same as that of the incident flux prescribed by the user (variable totflux, usually W m⁻² nm⁻¹ in real applications). Regarding the actinic flux, the convention in TARTES is to return values in the same unit as the incident flux, which is prescribed by the user in the variable totflux. The conversion from spectral irradiance W m⁻² nm⁻¹ to actinic flux photons s⁻¹ cm⁻² nm⁻¹ is left to the user. The irradiance and actinic flux profiles (Figs. 14b, c) show a series of near-linear decreasing trends in logarithm scale with varying slopes, which is equivalent to near-exponential decreases in the natural scale with the varying rate. This rate is approximately the asymptotic extinction (Libois et al., 2013), which can be deduced by approximating the AART extinction Eq. (41) as

$$k_e = \rho \sqrt{\frac{12B \chi \text{SSA}(1 - g^G)}{4\lambda \rho_{\text{ice}}}}. \quad (84)$$

This formulation is implemented in the Snowoptics package, function extinction_KZ04. Here we find $k_e = 9.14, 5.78$, and $8.09 m^{-1}$ for the layers from top to bottom, equivalent to e-folding depths ($l_e = 1/k_e$) of 11, 17, and 12 cm respectively. Calculating the vertical gradient of the irradiance logarithm from Fig. 14b (excluding 1 cm at the top and bottom of each layer) yields similar values: $k_e = 11.1, 5.78$, and $8.09$ cm respectively.

The behavior of the irradiance gradient in the top layer is affected by the proximity of the surface (where the direct radiation is progressively converted into diffuse radiation)
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Comparison of diffuse albedo calculated by different numerical models for different BC concentration. The snowpack is semi-infinite, SSA = 20 m² kg⁻¹, and density = 350 kg m⁻³.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Diffuse albedo at 550 nm as a function of BC particle radius, calculated by different numerical models. The snowpack is semi-infinite, BC concentration = 500 ng g⁻¹, SSA = 20 m² kg⁻¹, and density = 350 kg m⁻³.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Comparison of diffuse albedo spectra calculated for different dust types and sizes by two models. The snowpack is semi-infinite, dust concentration = 100 µg g⁻¹ for all types, SSA = 20 m² kg⁻¹, and density = 350 kg m⁻³.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Profiles of (a) absorption, (b) irradiance, and (c) actinic flux for a three-layer snowpack with SSA = 50, 20, and 20 m² kg⁻¹; density = 250, 250, and 350 kg m⁻³; and thickness = 5 cm, 20 cm, and ∞. The illumination is at 60° and the wavelength is 600 nm.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Spectral albedo (a) and irradiance at 20 cm depth, relative to the incident irradiance at the surface (b) for spheres and snow (i.e., $B = n^2$ and $g \in [0.8, 0.835]$ according to Robledano et al., 2023). The snowpack is semi-infinite with SSA = 20 m² kg⁻¹, and the illumination is at 60°.</figcaption></figure>

and the junction with the next layer that has different optical properties. This tends to bend the curve, meaning that the profile of irradiance is not exactly exponential. The actinic flux shows a similar behavior.

To conclude, this example illustrates that the profiles of irradiance and actinic flux can be approximated at first order by exponential decreases whose decay can be calculated from snow properties (density, SSA, $B_0$, and $g_0$) in each layer (Eq. 84)). However, near the surface and in the presence of contrasted layers, it is recommended to use a proper multi-layered radiative transfer model such as TARTES.

# 5 Discussion

TARTES was designed to perform simple radiative transfer calculations in a plane-parallel multi-layered snowpack, with the unique possibility of describing the shape of the particles using two parameters also used in the AART theory, namely the absorption enhancement parameter *B* and the asymmetry factor *g* (Libois et al., 2013). This choice is important and motivated by two reasons. First these parameters are the main factors controlling the influence of the shape on the absorption and scattering properties in a weakly ab-
sorbing medium. Second, these parameters can be calculated for two-phase porous media, without assuming that snow is a collection of particles with some given geometrical shape (Malinka, 2014; Robledano et al., 2023). While the physical meaning of *B* and *g* is often presented for single particles (e.g., Libois et al., 2013), the definition of these parameters is independent of the concept of particles. For this reason we qualify these parameters as descriptive of the “optical shape of snow” without requiring the medium to be actually composed of distinct particles (spheres, fractals, or cubes). Moreover, it was found that $B = n^2$ applies very well to snow (Robledano et al., 2023) and only *g* varies, but in a narrow range for snow compared to across common geometrical shapes (spheres, fractals, hexagonal plates). Furthermore, the value of $B = n^2$ and the values of *g* clearly indicate that snow does not behave as a collection of ice spheres. These results make TARTES inherently more suitable to snow than Mie-based models. Note that most simulations presented in this paper used the values of *B* and *g* for spheres for the sole purpose of comparison with the established Mie-based models. In practice, we do not recommend running TARTES in these conditions. Instead we recommend $B = n^2$ and $g = 0.82$ (at non-absorbed wavelengths), which is the middle of the range found by Robledano et al. (2023). These values are the defaults in TARTES v2.0. Figure 15 illustrates the difference in albedo and irradiance at depth (20 cm) when considering the default *B* and *g* values for snow or the values for spheres. The albedo is higher for snow than for spheres for a given SSA and density by 0.018 on average over the range 400–2000 nm and reaches 0.042 at 1400 nm. These values are significant for surface energy budget calculations, with a potent large impact in climate simulations (Räisänen et al., 2017). The irradiance at 20 cm depth is weaker for snow than for spheres, by about a factor 10 at 750 nm for instance. These differences can be explained by the strong forward scattering of spheres (high *g*) and the lower absorption enhancement parameter (low *B*), which tends to overestimate the penetration depth.

Despite this advantage, TARTES presents some limitations owing to its simplicity. It uses the conventional unpolarized radiative transfer, neglecting interferences, near-field and packing effects, and polarization effects. As a planet-parallel model, the surface is supposed to be perfectly flat, and surface roughness is neglected, which may impact simulations on rough terrain, especially at grazing angles. Likewise, the layers are perfectly smooth and horizontally semi-infinite. In practice, in areas where horizontal heterogeneity is strong, for instance as a result of snowdrift, this assumption might be inappropriate to simulate snow optical properties. Also, TARTES only considers snow excluding any other material that might be present in the snowpack and models its optical behavior as a homogeneous scattering medium. It means that the layers must be much thicker than the grain size (i.e., the free photon path). Regarding impurities, only their absorption is considered, and they are ran-

domly distributed. Consistent with the choice of only using the asymmetry factor *g* instead of requiring the full phase function, the two-stream approximation was selected to solve the radiative transfer equation in TARTES instead of a multi-stream approach as in DISORT-based models. As a direct consequence, TARTES cannot calculate the bidirectional reflectance distribution function (BRDF), which is essential for instance for satellite remote sensing applications. TARTES was initially designed for energy balance computations and only provides hemispherically averaged quantities, namely surface albedo, absorption in each layer, and profiles of upwelling and downwelling radiation flux.

As TARTES relies on the geometrical optics approximation, the allowed range of ice particle sizes (or more generally, the length scales in the microstructure) is limited. The particles must be significantly larger than the wavelength (typically > 5 µm in the solar range), which is usually valid for most snow types (Fierz et al., 2009; Walden et al., 2003). In contrast, this assumption is invalid for light-absorbing impurities such as BC that usually come as sub-wavelength-sized particles. It is still possible to account for the absorption of these small particles by neglecting their scattering. For very small particles, the Rayleigh approximation works well and allows a rigorous formulation of the absorption coefficient as a function of the particle complex refractive index and density. This approximation is acceptable for BC. For larger particles, such as dust, TARTES relies on tabulated values of mass absorption efficiency (MAE), which either can be obtained from direct measurements or can be estimated by offline Mie calculations as in SNICAR. We believe that this simple treatment of light-absorbing properties in snow is sufficient and of adequate complexity given the considerable uncertainties associated with the physical properties of light-absorbing impurities and the difficulty in measuring or simulating their concentration in snow.

Despite the differences between TARTES and the other models, the simulations of albedo presented in this paper (Sect. 4) generally show a good agreement, typically within 0.02. Notably, the errors between different models are typically lower than between the Python and FORTRAN versions of TARTES, which suggests that at this degree of agreement, most of the residual errors can result from implementation details and numerical issues rather than theoretical differences.

Possible future improvements in TARTES include the inclusion of terrain slopes, bubbly ice and slush layers, the extended impurity database, and internal mixture impurities.

## 6 Conclusions

TARTES and the ecosystem of tools developed around this radiative transfer model for snow allow accurate simulations of several snow optical properties, most notably the spectral albedo and irradiance profiles in the snowpack. TARTES is
intended to be user-friendly and easy to improve, thanks to the Python implementation. While technically the results of this paper demonstrate that TARTES performs equally well compared to other existing models when assuming snow as a collection of ice spheres, TARTES can handle a more general representation of snow that is more representative of natural snow than historical models based on idealized shapes to represent snow grains. For this particular reason, and despite the overall simplicity of the approximations implemented in TARTES, this model is able to accurately predict the optical properties of snow with a given SSA and density and is perfectly suited for implementation in atmospheric models, including climate models.

*Code and data availability.* The TARTES v2.0 model used in this study is available from https://doi.org/10.5281/zenodo.13950598 (Picardand Libois, 2024). The SnowTARTES web app is accessible from https://snow.univ-grenoble-alpes.fr/snowtartes/ (last access: 6 December 2024; Picard, 2021). No data sets were used in this article.

*Author contributions.* QL wrote the theoretical formulation of TARTES. QL and GP developed the TARTES Python version. GP ran the simulations and wrote the manuscript. Both authors commented on the manuscript.

*Competing interests.* The contact author has declared that neither of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

*Acknowledgements.* We thank Matthieu Lafaysse for providing the FORTRAN version that he implemented from the Python version of TARTES for integration in the open-source model Crocus. The authors are thankful to Laurent Arnaud, Marie Dumont, Alvaro Robledano, and François Tuzet for fruitful discussions on the model, as well as to Zhuang Jiang for detecting a last-minute bug in the actinic code.

*Financial support.* This research has been supported by the Agence Nationale de la Recherche project MiMESis-3D (grant no. ANR-19-CE01-0009).

*Review statement.* This paper was edited by Xiaohong Liu and reviewed by Cenlin He and Mark Flanner.

## References

Aoki, T., Kuchiki, K., Niwano, M., Kodama, Y., Hosaka, M., and Tanaka, T.: Physically based snow albedo model for calculating broadband albedos and the solar heating profile in snowpack for general circulation models, J. Geophys. Res., 116, D11114, https://doi.org/10.1029/2010jd015507, 2011.

Berk, A., Conforti, P., Kennett, R., Perkins, T., Hawes, F., and van den Bosch, J.: MODTRAN® 6: A major upgrade of the MODTRAN® radiative transfer code, in: 2014 6th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), Lausanne, Switzerland, 24–27 June 2014, https://doi.org/10.1109/whispers.2014.8077573, 2014.

Bisiaux, M. M., Edwards, R., McConnell, J. R., Curran, M. A. J., Van Ommen, T. D., Smith, A. M., Neumann, T. A., Pasteris, D. R., Penner, J. E., and Taylor, K.: Changes in black carbon deposition to Antarctica from two high-resolution ice core records, 1850–2000 AD, Atmos. Chem. Phys., 12, 4107–4115, https://doi.org/10.5194/acp-12-4107-2012, 2012.

Bohren, C. F.: Multiple scattering of light and some of its observable consequences, Am. J. Phys., 55, 524–533, https://doi.org/10.1119/1.15109, 1987.

Bohren, C. F. and Barkstrom, B. R.: Theory of the optical properties of snow, J. Geophys. Res., 79, 4527–4535, https://doi.org/10.1029/jc079i030p04527, 1974.

Bond, T. C. and Bergstrom, R. W.: Light Absorption by Carbonaceous Particles: An Investigative Review, Aerosol Sci. Tech., 40, 27–67, https://doi.org/10.1080/02786820500421521, 2006.

Box, J. E., Wehrlé, A., van As, D., Fausto, R. S., Kjeldsen, K. K., Dachauer, A., Ahlstrøm, A. P., and Picard, G.: Greenland Ice Sheet Rainfall, Heat and Albedo Feedback Impacts From the Mid-August 2021 Atmospheric River, Geophys. Res. Lett., 49, e2021GL097356, https://doi.org/10.1029/2021gl097356, 2022.

Brun, E., Martin, E., Simon, V., Gendre, C., and Coléou, C.: An energy and mass model of snow cover suitable for operational avalanche forecasting, J. Glaciol., 35, 333–342, 1989.

Caponi, L., Formenti, P., Massabó, D., Di Biagio, C., Cazaunau, M., Pangui, E., Chevaillier, S., Landrot, G., Andreae, M. O., Kandler, K., Piketh, S., Saeed, T., Seibert, D., Williams, E., Balkanski, Y., Prati, P., and Doussin, J.-F.: Spectral- and size-resolved mass absorption efficiency of mineral dust aerosols in the shortwave spectrum: a simulation chamber study, Atmos. Chem. Phys., 17, 7175–7191, https://doi.org/10.5194/acp-17-7175-2017, 2017.

Carmagnola, C. M., Domine, F., Dumont, M., Wright, P., Strellis, B., Bergin, M., Dibb, J., Picard, G., Libois, Q., Arnaud, L., and Morin, S.: Snow spectral albedo at Summit, Greenland: measurements and numerical simulations based on physical and chemical properties of the snowpack, The Cryosphere, 7, 1139–1160, https://doi.org/10.5194/tc-7-1139-2013, 2013.

Chandrasekhar, S.: Radiative transfer, New York, Dover, 416 pp., https://openlibrary.org/books/OL28022891M/Radiative_Transfer (last access: 17 December 2024), 1960.

Chevrollier, L.-A., Cook, J. M., Halbach, L., Jakobsen, H., Benning, L. G., Anesio, A. M., and Tranter, M.: Light absorption and albedo reduction by pigmented microalgae on snow and ice, J. Glaciol., 69, 333–341, https://doi.org/10.1017/jog.2022.64, 2022.

Colbeck, S. C.: An Overview of Seasonal Snow Metamorphism, Rev. Geophys., 20, 45–61, https://doi.org/10.1029/RG020i001p00045, 1982.
Cook, J. M., Hodson, A. J., Gardner, A. S., Flanner, M., Tedstone, A. J., Williamson, C., Irvine-Fynn, T. D. L., Nilsson, J., Bryant, R., and Tranter, M.: Quantifying bioalbedo: a new physically based model and discussion of empirical methods for characteris- ing biological influence on ice and snow albedo, The Cryosphere, 11, 2611–2632, https://doi.org/10.5194/tc-11-2611-2017, 2017.

Dang, C., Zender, C. S., and Flanner, M. G.: Intercompari-
son and improvement of two-stream shortwave radiative trans-
fer schemes in Earth system models for a unified treat-
ment of cryospheric surfaces, The Cryosphere, 13, 2325–2343,
https://doi.org/10.5194/tc-13-2325-2019, 2019.

Di Mauro, B., Garzonio, R., Ravasio, C., Orlandi, V., Baccolo, G., Gilardoni, S., Remias, D., Leoni, B., Rossini, M., and Colombo, R.: Combined effect of algae and dust on snow spectral and broadband albedo, J. Quant. Spectrosc. Ra., 316, 108906, https://doi.org/10.1016/j.jqsrt.2024.108906, 2024.

Dombrovsky, L. A., Kokhanovsky, A. A., and Randrianalisoa, J. H.: On snowpack heating by solar radiation: A computational model, J. Quant. Spectrosc. Ra., 227, 72–85, https://doi.org/10.1016/j.jqsrt.2019.02.004, 2019.

Domine, F., Salvatori, R., Legagneux, L., Salzano, R., Fily, M., and Casacchia, R.: Correlation between the specific surface area and the short wave infrared (SWIR) reflectance of snow, Cold Reg. Sci. Technol., 46, 60–68, https://doi.org/10.1016/j.coldregions.2006.06.002, 2006.

Domine, F., Albert, M., Huthwelker, T., Jacobi, H.-W., Kokhanovsky, A. A., Lehning, M., Picard, G., and Simpson, W. R.: Snow physics as relevant to snow photochemistry, Atmos. Chem. Phys., 8, 171–208, https://doi.org/10.5194/acp-8-171-2008, 2008.

Donahue, C., Skiles, S. M., and Hammonds, K.: Mapping liquid water content in snow at the millimeter scale: an intercompari- son of mixed-phase optical property models using hyperspectral imaging and in situ measurements, The Cryosphere, 16, 43–59, https://doi.org/10.5194/tc-16-43-2022, 2022.

Dumont, M., Arnaud, L., Picard, G., Libois, Q., Lejeune, Y., Nabat, P., Voisin, D., and Morin, S.: In situ continuous visible and near- infrared spectroscopy of an alpine snowpack, The Cryosphere, 11, 1091–1110, https://doi.org/10.5194/tc-11-1091-2017, 2017.

Dumont, M., Tuzet, F., Gascoin, S., Picard, G., Kutuzov, S., Lafaysse, M., Cluzet, B., Nheili, R., and Painter, T. H.: Accelerated Snow Melt in the Russian Caucasus Mountains After the Saharan Dust Outbreak in March 2018, J. Geophys. Res.-Earth, 125, e2020JF005641, https://doi.org/10.1029/2020jf005641, 2020.

Dumont, M., Flin, F., Malinka, A., Brissaud, O., Hagenmuller, P., Lapalus, P., Lesaffre, B., Dufour, A., Calonne, N., Rolland du Roscoat, S., and Ando, E.: Experimental and model-based investigation of the links between snow bidirectional reflectance and snow microstructure, The Cryosphere, 15, 3921–3948, https://doi.org/10.5194/tc-15-3921-2021, 2021.

Dunkle, R. V. and Bevans, J. T.: An Approximate Analy-
sis of the Solar Reflectance and Transmittance of a Snow
Cover, J. Meteorol., 13, 212–216, https://doi.org/10.1175/1520-
0469(1956)013<0212:aaaots>2.0.co;2, 1956.

Emde, C., Buras-Schnell, R., Kylling, A., Mayer, B., Gasteiger, J., Hamann, U., Kylling, J., Richter, B., Pause, C., Dowling, T., and Bugliaro, L.: The libRadtran software package for radia-

tive transfer calculations (version 2.0.1), Geosci. Model Dev., 9,
1647–1672, https://doi.org/10.5194/gmd-9-1647-2016, 2016.

Fair, Z., Flanner, M., Schneider, A., and Skiles, S. M.: Sensitivity of modeled snow grain size retrievals to solar geometry, snow particle asphericity, and snowpack impurities, The Cryosphere, 16, 3801–3814, https://doi.org/10.5194/tc-16-3801-2022, 2022.

Fierz, C., Armstrong, R. L., Durand, Y., Etchevers, P., Greene, E., McClung, D. M., Nishimura, K., Satyawali, P. K., and Sokratoğlu, S. A.: The international classification for seasonal snow on the ground, UNESCO/IHP, https://unesdoc.unesco.org/ark:/48223/pf0000186462 (last access: 17 December 2024), 2009.

Flanner, M. G. and Zender, C. S.: Snowpack radiative heating: Influence on Tibetan Plateau climate, Geophys. Res. Lett., 32, L06501, https://doi.org/10.1029/2004GL022076, 2005.

Flanner, M. G., Shell, K. M., Barlage, M., Perovich, D. K., and Tschudi, M. A.: Radiative forcing and albedo feedback from the Northern Hemisphere cryosphere between 1979 and 2008, Nat. Geosci., 4, 151–155, https://doi.org/10.1038/ngeo1062, 2011.

Flanner, M. G., Liu, X., Zhou, C., Penner, J. E., and Jiao, C.: Enhanced solar energy absorption by internally-mixed black carbon in snow grains, Atmos. Chem. Phys., 12, 4699–4721, https://doi.org/10.5194/acp-12-4699-2012, 2012.

Flanner, M. G., Arnheim, J. B., Cook, J. M., Dang, C., He, C., Huang, X., Singh, D., Skiles, S. M., Whicker, C. A., and Zender, C. S.: SNICAR-ADv3: a community tool for modeling spectral snow albedo, Geosci. Model Dev., 14, 7673–7704, https://doi.org/10.5194/gmd-14-7673-2021, 2021.

France, J. L., King, M. D., Frey, M. M., Erbland, J., Picard, G., Preunkert, S., MacArthur, A., and Savarino, J.: Snow optical properties at Dome C (Concordia), Antarctica: implications for snow emissions and snow chemistry of reactive nitrogen, Atmos. Chem. Phys., 11, 9787–9801, https://doi.org/10.5194/acp-11-9787-2011, 2011.

Gallet, J.-C., Domine, F., Zender, C. S., and Picard, G.: Measure-
ment of the specific surface area of snow using infrared re-
flectance in an integrating sphere at 1310 and 1550 nm, The
Cryosphere, 3, 167–182, https://doi.org/10.5194/tc-3-167-2009,
2009.

Gallet, J.-C., Domine, F., Arnaud, L., Picard, G., and Savarino, J.: 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, The Cryosphere, 5, 631–649, https://doi.org/10.5194/tc-5-631-2011, 2011.

Glendinning, J. H. G. and Morris, E. M.: Incorporation of spectral and directional radiative transfer in a snow model, Hydrol. Process., 13, 1761–1772, https://doi.org/10.1002/(sici)1099-1085(199909)13:12/13<1761::aid-hyp856>3.0.co;2-y, 1999.

Green, R. O., Dozier, J., Roberts, D., and Painter, T.: Spectral snow-reflectance models for grain-size and liquid-water fraction in melting snow for the solar-reflected spectrum, Ann. Glaciol., 34, 71–73, https://doi.org/10.3189/17275640278187987, 2002.

Grenfell, T. C. and Warren, S. G.: Representation of a nonspherical ice particle by a collection of independent spheres for scattering and absorption of radiation, J. Geophys. Res., 104, 31697–31710, https://doi.org/10.1029/1999JD900496, 1999.

Gueymard, C. A.: Parameterized transmittance model for direct beam and circumsolar spectral irradiance, Sol. Energy, 71, 325–346, https://doi.org/10.1016/s0038-092x(01)00054-8, 2001.
Hall, A.: The Role of Surface Albedo Feedback in Climate, J. Climate, 17, 1550–1568, https://doi.org/10.1175/1520-0442(2004)017<1550:trosaf>2.0.co;2, 2004.

He, C. and Flanner, M.: Snow Albedo and Radiative Transfer: Theory, Modeling, and Parameterization, pp. 67–133, Springer International Publishing, ISBN 9783030386962, https://doi.org/10.1007/978-3-030-38696-2_3, 2020.

Henyey, L. C. and Greenstein, J. L.: Diffuse radiation in the Galaxy, Astrophys. J., 93, 70–83, https://doi.org/10.1086/144246, 1941.

Hoffer, A., Gelencsér, A., Guyon, P., Kiss, G., Schmid, O., Frank, G. P., Artaxo, P., and Andreae, M. O.: Optical properties of humic-like substances (HULIS) in biomass-burning aerosols, Atmos. Chem. Phys., 6, 3563–3570, https://doi.org/10.5194/acp-6-3563-2006, 2006.

Joseph, J. H., Wiscombe, W. J., and Weinman, J. A.: The Delta-Eddington Approximation for Radiative Flux Transfer, J. Atmos. Sci., 33, 2452–2459, https://doi.org/10.1175/1520-0469(1976)033<2452:tdeavr>2.0.co;2, 1976.

Kaempfer, T. U., Schneebeli, M., and Sokratov, S. A.: A microstructural approach to model heat transfer in snow, Geophys. Res. Lett., 32, L21503, https://doi.org/10.1029/2005GL023873, 2005.

Kang, S., Zhang, Y., Qian, Y., and Wang, H.: A review of black carbon in snow and ice and its impact on the cryosphere, Earth-Sci. Rev., 210, 103346, https://doi.org/10.1016/j.earscirev.2020.103346, 2020.

King, M. D. and Simpson, W. R.: Extinction of UV radiation in Arctic snow at Alert, Canada (82°N), J. Geophys. Res., 106, 12499–12507, https://doi.org/10.1029/2001jd900006, 2001.

Kokhanovsky, A. A.: Light scattering media optics, Springer Praxis books in environmental sciences, 3. edn., Praxis Publ., Chichester, UK, ISBN 978-3-54021184-6, https://link.springer.com/book/9783540211846 (last access: 17 December 2024), 2004.

Kokhanovsky, A. A.: Light Scattering Reviews, Vol. 6: Light Scattering and Remote Sensing of Atmosphere and Surface, Springer Berlin Heidelberg, ISBN 9783642155314, https://doi.org/10.1007/978-3-642-15531-4, 2012.

Kokhanovsky, A. A.: Light penetration in snow layers, J. Quant. Spectrosc. Ra., 278, 108040, https://doi.org/10.1016/j.jqsrt.2021.108040, 2022.

Kokhanovsky, A. A. and Macke, A.: Integral light-scattering and absorption characteristics of large, nonspherical particles, Appl. Optics, 36, 8785, https://doi.org/10.1364/ao.36.008785, 1997.

Kokhanovsky, A. A. and Zege, E. P.: Scattering optics of snow, Appl. Optics, 43, 1589–1602, 2004.

Krol, Q. and Löwe, H.: Relating optical and microwave grain metrics of snow: the relevance of grain shape, The Cryosphere, 10, 2847–2863, https://doi.org/10.5194/tc-10-2847-2016, 2016.

Larue, F., Picard, G., Arnaud, L., Ollivier, I., Delcourt, C., Lamare, M., Tuzet, F., Revuelto, J., and Dumont, M.: Snow albedo sensitivity to macroscopic surface roughness using a new ray-tracing model, The Cryosphere, 14, 1651–1672, https://doi.org/10.5194/tc-14-1651-2020, 2020.

Lee-Taylor, J. and Madronich, S.: Calculation of ac tinic fluxes with a coupled atmosphere–snow radiative transfer model, J. Geophys. Res.-Atmos., 107, 4796, https://doi.org/10.1029/2002jd002084, 2002.

Leroux, C. and Fily, M.: Modeling the effect of sastrugi on snow reflectance, J. Geophys. Res., 103, 25779, https://doi.org/10.1029/98je00558, 1998.

Letcher, T., Parno, J., Courville, Z., Farnsworth, L., and Olivier, J.: A generalized photon-tracking approach to simulate spectral snow albedo and transmittance using X-ray microtomography and geometric optics, The Cryosphere, 16, 4343–4361, https://doi.org/10.5194/tc-16-4343-2022, 2022.

Libois, Q., Picard, G., France, J. L., Arnaud, L., Dumont, M., Carmagnola, C. M., and King, M. D.: Influence of grain shape on light penetration in snow, The Cryosphere, 7, 1803–1818, https://doi.org/10.5194/tc-7-1803-2013, 2013.

Libois, Q., Picard, G., Arnaud, L., Dumont, M., Lafaysse, M., Morin, S., and Lefebvre, E.: Summertime evolution of snow specific surface area close to the surface on the Antarctic Plateau, The Cryosphere, 9, 2383–2398, https://doi.org/10.5194/tc-9-2383-2015, 2015.

Malinka, A.: Stereological approach to radiative transfer in porous materials. Application to the optics of snow, J. Quant. Spectrosc. Ra., 295, 108410, https://doi.org/10.1016/j.jqsrt.2022.108410, 2023.

Malinka, A. V.: Light scattering in porous materials: Geometrical optics and stereological approach, J. Quant. Spectrosc. Ra., 141, 14–23, https://doi.org/10.1016/j.jqsrt.2014.02.022, 2014.

Manninen, T., Anttila, K., Jääskeläinen, E., Riihelä, A., Peltoniemi, J., Räisänen, P., Lahtinen, P., Siljamo, N., Thölix, L., Meinan-der, O., Kontu, A., Suokanerva, H., Pirazzini, R., Suomalainen, J., Hakala, T., Kaasalainen, S., Kaartinen, H., Kukko, A., Haute-coeur, O., and Roujean, J.-L.: Effect of small-scale snow surface roughness on snow albedo and reflectance, The Cryosphere, 15, 793–820, https://doi.org/10.5194/tc-15-793-2021, 2021.

Masson, V., Le Moigne, P., Martin, E., Faroux, S., Alias, A., Alkama, R., Belamari, S., Barbu, A., Boone, A., Bouyssel, F., Brousseau, P., Brun, E., Calvet, J.-C., Carrer, D., Decharme, B., Delire, C., Donier, S., Essaouini, K., Gibelin, A.-L., Giordani, H., Habets, F., Jidane, M., Kerdraon, G., Kourzeneva, E., Lafaysse, M., Lafont, S., Lebeaupin Brossier, C., Lemonsu, A., Mahfouf, J.-F., Marguinaud, P., Mokhtari, M., Morin, S., Pigeon, G., Salgado, R., Seity, Y., Taillefer, F., Tanguy, G., Tulet, P., Vincendon, B., Vionnet, V., and Voldoire, A.: The SURFEXv7.2 land and ocean surface platform for coupled or offline simulation of earth surface variables and fluxes, Geosci. Model Dev., 6, 929–960, https://doi.org/10.5194/gmd-6-929-2013, 2013.

Meador, W. E. and Weaver, W. R.: Two-Stream Approximations to Radiative Transfer in Planetary Atmospheres: A Unified Description of Existing Methods and a New Improvement, J. Atmos. Sci., 37, 630–643, https://doi.org/10.1175/1520-0469(1980)037<0630:tsatr>2.0.co;2, 1980.

Mie, G.: Beitraege zur Optik trueber Medien, speziell kolloidaler Metalloesungen, Annals of Physics, 330, 377–445, 1908.

Mishchenko, M. I., Travis, L. D., and Mackowski, D. W.: T-matrix computations of light scattering by nonspherical particles: A review, J. Quant. Spectrosc. Ra., 55, 535–575, https://doi.org/10.1016/0022-4073(96)00002-7, 1996.

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

Onuma, Y., Takeuchi, N., Tanaka, S., Nagatsuka, N., Niwano, M., and Aoki, T.: Physically based model of the contribu-
tion of red snow algal cells to temporal changes in albedo
in northwest Greenland, The Cryosphere, 14, 2087–2101,
https://doi.org/10.5194/tc-14-2087-2020, 2020.

Perovich, D. K.: Light reflection and transmission by a temperate snow cover, J. Glaciol., 53, 201–210, https://doi.org/10.3189/172756507782202919, 2007.

Picard, G.: snowtartes: a web application to perform snow albedo and irradiance profile simulation with the TARTES model, https://snow.univ-grenoble-alpes.fr/snowtartes/ (last access: 17 December 2024), 2021.

Picard, G. and Libois Q.: ghislainp/tartes: v2.0, Zenodo [code], https://doi.org/10.5281/zenodo.13950598, 2024.

Picard, G., Arnaud, L., Domine, F., and Fily, M.: Determining snow specific surface area from near-infrared reflectance measurements: Numerical study of the influence of grain shape, Cold Reg. Sci. Technol., 56, 10–17, https://doi.org/10.1016/j.coldregions.2008.10.001, 2009.

Picard, G., Domine, F., Krinner, G., Arnaud, L., and Lefebvre, E.: Inhibition of the positive snow-albedo feedback by precipitation in interior Antarctica, Nat. Clim. Change, 2, 795–798, https://doi.org/10.1038/nclimate1590, 2012.

Picard, G., Libois, Q., and Arnaud, L.: Refinement of the ice absorption spectrum in the visible using radiance profile measurements in Antarctic snow, The Cryosphere, 10, 2655–2672, https://doi.org/10.5194/tc-10-2655-2016, 2016.

Picard, G., Arnaud, L., Caneill, R., Lefebvre, E., and Lamare, M.: Observation of the process of snow accumulation on the Antarctic Plateau by time lapse laser scanning, The Cryosphere, 13, 1983–1999, https://doi.org/10.5194/tc-13-1983-2019, 2019.

Qu, X. and Hall, A.: What Controls the Strength of Snow-Albedo Feedback?, J. Climate, 20, 3971–3981, https://doi.org/10.1175/jcli4186.1, 2007.

Räisänen, P., Makkonen, R., Kirkevåg, A., and Debernard, J. B.: Effects of snow grain shape on climate simulations: sensitivity tests with the Norwegian Earth System Model, The Cryosphere, 11, 2919–2942, https://doi.org/10.5194/tc-11-2919-2017, 2017.

Réveillet, M., Dumont, M., Gascoin, S., Lafaysse, M., Nabat, P., Ribes, A., Nheili, R., Tuzet, F., Ménégoz, M., Morin, S., Picard, G., and Ginoux, P.: Black carbon and dust alter the response of mountain snow cover under climate change, Nat. Commun., 13, 5279, https://doi.org/10.1038/s41467-022-32501-y, 2022.

Ricchiazzi, P., Yang, S., Gautier, C., and Sowle, D.: SB-DART: A Research and Teaching Software Tool for Plane-Parallel Radiative Transfer in the Earth’s Atmosphere, B. Am. Meteorol. Soc., 79, 2101–2114, https://doi.org/10.1175/1520-0477(1998)079<2101:SARATS>2.0.CO;2, 1998.

Robledano, A., Picard, G., Arnaud, L., Larue, F., and Ollivier, I.: Modelling surface temperature and radiation budget of snow-covered complex terrain, The Cryosphere, 16, 559–579, https://doi.org/10.5194/tc-16-559-2022, 2022.

Robledano, A., Picard, G., Dumont, M., Flin, F., Arnaud, L., and Libois, Q.: Unraveling the optical shape of snow, Nat. Commun., 14, 3955, https://doi.org/10.1038/s41467-023-39671-3, 2023.

Schaepman-Strub, G., Schaepman, M., Painter, T., Dangel, S., and Martonchik, J.: Reflectance quantities in optical remote sensing – definitions and case studies, Remote Sens. Environ., 103, 27–42, https://doi.org/10.1016/j.rse.2006.03.002, 2006.

Shao, D., Xu, W., Li, H., Wang, J., and Hao, X.: Reconstruction of Remotely Sensed Snow Albedo for Quality Im-

provements Based on a Combination of Forward and Retrieval Models, IEEE T. Geosci. Remote, 56, 6969–6985, https://doi.org/10.1109/tgrs.2018.2846681, 2018.

Stamnes, K., Tsay, S. C., Jayaweera, K., and Wiscombe, W.: Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media, Appl. Optic., 27, 2502–2509, https://doi.org/10.1364/AO.27.002502, 1988a.

Stamnes, K., Tsay, S.-C., and Nakajima, T.: Computation of eigenvalues and eigenvectors for the discrete ordinate and matrix operator methods in radiative transfer, J. Quant. Spectrosc. Ra., 39, 415–419, https://doi.org/10.1016/0022-4073(88)90107-0, 1988b.

Stephens, G. L.: Radiation Profiles in Extended Water Clouds. II: Parameterization Schemes, J. Atmos. Sci., 35, 2123–2132, https://doi.org/10.1175/1520-0469(1978)035<2123:rpiewc>2.0.CO;2, 1978.

Toon, O. B., McKay, C. P., Ackerman, T. P., and Santhanam, K.: Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres, J. Geophys. Res.-Atmos., 94, 16287–16301, https://doi.org/10.1029/jd094id13p16287, 1989.

Tuzet, F., Dumont, M., Lafaysse, M., Picard, G., Arnaud, L., Voisin, D., Lejeune, Y., Charrois, L., Nabat, P., and Morin, S.: A multilayer physically based snowpack model simulating direct and indirect radiative impacts of light-absorbing impurities in snow, The Cryosphere, 11, 2633–2653, https://doi.org/10.5194/tc-11-2633-2017, 2017.

Tuzet, F., Dumont, M., Arnaud, L., Voisin, D., Lamare, M., Larue, F., Revuelto, J., and Picard, G.: Influence of light-absorbing particles on snow spectral irradiance profiles, The Cryosphere, 13, 2169–2187, https://doi.org/10.5194/tc-13-2169-2019, 2019.

Tuzet, F., Dumont, M., Picard, G., Lamare, M., Voisin, D., Nabat, P., Lafaysse, M., Larue, F., Revuelto, J., and Arnaud, L.: Quantification of the radiative impact of light-absorbing particles during two contrasted snow seasons at Col du Lautaret (2058 m a.s.l., French Alps), The Cryosphere, 14, 4553–4579, https://doi.org/10.5194/tc-14-4553-2020, 2020.

Usha, K. H., Nair, V. S., and Babu, S. S.: Modeling of aerosol induced snow albedo feedbacks over the Himalayas and its implications on regional climate, Clim. Dynam., 54, 4191–4210, https://doi.org/10.1007/s00382-020-05222-5, 2020.

van Dalum, C. T., van de Berg, W. J., Libois, Q., Picard, G., and van den Broeke, M. R.: A module to convert spectral to narrowband snow albedo for use in climate models: SNOWBAL v1.2, Geosci. Model Dev., 12, 5157–5175, https://doi.org/10.5194/gmd-12-5157-2019, 2019.

Veillon, F., Dumont, M., Amory, C., and Fructus, M.: A versatile method for computing optimized snow albedo from spectrally fixed radiative variables: VALHALLA v1.0, Geosci. Model Dev., 14, 7329–7343, https://doi.org/10.5194/gmd-14-7329-2021, 2021.

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snow-pack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791, https://doi.org/10.5194/gmd-5-773-2012, 2012.

Walden, V. P., Warren, S. G., and Tuttle, E.: Atmospheric Ice Crystals over the Antarctic Plateau in Winter, J. Appl.
Meteorol., 42, 1391-1405, https://doi.org/10.1175/1520-0450(2003)042<1391:AICOTA>2.0.CO;2, 2003.

Warren, S. G. and Brandt, R. E.: Optical constants of ice from the ultraviolet to the microwave: A revised compilation, J. Geophys. Res., 113, D14220, https://doi.org/10.1029/2007JD009744, 2008.

Warren, S. G. and Wiscombe, W. J.: A Model for the Spectral Albedo of Snow. II: Snow Containing Atmospheric Aerosols, J. Atmos. Sci., 37, 2734-2745, https://doi.org/10.1175/1520-0469(1980)037<2734:amftsa>2.0.co;2, 1980.

Warren, S. G., Brandt, R. E., and O'Rawe Hinton, P.: Effect of surface roughness on bidirectional reflectance of Antarctic snow, J. Geophys. Res., 103, 25789, https://doi.org/10.1029/98je01898, 1998.

Wiscombe, W. J.: The delta-Eddington approximation for a vertically inhomogeneous atmosphere, National Center for Atmospheric Research Boulder, Colorado, https://doi.org/10.5065/D65H7D6Z, 1977.

Wiscombe, W. J. and Warren, S. G.: A Model for the Spectral Albedo of Snow. I: Pure Snow, J. Atmos. Sci., 37, 2712-2733, https://doi.org/10.1175/1520-0469(1980)037<2712:AMFTSA>2.0.CO;2, 1980.

Xiong, C. and Shi, J.: Simulating polarized light scattering in terrestrial snow based on bicontinuous random medium and Monte Carlo ray tracing, J. Quant. Spectrosc. Ra., 133, 177-189, https://doi.org/10.1016/j.jqsrt.2013.07.026, 2014.