# Fast, and Accurate Radiative Transfer for Land Surface Models

Kazem Ardaneh, Fabienne Maignan, Sebastiaan Luyssaert, Philippe Peylin, Olivier Boucher

► To cite this version:

Kazem Ardaneh, Fabienne Maignan, Sebastiaan Luyssaert, Philippe Peylin, Olivier Boucher. Fast, and Accurate Radiative Transfer for Land Surface Models. Procedia Computer Science, 2025, 255, pp.3 - 12. 10.1016/j.procs.2025.02.255. hal-05019422

HAL Id: hal-05019422

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

Submitted on 3 Apr 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Proceedings of the Second EuroHPC user day

Fast, and Accurate Radiative Transfer for Land Surface Models

Kazem Ardaneh<sup>a,*</sup>, Fabienne Maignan<sup>b</sup>, Sebastiaan Luyssaert<sup>c</sup>, Philippe Peylin<sup>b</sup>, Olivier Boucher<sup>a</sup>

Abstract

Land surface models (LSMs) simulate processes occurring at the Earth's surface (including those related to vegetation, soil, and hydrology) and their interactions with the atmosphere. LSMs are crucial for environmental monitoring, weather forecasting, and climate studies. The radiative transfer (RT) through vegetation canopies is an important process that determines photosynthesis, and the land surface energy budget. Conventional multilayer iterative solvers for RT through the canopy are computationally demanding. Here, we develop a multilayer matrix-based RT solver for vegetation canopies. The results show that the solver matches the accuracy of existing models and significantly reduces the computational time for RT, highlighting its potential for practical applications.

© 2025 The Authors. Published by Elsevier B.V.

