# Spectrally refined unbiased Monte Carlo estimate of the Earth's global radiative cooling

Yaniss Nyffenegger-Péré, Raymond Armante, Mégane Bati, Stéphane Blanco, Jean-Louis Dufresne, Mouna El Hafi, Vincent Eymet, Vincent Forest, Richard Fournier, Jacques Gautrais, et al.

► To cite this version:

Yaniss Nyffenegger-Péré, Raymond Armante, Mégane Bati, Stéphane Blanco, Jean-Louis Dufresne, et al. Spectrally refined unbiased Monte Carlo estimate of the Earth's global radiative cooling. Proceedings of the National Academy of Sciences of the United States of America, 2024, 121 (5), 10.1073/pnas.2315492121. hal-04417426

HAL Id: hal-04417426

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

Submitted on 26 Jan 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_1.png" /></figure>

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

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

# Spectrally refined unbiased Monte Carlo estimate of the Earth's global radiative cooling

Yaniss Nyffenegger-Péré<sup>a,b,1</sup>, Raymond Armante<sup>c</sup>, Mégane Bati<sup>d</sup>, Stéphane Blanco<sup>a</sup>, Jean-Louis Dufresne<sup>c</sup>, Mouna El Hafi<sup>a,e</sup>, Vincent Eymet<sup>f</sup>, Vincent Forest<sup>f</sup>, Richard Fournier<sup>a</sup>, Jacques Gautrais<sup>g</sup>, Raphaël Lebrun<sup>c</sup>, Nicolas Mellado<sup>d</sup>, Nada Mourtaday<sup>a</sup>, and Mathias Paulin<sup>d</sup>

Edited by Mark Thiemens, University of California San Diego, La Jolla, CA; received October 9, 2023; accepted December 15, 2023

The Earth's radiative cooling is a key driver of climate. Determining how it is affected by greenhouse gas concentration is a core question in climate-change sciences. Due to the complexity of radiative transfer processes, current practices to estimate this cooling require the development and use of a suite of radiative transfer models whose accuracy diminishes as we move from local, instantaneous estimates to global estimates over the whole globe and over long periods of time (decades). Here, we show that recent advances in nonlinear Monte Carlo methods allow a paradigm shift: a completely unbiased estimate of the Earth's infrared cooling to space can be produced using a single model, integrating the most refined spectroscopic models of molecular gas energy transitions over a global scale and over years, all at a very low computational cost (a few seconds).

climate change | radiative forcing | line-by-line | Monte Carlo

The Earth radiative budget plays a fundamental role in the climate system. Its global value drives the global mean surface temperature changes, its latitudinal variation drives the atmosphere and ocean circulation, its value in the atmosphere drives the mean precipitation amount, etc. (1). Building an estimate of it requires the very demanding task of dealing jointly with linear transport physics and molecular spectroscopy, while handling millions of rotation/vibration lines of greenhouse gases inside a heterogeneous multiple scattering atmosphere containing clouds. One of the most widely publicized diagnostic is the instantaneous forcing of greenhouse gases (H₂O, CO₂, CH₄, etc.), defined as the sensitivity of the radiative flux at the top of the atmosphere or at the tropopause to a change in atmospheric greenhouse gas concentration, all other variables being kept constant. An accurate estimate of this sensitivity requires the same detailed physical modeling as for remote sensing, i.e., to include all the statistical physics and quantum mechanics responsible for the intensity and frequency shape of all the spectral lines. Estimates for climate studies require integration over a climatic period [t₁, t₂], over the whole globe, over all frequencies v, and over the whole height of the atmosphere (Fig. 1).

The long-term average of top-of-atmosphere flux emitted by the atmosphere toward space, assuming local thermodynamic equilibrium, is given by integrals in space, time, and electromagnetic frequency:

$$E_{\text{LW}} = \int_{t_1}^{t_2} dt \int_{0}^{+\infty} dv \int_{\text{Globe}} dx \int_{0}^{Z_T} dz \int_{\Gamma(\mathbf{x})} D\gamma 4\pi k_{d,v}(\mathbf{x}, z, t) I_v^{eq}(T(\mathbf{x}, z, t)) T(\gamma), \quad [1]$$

