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

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

## Abstract

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).

### Sign up for PNAS alerts.

Get alerts for new articles, or get an alert when an article is cited.

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 (${\text{H}}_{2}\text{O}$, ${\text{CO}}_{2}$, ${\text{CH}}_{4}$, 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}_{1}$,${t}_{2}]$, over the whole globe, over all frequencies $\nu $, and over the whole height of the atmosphere (Fig. 1).

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:where $\mathit{x}$ is the longitudinal/latitudinal location, $z$ the altitude, and $\mathrm{\Gamma}$ the space of all sample paths $\gamma $ starting at this location. $4\pi {k}_{a,\nu}{I}_{\nu}^{\mathit{eq}}$ is the local atmospheric emission, where ${k}_{a,\nu}$ is the monochromatic absorption coefficient and ${I}_{\nu}^{\mathit{eq}}$ the Planck function at the local temperature. $\mathcal{T}$ is either $1$ if the path reaches the top of the atmosphere at altitude ${Z}_{T}$ 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 $\mathcal{D}\gamma $ 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 $\mu $ on $[-1,1]$ and the path length $l$ on $[0,+\infty [$. For example, arbitrarily extending ${k}_{a,\nu}$ below the ground and above the atmosphere to represent complete absorption when reaching ground or space, $D\gamma $ would typically be written as:

$${E}_{LW}={\displaystyle {\int}_{{t}_{1}}^{{t}_{2}}dt}{\displaystyle {\int}_{0}^{+\infty}d\nu}{\displaystyle {\int}_{Globe}dx}{\displaystyle {\int}_{0}^{{Z}_{T}}dz}{\displaystyle {\int}_{\Gamma (x)}\mathcal{D}\gamma 4\pi {k}_{a,\nu}(x,z,t){I}_{\nu}^{eq}(T(x,z,t))\mathcal{T}(\gamma ),}$$

[1]

$$\begin{array}{c}\hfill \mathcal{D}\gamma =\frac{1}{2}{k}_{a,\nu}(\mathit{x},z,t)exp\left(-{\int}_{0}^{l}{k}_{a,\nu}(\mathit{x},z+\mu {l}^{\prime},t)d{l}^{\prime}\right)d\mu dl.\end{array}$$

[2]

The total absorption coefficient:

$$\begin{array}{c}\hfill {k}_{a,\nu}(\mathit{x},z,t)={\displaystyle \sum _{i=1}^{{N}_{t}}}{h}_{a,\nu ,i}(\mathit{x},z,t),\end{array}$$

[3]

is the sum of absorption due to many spectral lines, ${h}_{a,\nu ,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.

Fig. 2.

Why haven’t such Monte Carlo simulations been performed before? The reason lies essentially in the absorption coefficient ${k}_{a,\nu}$ appearing inside the exponential of Beer’s extinction law in Eq.

**2**. If ${k}_{a,\nu}$ 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,\nu}$ is actually a computationally expensive sum over molecular state transitions that indeed vary sharply on the vertical axis, because of absorbing gas concentrations (${\text{O}}_{3}$, ${\text{H}}_{2}\text{O}$, …), pressure, and temperature. So the very computation of ${k}_{a,\nu}$ 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,\nu}$, which vary with location, to be replaced by a uniform collision coefficient, whose value is set to ${\widehat{k}}_{a,\nu}$, 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$, $\nu $, $\mathit{x}$, $z$ in Eq.

**1**and $\mu $ in Eq.**2**) have been sampled, we sample a path-length as if the photon were traveling through this upper-bound field ${\widehat{k}}_{a,\nu}$ and move it to the collision location. At this stage, a standard null-collision algorithm would require defining a collision-rejection probability ${P}_{n}=({\widehat{k}}_{a,\nu}-{k}_{a,\nu})/{\widehat{k}}_{a,\nu}$, where ${k}_{a,\nu}$ is a function of the local thermodynamic state of the atmosphere at this location. However, the calculation of ${k}_{a,\nu}$ 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,\nu}$ 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 ${\widehat{k}}_{a,\nu}$ 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\phantom{\rule{0.166667em}{0ex}}\mathrm{W}\phantom{\rule{0.166667em}{0ex}}{\mathrm{m}}^{-2}$ for a doubling of the ${\text{CO}}_{2}$ concentration by doing the difference between two mean fluxes, consistent with their mean estimate of $2.71\phantom{\rule{0.166667em}{0ex}}\mathrm{W}\phantom{\rule{0.166667em}{0ex}}{\mathrm{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\phantom{\rule{0.166667em}{0ex}}\mathrm{W}\phantom{\rule{0.166667em}{0ex}}{\mathrm{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. 2*A*). Likewise, for a given spectral domain and a given relative uncertainty, the computation time does not increase when increasing the spatial (Fig. 2*B*) and the time (Fig. 2*C*) 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),

$$\begin{array}{c}\hfill \mathit{u}.\mathbf{\nabla}I=-kI+k{I}^{\mathit{eq}},\end{array}$$

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

$$\begin{array}{c}\hfill \mathit{u}.\mathbf{\nabla}I=-\widehat{k}I+\widehat{k}\left[(1-\omega ){I}^{\mathit{eq}}+\omega {\int}_{4\pi}\delta (\mathit{u}-{\mathit{u}}^{\mathbf{\prime}}){I}^{\prime}d\mathit{u}\right]\end{array}$$

with ${I}^{\prime}\equiv I(\mathit{x},z,{\mathit{u}}^{\prime})$ and single scattering albedo $\omega =\frac{\widehat{k}-k}{\widehat{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 ${\widehat{h}}_{i}>{h}_{i}$ to write:

$$\begin{array}{c}\hfill \mathit{u}.\mathbf{\nabla}I=-\widehat{k}I+\widehat{k}\sum {P}_{i}\left[(1-{\omega}_{i}){I}^{\mathit{eq}}+{\omega}_{i}{\int}_{4\pi}\delta (\mathit{u}-{\mathit{u}}^{\mathbf{\prime}}){I}^{\prime}d\mathit{u}\right]\end{array}$$

with $\widehat{k}=\sum {\widehat{h}}_{i}$, ${P}_{i}=\frac{{\widehat{h}}_{i}}{\widehat{k}}$, and ${\omega}_{i}=\frac{{\widehat{h}}_{i}-{h}_{i}}{{\widehat{h}}_{i}}$. As for standard Monte Carlo algorithms involving mixtures, $\widehat{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”.

### 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.

### Competing interests

The authors declare no competing interest.

## Supporting Information

Appendix 01 (PDF)

- Download
- 129.67 KB

## References

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. Jacquinet-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.## Information & Authors

### Information

#### Published in

#### Classifications

#### Copyright

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

#### Data, Materials, and Software Availability

Code data have been deposited in RadForcE (https://gitlab.com/yanissnp/RadForcE) (15).

#### Submission history

**Received**: October 9, 2023

**Accepted**: December 15, 2023

**Published online**: January 22, 2024

**Published in issue**: January 30, 2024

#### Keywords

#### 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”.

##### 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.

##### Competing Interests

The authors declare no competing interest.

### Authors

## Metrics & Citations

### Metrics

#### Citation statements

#### Altmetrics

### Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

## View Options

### View options

#### PDF format

Download this article as a PDF file

DOWNLOAD PDF### Get Access

#### Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login#### Recommend to a librarian

Recommend PNAS to a Librarian#### Purchase options

Purchase this article to get full access to it.