This is an open access article under the CC BY 4.0 license (https://creativecommons.org/licenses/by/4.0)

Peer-review under responsibility of the scientific committee of the Proceedings of the Second EuroHPC user day

Keywords: Land surface modeling; ORCHIDEE; Vegetation canopies; Radiative transfer; Two-stream model; Matrix-based solver

# 1. Introduction

Climate models help to understand and predict the Earth's future climate by simulating the interactions between the atmosphere, oceans, and land surfaces. Within these models, land surface models (LSMs) represent the exchange of carbon, water, and energy, between the land and the atmosphere. ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems) [1] is a process-based LSM developed at the Institut Pierre-Simon Laplace (IPSL). This model is widely used in Earth System Modeling and integrates experimental, inventory, and monitoring data from field sites, atmospheric stations, and satellites to refine predictions of land-atmosphere exchange and feedback in different climates and ecosystems.

Radiative transfer (RT) through the vegetation canopy is important when modeling the land surface as it influences photosynthesis, the energy balance of the land surface, and the climate system. The vegetation canopy interacts with

<footnote><sup>a</sup>Institut Pierre-Simon Laplace (IPSL), Sorbonne Université / CNRS, Paris, France</footnote>

<footnote><sup>b</sup>Laboratoire des sciences du climat et de l'environnement (LSCE), IPSL, CEA/CNRS/UVSQ, Gif-sur-Yvette, France</footnote>

<footnote><sup>c</sup>Section Systems Ecology, Amsterdam Institute for Life and Environment (A-LIFE), Vrije Universiteit Amsterdam, the Netherlands</footnote>

<footnote>* Corresponding author</footnote>

<footnote>E-mail address: kardaneh@ipsl.fr</footnote>
incoming solar radiation through absorption, reflection, and transmission by the vegetation elements consisting of leaves, branches, and stems and their background consisting of litter and soils. A two-stream approximation in the RT model is widely used to simulate radiation's interaction with vegetation canopy. It simplifies the RT equation by considering only two radiation directions: direct (downward) and diffuse (upward and downward). In spite of the simplification, such an approach captures the basic characteristics of radiation interactions within the vegetation canopy while remaining simple.

In considering existing approaches to RT within vegetation canopies, Ref. [2] provided analytical solutions of the two-stream equations for vegetation canopies as presented in Ref. [3], which is used in the Community Land Model (CLM) as detailed in Ref. [4]. The Sellers model, although pioneering, is limited by its single-layer structure. Additionally, Ref. [5] presented a single-layer two-stream model incorporating the black background equations of Ref. [6] and considering successive reflections from the surface. The multilayer version of Ref. [5], presented in Ref. [7], addresses the limits of the single-layer approach by considering the vertical heterogeneity of the vegetation canopies. However, it relies on an iterative approach to find the correct reflection, transmission, and absorption for each layer, making it computationally demanding as it currently accounts for almost 30% of the computing time of the ORCHIDEE LSM.

In this study, we developed a multilayer matrix-based two-stream RT solver for vegetation canopy based on Ref. [8]. We formulate the RT model as a system of linear equations, which can then be solved using efficient matrix operations. This allows us to harness the power of linear algebra algorithms and exploit parallelism more effectively than traditional iterative methods. The solver is validated against the models of Refs. [2, 5, 7]. While producing almost the same results, the new solver considerably reduces the computational cost for RT in the ORCHIDEE code and argues in favor of adopting matrix methods in practical applications.

## 2. Matrix-based RT

The key to our approach is the matrix-based RT, which deals with the interactions of radiative fluxes within and between multiple layers. This approach is particularly well suited for complex vegetation canopies for which their vertical heterogeneity can have a significant impact on RT. In the two-stream approximation, the basic equations for the upward ($F^\uparrow$) and downward ($F^\downarrow$) diffuse fluxes are given by [8]:

$$ \begin{cases} \frac{dF^\uparrow}{d\tau} = \gamma_1 F^\uparrow - \gamma_2 F^\downarrow - \gamma_3 \omega \pi F_\odot \exp(-\tau/\mu_0) \\ \\ \frac{dF^\downarrow}{d\tau} = \gamma_2 F^\uparrow - \gamma_1 F^\downarrow + \gamma_4 \omega \pi F_\odot \exp(-\tau/\mu_0) \end{cases} \qquad (1) $$

The gamma coefficients ($\gamma_1, \gamma_2, \gamma_3$, and $\gamma_4 = 1 - \gamma_3$) are key parameters in the two-stream approximation that prescribe the scattering and absorption properties of the medium and are thus essential in determining the radiative fluxes within each layer. In particular, $\gamma_1$, and $\gamma_2$ relate to the interactions of upward and downward fluxes, while $\gamma_3$ and $\gamma_4$ represent the source terms associated with scattering and absorption processes. The values of these coefficients can be derived from various approximation methods such as the Eddington approximation, quadrature, and $\delta$-methods [6]. In Eqs. (1), $\tau$ is the optical depth, $\omega$ is the single-scattering albedo, $\pi F_\odot$ is the incident flux at the top of the canopy, and $\mu_0$ is the cosine of the solar zenith angle.

When in Eqs. (1) $\gamma_1$ is substituted with $[1 - \omega(1 - \beta)]/\bar{\mu}$, $\gamma_2$ is replaced by $\omega\beta/\bar{\mu}$, $\gamma_3$ is taken as $\beta_0$, $\mu_0$ is redefined as $1/K$, $\pi F_\odot$ is set equal to $K$, and $\tau$ is replaced with $L$ (leaf area index) the two-stream RT equations for the vegetation canopies are obtained [2, 3, 5]:

$$ \begin{cases} \frac{dF^\uparrow}{dL} = \frac{[1 - \omega(1 - \beta)]}{\bar{\mu}} F^\uparrow - \frac{\omega\beta}{\bar{\mu}} F^\downarrow - \beta_0 \omega K \exp(-KL) \\ \\ \frac{dF^\downarrow}{dL} = \frac{\omega\beta}{\bar{\mu}} F^\uparrow - \frac{[1 - \omega(1 - \beta)]}{\bar{\mu}} F^\downarrow + (1 - \beta_0) \omega K \exp(-KL) \end{cases} \qquad (2) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. RT through a multilayer medium. The incoming solar flux $F_\odot$ is incident at an angle $\theta_0$ with the surface normal. The medium is split into N layers, with optical depths of $\Delta\tau_l$. The total optical depth above the lth layer is $\tau_c$, while $\tau$ gives the optical depth within a layer measured from the top of that layer. For each layer l, the downward diffuse flux $F_l^\downarrow$ is partially scattered and absorbed, giving an upward diffuse flux $F_l^\uparrow$.</figcaption></figure>

Extending the single-layer two-stream model [Eqs. (1)] to multiple layers involves calculating the fluxes within each layer, considering the cumulative optical depths ($\tau_c$). In a given layer $l$, the fluxes are expressed as follows:

$$ 
\left\{
\begin{aligned}
F_l^\downarrow(\tau) &= \Gamma_l B_{l2} \exp(\lambda_l \tau) + B_{l1} \exp(-\lambda_l \tau) + C_l^\downarrow \exp[-(\tau + \tau_c)/\mu_0] \\
F_l^\uparrow(\tau) &= B_{l2} \exp(\lambda_l \tau) + \Gamma_l B_{l1} \exp(-\lambda_l \tau) + C_l^\uparrow \exp[-(\tau + \tau_c)/\mu_0] \\
\lambda_l &= (\gamma_{l1}^2 - \gamma_{l2}^2)^{1/2} \\
\Gamma_l &= \frac{\gamma_{l1} - \lambda_l}{\gamma_{l2}} \\
C_l^\uparrow &= \frac{\omega_0 \pi F_\odot [(\gamma_{l1} - 1/\mu_0)\gamma_{l3} + \gamma_{l4}\gamma_{l2}]}{(\lambda_l^2 - 1/\mu_0^2)} \\
C_l^\downarrow &= \frac{\omega_0 \pi F_\odot [(\gamma_{l1} + 1/\mu_0)\gamma_{l4} + \gamma_{l2}\gamma_{l3}]}{(\lambda_l^2 - 1/\mu_0^2)}
\end{aligned}
\right.
\qquad (3)
$$
Equations (3) gives the analytical solution for Eqs. (1), where $\tau$ is the optical depth within each layer, $\tau_c$ is the cumulative optical depth at the top of each layer (Fig. 1). The coefficients $B_{l1}$ and $B_{l2}$ are derived from the boundary conditions and the continuity at the interfaces between the layers. For numerical stability and to keep all arguments to the exponential terms negative, the variables $A_l$ and $B_l$ as in Eqs. (4) are introduced. These variables represent the modified flux amplitudes within each layer, derived such as to maintain consistent sign conventions between the exponential terms.

$$ \left\{ \begin{aligned} A_l &= \frac{B_{l1} \exp(\lambda_l \Delta \tau_l) + B_{l2}}{2} \\ B_l &= \frac{B_{l1} \exp(\lambda_l \Delta \tau_l) - B_{l2}}{2} \end{aligned} \right. \qquad (4) $$

In these equations, $\Delta\tau$ represents the total optical depth of the layer (Fig. 1). The downward and upward fluxes then read:

$$ \left\{ \begin{aligned} F_l^\uparrow(\tau) &= A_l \{\exp[-\lambda_l(\Delta\tau_l - \tau)] + \Gamma_l \exp(-\lambda_l\tau)\} + B_l \{\exp[-\lambda_l(\Delta\tau_l - \tau)] - \Gamma_l \exp(-\lambda_l\tau)\} + C_l^\uparrow \exp[-(\tau + \tau_c)/\mu_0] \\ F_l^\downarrow(\tau) &= A_l \{\Gamma_l \exp[-\lambda_l(\Delta\tau_l - \tau)] + \exp(-\lambda_l\tau)\} + B_l \{\Gamma_l \exp[-\lambda_l(\Delta\tau_l - \tau)] - \exp(-\lambda_l\tau)\} + C_l^\downarrow \exp[-(\tau + \tau_c)/\mu_0] \end{aligned} \right. \qquad (5) $$

In multilayer RT, flux continuity and boundary conditions are maintained using a matrix of coefficients. The continuity equations for upward and downward fluxes at the interfaces of layers $l$ and $l + 1$ are given by:

$$ \left\{ \begin{aligned} F_l^\uparrow(\tau = \Delta\tau_l) &= F_{l+1}^\uparrow(\tau = 0), & \text{for } l = 1, 2, \dots, (N-1) \\ F_l^\downarrow(\tau = \Delta\tau_l) &= F_{l+1}^\downarrow(\tau = 0), & \text{for } l = 1, 2, \dots, (N-1) \end{aligned} \right. \qquad (6) $$

Using Eqs. (5), Eqs. (6) can be expressed as a system of linear equations:

$$ \left\{ \begin{aligned} e_l A_l + d_l B_l - f_{l+1} A_{l+1} + g_{l+1} B_{l+1} &= C_{l+1}^{\uparrow}(0) - C_l^{\uparrow}(\Delta\tau_l) \\ f_l A_l + g_l B_l - e_{l+1} A_{l+1} + d_{l+1} B_{l+1} &= C_{l+1}^{\downarrow}(0) - C_l^{\downarrow}(\Delta\tau_l) \end{aligned} \right. \qquad (7) $$

where the coefficients $e_l$, $d_l$, $f_l$, and $g_l$ are defined as:

$$ \left\{ \begin{aligned} e_l &= 1 + \Gamma_l \exp(-\lambda_l \Delta \tau_l) \\ d_l &= 1 - \Gamma_l \exp(-\lambda_l \Delta \tau_l) \\ f_l &= \Gamma_l + \exp(-\lambda_l \Delta \tau_l) \\ g_l &= \Gamma_l - \exp(-\lambda_l \Delta \tau_l) \end{aligned} \right. \qquad (8) $$

We derive a tridiagonal system of equations to reduce computational cost. The following steps outline the derivation process. Multiply the first equation in Eqs. (7) by $d_{l+1}$ and the second by $g_{l+1}$, subtract the resulting equations to eliminate $B_{l+1}$, leading to a new equation that relates $A_l$, $B_l$, and $A_{l+1}$. Multiply the second equation in Eqs. (7) by $d_{l+1}e_l - g_{l+1}f_l$ and first by $f_l$, subtract these results to eliminate $A_{l+1}$, leading to another new equation that relates $A_l$, $B_l$, and $B_{l+1}$. The final system of equations reads:

$$ \left\{ \begin{aligned} & \underbrace{\alpha_{l1}}_{[d_{l+1} e_l - f_l g_{l+1}]} A_l + \underbrace{\beta_{l1}}_{[d_{l+1} d_l - g_l g_{l+1}]} B_l + \underbrace{\gamma_{l1}}_{[e_{l+1} g_{l+1} - f_{l+1} d_{l+1}]} A_{l+1} = \\ & \qquad \underbrace{[C_{l+1}^{\uparrow}(0) - C_l^{\uparrow}(\Delta\tau_l)] d_{l+1} + [C_l^{\downarrow}(\Delta\tau_l) - C_{l+1}^{\downarrow}(0)] g_{l+1}}_{\chi_{l1}} \\ & \underbrace{\alpha_{l2}}_{[d_l f_l - e_l g_l]} B_l + \underbrace{\beta_{l2}}_{[e_{l+1} e_l - f_l f_{l+1}]} A_{l+1} + \underbrace{\gamma_{l2}}_{[f_l g_{l+1} - e_l d_{l+1}]} B_{l+1} = \\ & \qquad \underbrace{[C_{l+1}^{\uparrow}(0) - C_l^{\uparrow}(\Delta\tau_l)] f_l + [C_l^{\downarrow}(\Delta\tau_l) - C_{l+1}^{\downarrow}(0)] e_l}_{\chi_{l2}} \end{aligned} \right. \qquad (9) $$
where the coefficients $\alpha_{l1}, \beta_{l1}, \gamma_{l1}$ and $\alpha_{l2}, \beta_{l2}, \gamma_{l2}$ are derived from the original coefficients $e_l, d_l, f_l, g_l$ and the source terms $C_l^\uparrow, C_l^\downarrow$.

The downward diffuse flux at the top of the multilayer structure is equal to the incident diffuse flux. In addition, the upward flux at the bottom of the multilayer structure is equal to the product of the downward flux and the reflectance of the surface ($R_S$), i.e. the reflected incident flux at the surface is all diffuse. Hence, boundary conditions at the top and bottom boundaries are as follows:

$$ 
\left\{ 
\begin{aligned}
F_1^\downarrow(\tau = 0) &= A_1\{\Gamma_1 \exp[-\lambda_1(\Delta\tau_1)] + 1\} + B_1\{\Gamma_1 \exp[-\lambda_1(\Delta\tau_1)] - 1\} + C_1^\downarrow(0) \\
&= \text{Downward diffuse flux}
\end{aligned}
\right. 
\qquad (10) \\
F_N^\uparrow(\tau = \Delta\tau_N) = A_N\{1 + \Gamma_N \exp(-\lambda_N\Delta\tau_N)\} + B_N\{1 - \Gamma_N \exp(-\lambda_N\Delta\tau_N)\} + C_N^\uparrow(\Delta\tau_N) \\
= R_S F_N^\downarrow(\tau = \Delta\tau_N) + R_S F_\odot(\tau = \Delta\tau_N)
$$

We need to solve a tridiagonal system ($2N - 2$ equations) to determine $A_l, B_l, A_{l+1}, B_{l+1}$ [Eqs. (9)] and the related variables such as $F_l^\downarrow$ and $F_l^\uparrow$. We use the standard method of tridiagonal solving, known as the Thomas algorithm. This algorithm is efficient for solving tridiagonal systems, as it reduces computational complexity compared with general matrix solvers.

Several key parameters need to be precalculated to model RT in vegetation. These are $G(\mu)$, $K(\mu)$, $\bar{\mu}$, $\omega\beta$, and $\omega\beta_0$. The function $G(\mu)$, known as the asymmetry factor of the phase function, is given by [9]:

$$ G(\mu) = \frac{1}{2\pi} \int_{\Omega'} d\hat{\Omega}' g(\theta')h(\phi') |\hat{\Omega} \cdot \hat{\Omega}'|, \qquad (11) $$

where $g(\theta)$ and $(1/2\pi)h(\phi)$ are the probability density functions of leaf normal inclination and azimuth, respectively. Several common functions exist for the probability distribution $g(\theta) \sin\theta$. For instance, the function is $\sin\theta$ in the spherical case while the uniform function is $2/\pi$. Once the $G$ function has been determined, the other parameters can be directly derived. These are:

$$ 
\left\{ 
\begin{aligned}
K(\mu) &= \frac{G(\mu)}{\mu} \\
\bar{\mu} &= \int_0^1 d\mu' \frac{\mu'}{G(\mu')} \\
\omega &= (r + t) \\
\delta &= (r - t) \\
\omega\beta &= \frac{1}{2} \left( \omega + \delta \int_0^{\pi/2} d\theta \sin\theta \cos^2\theta g(\theta) \right) \\
\omega\beta_0 &= \frac{1}{2} \left( \omega + \frac{\mu_0}{G(\mu_0)} \delta \int_0^{\pi/2} d\theta \sin\theta \cos^2\theta g(\theta) \right)
\end{aligned}
\right. 
\qquad (12) \\
\text{where } r \text{ is the reflectivity and } t \text{ is the transmissivity of vegetation (leaf). The last two equations in Eqs. (12) represent the back-scattering parameters } (\beta, \beta_0) \text{ for the diffuse and direct beams, respectively, as defined in Ref. [5]. However, these two parameters need to be adjusted for other cases. In Refs. [3, 2, 4], the back-scattering parameter for the diffuse beam is defined as}
$$

$$ \omega\beta = \frac{1}{2} (\omega + \delta \cos^2 \Theta) \qquad (13) $$

where $\Theta$ is the mean leaf inclination angle relative to the horizontal surface. The back-scattering parameter for direct beam, $\beta_0$ in Refs. [3, 2, 4] reads:

$$ \omega\beta_0 = \frac{1 + \bar{\mu} K}{\bar{\mu} K} \alpha_{ss}(\mu), \qquad (14) $$

where $\alpha_{ss}(\mu)$, the single scattering albedo of the canopy, is given by

$$ \alpha_{ss}(\mu) = \omega \int_0^1 d\mu' \frac{\mu' G(\mu)}{\mu G(\mu') + \mu' G(\mu)}. \qquad (15) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Comparison of visible albedo ($\alpha_{VIS}$) and near-infrared albedo ($\alpha_{NIR}$) as a function of the solar zenith angle for different leaf area index (LAI). The blue circles represent the matrix-based solver results and the red circles represent the results of the model by Ref. [2].</figcaption></figure>

### 3. Verification and validation

The matrix-based RT solver has been verified to ensure its accuracy and reliability. The verification process confirms that the solver produces the same results for single-layer and multilayer configurations, provided the same properties are applied to all layers. The solver was validated to reproduce the results presented in Ref. [8]. Specifically, the results of the solver were compared to the data in Tables (4), (5), and (6) of the referenced article and were found to be in exact agreement. We further evaluated the solver against Ref. [2], Ref. [5], and the multilayer iterative radiation solver in the ORCHIDEE model as presented in Ref. [7].

Fig. 2 shows the comparison between the matrix solver and the canopy reflectance model developed by Ref. [2] for different leaf area index (LAI) values. We note that an LAI value of 100 is excessive and was chosen to represent the infinite case in Ref. [2]. The top row shows the visible albedo ($\alpha_{VIS}$) while the bottom row shows the near-infrared albedo ($\alpha_{NIR}$) for solar zenith angles ranging from 0° to 90°. For each LAI value, the matrix solver results (blue circles) are plotted next to the Sellers model results (red circles). The matrix solver tends to follow the general trend of the Sellers model with a maximum relative difference of 8%.

Fig. 3 illustrates the predicted absorption of photosynthetically active radiation (PAR) for a growing wheat canopy, compared between the matrix-based solver and the Sellers model for different solar angles (10°, 60°, and 80°). Absorption percentages are plotted against LAI. At lower solar zenith angles, both models show high absorption rates, which saturate as the LAI increases. At higher solar zenith angles, the absorption decreases with increasing LAI. There is a maximum relative difference of 5% between the Sellers model (red circles) and the matrix solver (blue circles).
<figure><img src="image_4.png" /><figcaption>Fig. 3. Predicted absorption of photosynthetically active radiation (PAR) for a growing wheat canopy, compared between the matrix-based solver and Ref. [2] model at solar zenith angles of 10°, 60° and 80°. The blue circles represent the matrix-based solver results while the red circles represent the results of the Ref. [2] model.</figcaption></figure>

Fig. 4 shows a comparison of the directional hemispherical reflectance (DHR) calculated using the matrix-based solver against the results of Ref. [5]. The comparison is performed for different vegetation densities (sparse, medium, dense) and spectral bands (near-infrared, red, snow cover condition in the near-infrared), with each row corresponding to a different vegetation density and each column representing a different spectral band. DHR values are plotted against the solar zenith angle, ranging from 0° to 75°.

In the top row, for sparse vegetation, the near-infrared band shows increasing DHR with the solar zenith angle for both solvers although the matrix solver maximally gives 2% higher DHR at higher angles. For the red band, the DHR decreases with an increasing solar zenith angle, with a maximum difference of 3% between the two approaches. For snow cover, the DHR remains high and relatively constant, with the matrix solver results being at most 2% higher than those in Ref. [5], particularly at lower solar zenith angles.

The middle row, representing the average density of vegetation, shows a similar trend in the near-infrared band with an increase in DHR with a maximum difference of 2% between the two solvers. The red band also shows a decrease in DHR with increasing solar zenith angle, with both methods producing similar results, with a maximum difference of 10% at higher angles. For snow cover, the DHR values are lower than for sparse vegetation but remain relatively stable, with the matrix solver results being at most 5% higher at lower solar zenith angles.

In the bottom row, for dense vegetation, the near-infrared band shows a sharper increase in DHR with solar zenith angle compared to sparse and medium vegetation, with the matrix solver closely matched to Ref. [5], although with a maximum difference of 3% at higher angles. The red band shows an almost stable DHR with increasing solar zenith angle, and the matrix solver results have a maximum difference of 7% relative to those in Ref. [5]. For snow cover, the DHR values also show a maximum difference of 5% between the two solvers.

The matrix-based RT solver has been further compared with the ORCHIDEE iterative solver [7] across all land points, Plant Functional Types (PFTs), spectral bands, and for different simulation durations. For the meteorological forcing data, we used the CRUJRA dataset [10], which provides 6-hourly, gridded data for 1901–2018. This dataset integrates the observational data from the Climatic Research Unit (CRU) with the Japanese Reanalysis (JRA). It includes variables such as temperature, and precipitation at a 0.5° spatial resolution. The comparison was configured following a default ORCHIDEE set-up called FG_CRUJRA_SPIN. This setup uses 2° ×2° CRUJRA forcing, repeatedly cycling over 1901–1910 for a 320-year long-term spinup. It begins from scratch, with 15 PFTs, no land cover changes, fixed nitrogen input (1850, 1860, or 1900 depending on nitrogen species), and CO₂ concentrations fixed to the year 1860. Following spinup the simulation continues with a default transient set-up called FG_CRUJRA_TRANS. FG_CRUJRA_TRANS uses CRUJRA forcing data over 1861–1900, restarting from December 31st of the last year
<figure><img src="image_5.png" /><figcaption>Fig. 4. Comparison of the directional hemispheric reflectance (DHR) for three spectral bands (Near-Infrared, Red, and snow cover condition in the near-infrared) and three vegetation densities (sparse, medium, and dense) across varying solar zenith angles. Blue circles represent results from the matrix-based solver while red circles represent results from Ref. [5].</figcaption></figure>

of FG_CRUJRA_SPIN. It includes 15 PFTs, fixed land cover, nitrogen input for 1900, and annually varying CO₂ concentrations. The simulation is concluded with a historical default run. For ORCHIDEE this set-up is called FG_CRUJRA_HIST and applies the annual CRUJRA forcing from 1901 to 2010. The set-up starts from the last day of the FG_CRUJRA_TRANS, it includes annual updates for land cover, nitrogen input, and CO₂ concentrations, reflecting dynamic historical conditions.

The comparison is based on three main parameters: surface albedo, absorption, and transmission to the ground.
Fig. 5 shows α<sub>VIS</sub> [panels (a) and (b)] and α<sub>NIR</sub> [panels (c) and (d)], for the iterative (top row) and the matrix solver
<figure><img src="image_6.png" /><figcaption>Fig. 5. Visible albedo $\alpha_{VIS}$ [(a) and (b)] and near-infrared albedo $\alpha_{NIR}$ [(c) and (d)], for the iterative solver (top row) and the matrix solver (middle row) at the snapshot of 04/30/1902 for a two-year FG_CRUJRA_HIST simulation, the density plot (logarithm base 10) for the absolute and relative differences between the albedo of the matrix solver and the iterative one on all data points and the full two-year FG_CRUJRA_HIST simulation [(e)], comparison of the computation speed between the two solvers [(f)].</figcaption></figure>
(middle row) at the snapshot time of 30/04/1902 for a two-year FG_CRUJRA_HIST simulation. Both solvers give almost the same results. Fig. 5e shows the density plot of relative and absolute differences for comparisons between the albedos of the two solvers (all land points, all spectral bands, and all PFTs accumulated over two years). The scatterplot highlights the presence of a dense cluster of points near the origin, indicating that in many cases both the absolute and relative differences are small and the two solvers give very similar results. However, there are also distinct bands and peaks where the absolute or relative differences are more pronounced (representing approximately one percent of data points). The behavior of these differences varies depending on the magnitude of the values being compared. The relative difference is important when comparing two very small numbers, while the absolute difference is important when comparing two large numbers.

The tests presented in this study were run on a Leonardo Booster machine and compiled using the NVIDIA NVHPC compiler with the options -i4 -r8 -O2 -Kieee -Ktrap=fp, the code was run in a fully MPI setup, and on CPUs. The matrix method is approximately six times faster than the iterative (Fig. 5f). For performance comparison, the RT solver runtime is summed over two years. The matrix-based solver is faster because it calculates the interactions between layers using derived analytical equations, thus avoiding the repeated calculations required in iterative methods.

## 4. Conclusion

The presented results indicate that the matrix-based solver generally provides results in good agreement with other LSM RT solvers, thus validating its use for modeling under various vegetation covers. The substantial speedup provided by the matrix approach suggests it is advantageous for large-scale simulations, for which computational resources and execution times are critical factors. The results argue for the adoption of the matrix method in the standard version of the LSMs, which could lead to more efficient and faster calculations.

## Acknowledgements

We acknowledge CINECA and EuroHPC JU for awarding access to the Leonardo supercomputing hosted at CINECA (project ID: EHPC-DEV-2023D10-042).

## References

[1] Boucher, Olivier, et al. "Presentation and evaluation of the IPSL-CM6A-LR climate model." *Journal of Advances in Modeling Earth Systems* **12**.7 (2020): e2019MS002010. DOI: 10.1029/2019MS002010.

[2] Sellers, P. J. (1985) "Canopy reflectance, photosynthesis and transpiration." *International Journal of Remote Sensing* **6** (8): 1335-1372. DOI: 10.1080/01431168508948283.

[3] Dickinson, R. E. (1983) "Land surface processes and climate—Surface albedos and energy balance." *Advances in Geophysics* **25**: 305-353. Elsevier. DOI: 10.1016/S0065-2687(08)60176-4.

[4] Thornton, E. (2010) "Technical Description of version 4.0 of the Community Land Model (CLM)." NCAR, Climate and Global. No DOI available.

[5] Pinty, B., et al. (2006) "Simplifying the interaction of land surfaces with radiation for relating remote sensing products to climate models." *Journal of Geophysical Research: Atmospheres* **111** (D2). DOI: 10.1029/2005JD005952.

[6] Meador, W. E., and W. R. Weaver. (1980) "Two-stream approximations to radiative transfer in planetary atmospheres: A unified description of existing methods and a new improvement." *Journal of Atmospheric Sciences* **37** (3): 630-643. DOI: 10.1175/1520-0469(1980)037:0630.2.0.CO;2.

[7] McGrath, M. J., Ryder, J., Pinty, B., Otto, J., Naudts, K., Valade, A., et al. (2016) "A multi-level canopy radiative transfer scheme for ORCHIDEE (SVN r2566), based on a domain-averaged structure factor." *Geoscientific Model Development Discussions* **2016**: 1-22. DOI: 10.5194/gmd-9-3417-2016.

[8] Toon, O. B., et al. (1989) "Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres." *Journal of Geophysical Research: Atmospheres* **94** (D13): 16287-16301. DOI: 10.1029/JD094iD13p16287.

[9] Myneni, Ranga B., Juhan Ross, and Ghassem Asrar. "A review on the theory of photon transport in leaf canopies." *Agricultural and Forest Meteorology* **45** (1-2): 1-153. DOI: 10.1016/0168-1923(89)90054-7.

[10] University of East Anglia Climatic Research Unit; Harris, I.C. "CRU JRA: Collection of CRU JRA forcing datasets of gridded land surface blend of Climatic Research Unit (CRU) and Japanese reanalysis (JRA) data." *Centre for Environmental Data Analysis*, 2019.