where x is the longitudinal/latitudinal location, z the altitude, and Γ the space of all sample paths γ starting at this location. 4πk<sub>d,v</sub>I<sub>v</sub><sup>eq</sup> is the local atmospheric emission, where k<sub>d,v</sub> is the monochromatic absorption coefficient and I<sub>v</sub><sup>eq</sup> the Planck function at the local temperature. T is either 1 if the path reaches the top of the atmosphere at altitude Z<sub>T</sub> and contributes to the top of atmosphere flux, or 0 if the path ends with absorption either by the atmosphere or the ground. The notation Dγ is used for the probabilistic measure of the multiple scattering, multiple reflection process describing the way photons interact with both the atmosphere and the ground. For a nonreflecting ground and a nonscattering stratified atmosphere, the statistics could be reduced to those of the azimuth angle cosine μ on [-1, 1] and the path length l on [0, +∞]. For example, arbitrarily extending k<sub>d,v</sub> below the ground and above the atmosphere to represent complete absorption when reaching ground or space, Dγ would typically be written as:

$$D\gamma = \frac{1}{2} k_{d,v}(\mathbf{x}, z, t) \exp \left( - \int_{0}^{l} k_{d,v}(\mathbf{x}, z + \mu l', t) dl' \right) d\mu dl. \quad [2]$$

Author affiliations: <sup>a</sup>Laboratoire Plasma et Conversion d'Energie, Université de Toulouse, CNRS, Institut national polytechnique de Toulouse (INPT), Université Paul Sabatier (UPS), Toulouse 31062, France; <sup>b</sup>Instituto de Astrofísica de Andalucía, Consejo Superior de Investigaciones Científicas (CSIC), Glorieta de la Astronomía s/n, Granada 18008, Spain; <sup>c</sup>Laboratoire de Météorologie Dynamique-Institut Pierre Simon Laplace, Sorbonne Université, École normale supérieure (ENS), Université Paris Sciences et Lettres (PSL), École polytechnique, Institut Polytechnique de Paris, CNRS, Paris 75005, France; <sup>d</sup>Institut de recherche en informatique de Toulouse, Université de Toulouse, CNRS, Toulouse 31062, France; <sup>e</sup>Institut Mines-Télécom Mines Albi, Centre de Recherche d'Albi en génie des Procédés des Solides Divisés de l'Energie et de l'Environnement, CNRS, Campus Jarlard, Albi 81013, France; <sup>f</sup>Méso-Star, Toulouse 31200, France; and <sup>g</sup>Centre de Recherches sur la Cognition Animale, Centre de Biologie Intégrative, Université de Toulouse, CNRS, Toulouse 31062, France

Author contributions: S.B., J.-L.D., R.F., N. Mellado, and M.P. designed research; Y.N.-P., R.A., M.B., J.-L.D., V.E., V.F., R.L., N. Mellado, N. Mourtaday, and M.P. performed research; Y.N.-P., R.A., M.B., J.-L.D., V.E., V.F., R.L., N. Mellado, N. Mourtaday, and M.P. contributed new reagents/analytic tools; and Y.N.-P., R.A., M.B., S.B., J.-L.D., M.E.H., V.E., V.F., R.F., J.G., R.L., N. Mellado, N. Mourtaday, and M.P. wrote the paper.

The authors declare no competing interest.

Copyright © 2024 the Author(s). Published by PNAS. This open access article is distributed under Creative Commons Attribution License 4.0 (CC BY).

¹To whom correspondence may be addressed. Email: nyffenegger@laplace.univ-tlse.fr.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas. 2315492121/-/DCSupplemental.

Published January 22, 2024.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Complexity of the integration domain. The flux density at top of atmosphere (TOA), at a given location and a given time, is a double integral: over all frequencies and over the space of multiple-scattering paths. This flux density is integrated over all TOA locations and over a climatic period. The local thermodynamic state of the atmosphere (temperature T, total pressure P, partial pressure P_t for species i) is different at each time during this climatic period, which is essentially translated into a time-varying line spectrum, associated with molecular states' transitions broadened by frequency shifts and collisions, with Doppler line shapes in the upper atmosphere (no collision) that are very significantly sharper than Lorentz shapes close to the surface where pressure-broadening dominates.</figcaption></figure>

The total absorption coefficient:

$$k_{a,v}(\mathbf{x}, z, t) = \sum_{i=1}^{N_t} h_{a,v,i}(\mathbf{x}, z, t), \quad [3]$$

is the sum of absorption due to many spectral lines, $h_{a,v,i}$, each of which corresponds to a given quantum mechanical transition provided by spectroscopic databases such as Gestion et Etude des Informations Spectroscopiques Atmosphériques (GEISA) (2) or high-resolution transmission molecular absorption database (HITRAN) (3).

Computing the multiple integrals and the transition sums in Eqs. 1–3 with a standard discretization approach is a routine practice but computation times lengthen significantly when summing over many lines at each frequency. Alternative strategies have therefore been adopted by the community by reducing climate to a limited number of representative atmospheric profiles and/or by simplifying the spectroscopic description of gas absorption (1, 4, 5).

We report here that nonlinear Monte Carlo strategies now offer a practical way of bypassing these approximations and working directly with the most refined radiative transfer models, without sacrificing any spectral, spatial, or temporal dimensions. We even show that this Monte Carlo integration is fully insensitive to each of these integration dimensions and that to estimate the global mean radiative flux at the top of the atmosphere does not require more than a few seconds on a standard laptop (Fig. 2). This is a breakthrough when compared to conventional approaches to calculating radiative forcing and opens up research avenues in atmospheric radiative heat transfer.

Why haven't such Monte Carlo simulations been performed before? The reason lies essentially in the absorption coefficient $k_{a,v}$ appearing inside the exponential of Beer's extinction law in Eq. 2. If $k_{a,v}$ were known and did not vary along the vertical axis, Monte Carlo integration in time, space, and frequency Eq. 1 would present no particular difficulty. However, $k_{a,v}$ is actually a computationally expensive sum over molecular state transitions

<figure><img src="image_5.png" /><figcaption>Fig. 2. No sensitivity of computation time to the extension of the integration domains (● without and ◦ with scattering by 1D clouds). Time to perform a Monte Carlo computation of the radiative flux at the top of the atmosphere, with an uncertainty of 0.1%, on a personal computer with 12 CPUs, as a function of the size of the main integration domains: (A) a single time and a single column but different spectral integration ranges; (B) a single time, all infrared, but different spatial integration ranges (number of columns); (C) all infrared, the whole Earth, but different climatic period integration ranges. The sky is either clear or with an idealized high level cloud to show the very limited impact of including some scattering. Note that while the variance of the estimates is unaffected by the widening of the integration domains, the estimates of radiative flux are themselves affected. Only the last green dot represents the computation time required for a climatically meaningful quantity: the mean radiative flux over all thermal infrared domain, all the globe and a 10-y period.</figcaption></figure>

that indeed vary sharply on the vertical axis, because of absorbing gas concentrations (O_3, H_2O, ...), pressure, and temperature. So the very computation of $k_{a,v}$ needs to be nested into this first-level integral through a nonlinear function (integrating the exponential of an integral). Unfortunately, when it comes to nonlinear combinations of integration spaces, Monte Carlo methods are notoriously inefficient, if not impossible to use (6, 7).

Breakthroughs have been reported in refs. 8 and 9. Here, we essentially retain the null-collision concept as a way to bypass the nonlinearity of the exponential function (10, 11). In short, the idea is to add virtual colliders that do not modify the solution of the equations such that their addition allows the true absorption coefficients $k_{a,v}$, which vary with location, to be replaced by a uniform collision coefficient, whose value is set to $\hat{k}_{a,v}$, the upper bound of the true ones. This virtual value is then used to sample the distance covered before the next collision. This distance is inevitably underestimated, and a rejection technique is used to ignore the virtual colliders. In such cases, the algorithm continues, sampling a new path length and a new probability of being absorbed, and so on until absorption.

In practice, once time, frequency, position on the globe and direction of propagation (t, v, x, z in Eq. 1 and μ in Eq. 2) have been sampled, we sample a path-length as if the photon were traveling through this upper-bound field $\hat{k}_{a,v}$ and move it to the collision location. At this stage, a standard null-collision algorithm would require defining a collision-rejection probability $P_n = (\hat{k}_{a,v} - k_{a,v})/\hat{k}_{a,v}$, where $k_{a,v}$ is a function of the local thermodynamic state of the atmosphere at this location. However, the calculation of $k_{a,v}$ represents by itself a very large computational cost due to the large number $N_t$ of molecular transitions that need to be taken into account in
Eq. 3. This is often replaced in line-by-line models by a costly pre-computation of look-up tables that are used to interpolate the value of $k_{a,v}$ at each location. In order to avoid this pre-computation step and obtain direct access to the spectroscopic data, we push the framework further to include the calculation of this sum in the Monte Carlo integration itself: In ref. 12, it was shown that when $P_n$ is an expectation, the corresponding random variable can be sampled and null collisions can be decided from the sample without ever estimating $P_n$. Here, the corresponding idea is to sample one transition and decide null collisions using this single transition. In this way, we have designed a fully unbiased spectrally refined Monte Carlo estimate: by enriching the path-space structure with orthogonal random visiting of energy transitions (Material and Methods). We also made use of machine learning techniques to further accelerate the computation of both $\hat{k}_{a,v}$ and the sampling of transitions. Here, the machine learning induces no uncontrolled source of uncertainty: These accelerations have been constructed in a strictly consistent way to ensure that the Monte Carlo estimate remains rigorously unbiased (SI Appendix).

We validated the fluxes against the state-of-the-art 4A-Flux line-by-line model (13) for a set of standard atmospheric columns. Using the 100 atmospheric profiles proposed by Pincus et al. (5), we obtained a radiative forcing of $2.73 \pm 0.06$ W m$^{-2}$ for a doubling of the CO$_2$ concentration by doing the difference between two mean fluxes, consistent with their mean estimate of $2.71$ W m$^{-2}$. We then went one step further by performing a computation over the whole globe and over climatic periods. To achieve this, the atmospheric columns were first randomly sampled from three hourly outputs of a 10-y simulation of the Institut Pierre-Simon Laplace (IPSL) climate model (14) before the radiative calculation was performed. We obtained a global average clear-sky radiative forcing of $2.56 \pm 0.06$ W m$^{-2}$.

The calculations are very fast as it takes only a few seconds on a standard personal computer to estimate this flux, with a statistical error of 0.1%, obtained by launching about $N = 10^7$ samples. As expected for a converging Monte Carlo code, the statistical error decreases with the square root of $N$ and the computation time cost lies in the integration including the largest source of variance and does not increase when extending the size of the other integration domains (Fig. 2). For a given atmospheric column and a given relative uncertainty, the computation time required to estimate the radiative flux does not increase as the spectral range of integration varies from $6$ cm$^{-1}$ to the entire

spectral range (Fig. 2A). Likewise, for a given spectral domain and a given relative uncertainty, the computation time does not increase when increasing the spatial (Fig. 2B) and the time (Fig. 2C) integration domains (SI Appendix).

We have thus devised a line-by-line Monte Carlo integration of TOA radiative flux at the scale of the Earth and for climatic periods, which is insensitive to the widening of frequency, spatial, and temporal integration domains, and which is valid for a scattering or nonscattering medium. Freeing the climate science community from the need for pre-computed approximation schemes, these estimates are now available in a matter of seconds on a desktop computer, thus opening up a wide range of perspectives in climate science (SI Appendix).

## Materials and Methods

Starting from the radiative transfer equation (with pure absorption for the sake of simplicity),

$$ \mathbf{u} \cdot \nabla l = -kl + kl^{eq}, $$

where $l \equiv l(\mathbf{x}, z, \mathbf{u})$ is the specific intensity in direction $\mathbf{u}$, the simplest viewpoint on null-collision (11) is to consider virtual colliders as pure forward scatterers:

$$ \mathbf{u} \cdot \nabla l = -\hat{k}l + \hat{k} \left[ (1-\omega)l^{eq} + \omega \int_{4\pi} \delta(\mathbf{u}-\mathbf{u}')l'd\mathbf{u} \right] $$

with $l' \equiv l(\mathbf{x}, z, \mathbf{u}')$ and single scattering albedo $\omega = \frac{\hat{k}-k}{k}$, the Kronecker symbol $\delta$ representing the pure forward scattering phase function. Looking at each transition $i$ as if it were an independent species, with $k = \sum h_i$, virtual colliders can be added to each transition, with $\hat{h}_i > h_i$ to write:

$$ \mathbf{u} \cdot \nabla l = -\hat{k}l + \hat{k} \sum P_i \left[ (1-\omega_i)l^{eq} + \omega_i \int_{4\pi} \delta(\mathbf{u}-\mathbf{u}')l'd\mathbf{u} \right] $$

with $\hat{k} = \sum \hat{h}_i$, $P_i = \frac{\hat{h}_i}{k}$, and $\omega_i = \frac{\hat{h}_i-\hat{h}_i}{\hat{h}_i}$. As for standard Monte Carlo algorithms involving mixtures, $\hat{k}$ is used to find the next collision location, $P_i$ to decide which transition was responsible for the collision, and $\omega_i$ to decide whether the collision is an absorption or a scattering event (here a rejection).

**Data, Materials, and Software Availability.** Code data have been deposited in RadForce (https://gitlab.com/yanissnp/RadForce) (15).

**ACKNOWLEDGMENTS.** This work was supported by the "Monte-Carlo Global Radiative Forcings Computation" (ANR-18-CE46-0012) project and by the "Centre national d'études spatiales".

1. P. Forster et al., "The Earth's energy budget, climate feedbacks, and climate sensitivity" in *Climate Change 2021: The Scientific Basis*, V. Masson-Delmotte et al., Eds. (Cambridge University Press, 2021). chap. 7, pp. 923-1054.

2. N. Jacquet-Husson et al., The 2015 edition of the GEISA spectroscopic database. *J. Molec. Spectrosc.* **327**, 1-72 (2016).

3. I. E. Gordon et al., The HITRAN2016 molecular spectroscopic database. *J. Quant. Spectrosc. Radiat. Transf.* **203**, 3-69 (2017).

4. G. Myhre, E. J. Highwood, K. P. Shine, F. Stordal, New estimates of radiative forcing due to well mixed greenhouse gases. *Geophys. Res. Lett.* **25**, 2715-2718 (1998).

5. R. Pincus et al., Benchmark calculations of radiative forcing by greenhouse gases. *J. Geophys. Res. Atm.* **125**, e2020JD033483 (2020).

6. J. H. Curtiss, "Monte Carlo" methods for the iteration of linear operators. *J. Math. Phys.* **32**, 209-232 (1953).

7. M. H. Kalos, P. A. Whitlock, *Monte Carlo Methods* (John Wiley & Sons, 2009).

8. I. T. Dimov, *Monte Carlo Methods for Applied Scientists* (World Scientific, 2008).

9. J. Dauchet et al., Addressing nonlinearities in Monte Carlo. *Sci. Rep.* **8**, 1-11 (2018).

10. M. Galtier et al., Integral formulation of null-collision Monte Carlo algorithms. *J. Quant. Spectrosc. Radiat. Transf.* **125**, 57-68 (2013).

11. M. El Hafi et al., Three viewpoints on null-collision Monte Carlo algorithms. *J. Quant. Spectrosc. Radiat. Transf.* **260**, 107402 (2021).

12. G. Terrée et al., Addressing the gas kinetics Boltzmann equation with branching-path statistics. *Phys. Rev. E* **105**, 025305 (2022).

13. Y. Tellier, C. Crevoisier, R. Armante, J.-L. Dufresne, N. Meilhac, Computation of longwave radiative flux and vertical heating rate with 4A-Flux v1.0 as an integral part of the radiative transfer code 4A/OP v1.5. *Geosci. Model Dev.* **15**, 5211-5231 (2022).

14. O. Boucher et al., Presentation and evaluation of the IPSL-CM6A-LR climate model. *J. Adv. Model. Earth Syst.* **12**, e2019MS002010 (2020).

15. Y. Nyffenegger-Péré et al., Radiative Forcings Estimator. RadForce. https://gitlab.com/yanissnp/RadForce. Deposited 1 July 2023.