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

Based on advanced spectroscopic databases, line-by-line and layer-by-layer radiative transfer codes numerically solve the radiative transfer equation with very high accuracy. Taking advantage of its pre-calculated optical depth lookup table, the fast and accurate radiative transfer model Automatized Atmospheric Absorption Atlas OPerational (4A/OP) calculates the transmission and radiance spectra for a user-defined layered atmospheric model. Here, we present a module called 4A-Flux, which is developed and implemented into 4A/OP in order to include the calculation of the clear-sky longwave radiative flux profiles and heating rate profiles at a very high spectral resolution. Calculations are performed under the assumption of local thermodynamic equilibrium, a plane-parallel atmosphere, and specular reflection on the surface. The computation takes advantage of pre-tabulated exponential integral functions that are used instead of a classic angular quadrature. Furthermore, the sub-layer variation of the Planck function is implemented to better represent the emission of layers with a high optical depth. Thanks to the implementation of 4A-Flux, 4A/OP models have participated in the Radiative Forcing Model Intercomparison Project (RFMIP-IRF) along with other state-of-the-art radiative transfer models. 4A/OP hemispheric flux profiles are compared to other models over the 1800 representative atmospheric situations of RFMIP, yielding an outgoing longwave radiation (OLR) mean difference between 4A/OP and other models of −0.148 W m⁻² and a standard deviation of 0.218 W m⁻², showing a good agreement between 4A/OP and other models. 4A/OP is applied to the Thermodynamic Initial Guess Retrieval (TIGR) atmospheric database to analyze the response of the OLR and vertical heating rate to several perturbations of temperature or gas concentration. This work shows that 4A/OP with 4A-Flux module can successfully be used to simulate accurate flux and heating rate profiles and provide useful sensitivity studies including sensitivities to minor trace gases such as HFC134a, HCFC22, and CFC113. We also highlight the interest for the modeling community to extend intercomparison between models to comparisons between spectroscopic databases and modeling to improve the confidence in model simulations.

1 Introduction

Atmospheric radiative transfer is the main driver of the climate system and plays a central role in many atmospheric processes (Hartmann et al., 1986; Trenberth et al., 2009). Accurate models that calculate radiative variables such as optical depths, transmittance, radiance, radiative fluxes (irradiances), or vertical heating rates are currently used for multiple purposes. A certain family of algorithms called line-by-line models are commonly used as forward models in the retrieval of geophysical parameters from indirect measurements. Such codes are also needed in the development and exploitation of Earth observation satellite missions. Another family of algorithms called parameterized models are central to general circulation models (GCMs) as the energy balance between net shortwave and net longwave radiation fundamentally drives the climate system. Thus, radiative transfer models participate in the improvement of our understanding of the atmosphere and the climate.

By resolving the individual spectral lines for every absorbing species at every layer of the atmosphere, line-by-line and layer-by-layer radiative transfer codes offer a very accurate way to compute radiative transfer. Their precision relies on spectroscopic databases such as GEISA¹ (Jacquinet-Husson et al., 2016, 1999) or HITRAN² (Rothman et al., 2013) that contain the parameters of absorption lines. Being resource-intensive, these line-by-line calculations are not used in GCMs that instead use fast models based on a radiation parameterization that approximates the radiative transfer. However, thanks to its high accuracy, line-by-line radiative transfer codes can be, and have been, used as a reference to improve the performances of GCMs by improving their radiation parameterization (Buehler et al., 2006).

Several line-by-line models have recently participated in numerical experiments within the Radiative Forcing Model Intercomparison Project (RFMIP, see Pincus et al., 2016), endorsed by the sixth phase of the Coupled Model Intercomparison Project (CMIP6, see Eyring et al., 2016), that aimed at characterizing the clear-sky GCM radiation parameterization error. Pincus et al. (2020) have presented the clear-sky instantaneous radiative forcing (IRF) by greenhouse gases (GHGs) computed with these six benchmark models including the line-by-line radiative transfer code 4A/OP.

Many line-by-line radiative transfer models have been developed to provide high-accuracy radiative flux and heating rates (Buehler et al., 2006, and references therein). Among these codes we can cite LBLRTM (the Line-by-line Radiative Transfer Model – Clough et al., 2005; Clough and Iacono, 1995), ARTS (the Atmospheric Radiative Transfer Simulator – Buehler et al., 2018; Eriksson et al., 2011), and RFM (the Reference Forward Model – Dudhia, 2017). These codes have been used for many purposes such as, among many other applications, retrieval of atmospheric and surface parameters from indirect measurements (Dogniaux et al., 2020; Alvarado et al., 2013) and sensitivity studies (Clough et al., 1992; Huang et al., 2007).

The Automatized Atmospheric Absorption Atlas OPerational (4A/OP) is a fast and accurate line-by-line radiative transfer model particularly efficient in the infrared region of the spectrum. It offers the calculation of the transmittance, the Jacobians, and the radiance from a user-defined atmospheric and surface description. The concept of 4A/OP (Scott and Chedin, 1981; Cheruy et al., 1995) relies on a compressed lookup table of optical depths, called the atlas, that are computed from reference atmospheres using STRANSAC, the complete line-by-line and layer-by-layer model developed at the Laboratoire de Météorologie Dynamique (Scott, 1974; Tournier et al., 1995). 4A/OP interpolates the atlas to compute the transmittance and the radiance at the correct layer temperature and pressure levels as well as for the desired absorber composition and viewing angle.

4A/OP has a long history of validation within the framework of the international radiative transfer community. Numerous intercomparison exercises have taken place, in particular in the framework of the Intercomparison of Transmittance and Radiance Algorithms working groups of the International Radiation Commission (Chédin et al., 1988) and during the ICRCCM (Intercomparison of Radiation Codes in Climate Models) campaigns (Luther et al., 1988). However, these intercomparisons did not include the modeling of outgoing longwave radiation (OLR) and heating rates. The first developments of the computation of OLR and vertical heating rates using 4A/OP were performed in the 1990s (Chéruy et al., 1996; Chevallier et al., 1998, 2000). The computation took advantage of the readily available radiances computed by 4A/OP to calculate the radiative quantities by a spectral and an angular integration. At the time, for computation efficiency reasons, the angular quadratures were often performed after a spectral integration of the radiance spectra, which lead to errors in the modeled fluxes. The angular integration was performed with either a Gaussian quadrature or at a single angle under the diffuse approximation. Here, different integration approaches have been used. Firstly, the spectral integration is systematically performed at the finest resolution available. And secondly, exponential integral functions applied to optical depth to directly compute the fluxes have been chosen instead of the quadrature of the radiances.

The main objective of this article is to describe the new implementation of the calculation of the clear-sky longwave radiative flux profiles and vertical cooling rates into 4A/OP radiative transfer code. This module, called 4A-Flux, takes advantage of the efficiency and accuracy of 4A/OP to compute the radiative quantities. Here, we first describe the implementation of 4A-Flux into 4A/OP. The radiative flux profiles are then compared to the outputs of the five other radiative transfer models that have contributed to RFMIP. Finally, we have applied the 4A-Flux module to the computation of the outgoing longwave radiation (OLR) and vertical heating rate sensitivities to several surface and atmospheric parameters.

2 Implementation of 4A-Flux into the 4A/OP radiative transfer model

2.1 General approach to flux computation

The 4A-Flux module computes the clear-sky longwave downwelling, upwelling, and net radiative fluxes (or irradiance) as well as the vertical heating rate at every pressure level of the user-defined atmospheric model and at the finest spectral resolution offered by the line-by-line radiative transfer model 4A/OP of 5 × 10⁻⁴ cm⁻¹ (Scott and Chedin, 1981). The 4A/OP model offers a standard atlas that extends from 0 to 14 000 cm⁻¹. As the model can handle any atlas, this spectral range can be extended by generating a new atlas in the required regions. In the context of the development of 4A-Flux that computes radiative flux, OLR, and heating rates without the contribution of the solar flux, we have concentrated on the spectral range from 10 to 3250 cm⁻¹.

Figure 1 4A/OP and 4A-Flux flowchart. P, T, and ρ are the pressure grid, the temperature profile, and the profiles of gas concentrations, respectively.

Figure 1 presents the flowchart of 4A/OP with its new module 4A-Flux. The 4A/OP model provides a suitable basis for the computation of the radiative fluxes as it computes the vertical distribution of total optical depths at the spectral resolution of absorption lines. The inputs of 4A/OP are the atlases, a user-defined description of the pressure grid, the vertical distribution of temperature and gas concentrations, and surface properties such as surface temperature and spectral emittance. In this work, the spectroscopic database used in the calculation of the atlases is GEISA-2015 (Jacquinet-Husson et al., 2016); the continuum used is the MT_CKD 3.2 continuum (Mlawer et al., 2012), and the CO₂ line mixing is the Lamouroux et al. (2015). Note that Lamouroux et al. (2015) is based on the parameters of HITRAN2012. Ongoing efforts will soon lead to a CO₂ line mixing that will be based on the 2020 edition of the GEISA database (Delahaye et al., 2021). The pressure grid, the vertical distribution of temperature and gas concentrations, and surface properties such as surface temperature and spectral emittance can be specified by the user. Based on this description of the atmosphere and optical properties of gases, 4A-Flux computes the upwelling and downwelling monochromatic radiative flux profiles, also referred to as hemispherical monochromatic flux profiles [unit: $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}({\mathrm{cm}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$]. This computation is based on an exponential integral evaluation to perform the angular integration and a vertical integration. In the distributed version of 4A/OP, the vertical integration of the radiative transfer equation is performed under the assumption of isothermal layers (emission at the mean layer temperature), whereas in 4A-Flux, to better represent the real emission of each layer, we have implemented a sub-layer Planck function that varies linearly with the optical depth. 4A-Flux also offers the option of computing the monochromatic outgoing longwave radiation (OLR) instead of the full vertical description of fluxes, allowing a gain of computational time whenever the user does not require the complete vertical description of fluxes. The output spectral resolution is set by default to 1 cm⁻¹ but can be specified by the user at any other value. The net monochromatic flux profile and the monochromatic heating rate vertical profile are then derived from the hemispheric profiles. The following section defines the radiative quantities and all the assumptions that have been made. We first focus on the definition of the downwelling flux, and then we treat the upwelling flux; finally, all the derived quantities that are computed by 4A-Flux are defined.

2.2 The computation of fluxes and heating rate profiles

The monochromatic³ hemispheric fluxes [unit: $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}({\mathrm{cm}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$], denoted F^↓(σ) for the downwelling flux and F^↑(σ) for the upwelling flux, are represented here by functions of σ, the monochromatic optical depth [unit: dimensionless] for all accounted gases in the model, as seen from the top of the atmosphere (TOA), where σ is null, to the bottom of the atmosphere (BOA) where it reaches its maximum value σ_BOA. Here, σ serves as a vertical coordinate as it increases monotonically with the pressure. The hemispheric fluxes are defined as the integrals of the monochromatic hemispheric radiances $I^{\downarrow }(\mathit{\sigma },\mathit{\mu },\mathit{\varphi })$ and $I^{\uparrow }(\mathit{\sigma },\mathit{\mu },\mathit{\varphi })$ [unit: $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{sr}}^{-\mathrm{1}}({\mathrm{cm}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$] over the respective hemispheres (Stamnes et al., 2017; Liou, 2002). These are defined as functions of the vertical coordinate σ, the azimuthal angle ϕ, and the direction μ=|cos (θ)|, which is the absolute value of the cosine of the zenith angle θ, where $\mathit{\theta }>\mathit{\pi }/\mathrm{2}$ for the downwelling flux and $\mathit{\theta }\le \mathit{\pi }/\mathrm{2}$ for the upwelling flux:

$$\begin{array}{}\text{(1)}& F^{\downarrow \uparrow }\left(\mathit{\sigma }\right)=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathit{\mu }{I}^{\downarrow \uparrow }(\mathit{\sigma },\mathit{\mu },\mathit{\varphi })\mathrm{d}\mathit{\mu }\mathrm{d}\mathit{\varphi }.\end{array}$$

As our computations are performed on the longwave spectrum and for clear-sky conditions, we assume the atmosphere to be non-scattering, and we neglect the radiance of the sun. We also assume an isotropic thermal emission of gases. Under such assumptions, the radiance is independent of the azimuthal angle, and thus we replace $I^{\downarrow \uparrow }(\mathit{\sigma },\mathit{\mu },\mathit{\varphi })$ by $I^{\downarrow \uparrow }(\mathit{\sigma },\mathit{\mu })$. By integration with respect to ϕ we obtain

$$\begin{array}{}\text{(2)}& F^{\downarrow \uparrow }\left(\mathit{\sigma }\right)=\mathrm{2}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathit{\mu }{I}^{\downarrow \uparrow }(\mathit{\sigma },\mathit{\mu })\mathrm{d}\mathit{\mu }.\end{array}$$

2.2.1 Downwelling flux

We focus on the downwelling flux in this section. Under the hypotheses of local thermodynamic equilibrium and local stratification of the atmosphere (plane-parallel geometry), the radiative transfer equation can be solved for the radiance. As we focus on the longwave only, the radiative contribution of deep space (less than 3 K) to the downwelling flux is negligible compared to the atmospheric contribution. Thus, the only remaining term in the solution of the radiative transfer equation is the atmospheric component, commonly referred to as the Schwarzschild equation:

$$\begin{array}{}\text{(3)}& I^{\downarrow }(\mathit{\sigma },\mathit{\mu })=\underset{\mathrm{0}}{\overset{\mathit{\sigma }}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\mu }}}B\left[T\right({\mathit{\sigma }}^{\prime }\left)\right]\exp [-(\mathit{\sigma }-{\mathit{\sigma }}^{\prime })/\mathit{\mu }]\mathrm{d}{\mathit{\sigma }}^{\prime }.\end{array}$$

Thus, by inserting Eq. (3) into Eq. (2), we obtain

$$\begin{array}{}\text{(4)}& F^{\downarrow }\left(\mathit{\sigma }\right)=\mathrm{2}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\sigma }}{\int }}B\left[T\right({\mathit{\sigma }}^{\prime }\left)\right]\exp [-(\mathit{\sigma }-{\mathit{\sigma }}^{\prime })/\mathit{\mu }]\mathrm{d}{\mathit{\sigma }}^{\prime }\mathrm{d}\mathit{\mu },\end{array}$$

where B[T(σ^′)] is the Planck function [unit: $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{sr}}^{-\mathrm{1}}({\mathrm{cm}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$] evaluated at wavenumber ν (omitted in the notation) and at temperature T(σ^′) [unit: K]. By switching the integration order, we can rewrite the downwelling flux as

$$\begin{array}{}\text{(5)}& F^{\downarrow }\left(\mathit{\sigma }\right)=\mathrm{2}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathit{\sigma }}{\int }}B\left[T\right({\mathit{\sigma }}^{\prime }\left)\right]E_{\mathrm{2}}(\mathit{\sigma }-{\mathit{\sigma }}^{\prime })\mathrm{d}{\mathit{\sigma }}^{\prime }.\end{array}$$

Here, we have used the exponential integral of order n defined for n>0 and x≥0 as

$$\begin{array}{}\text{(6)}& E_n\left(x\right)=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\mu }}^{n-\mathrm{2}}{e}^{-x/\mathit{\mu }}\mathrm{d}\mathit{\mu }.\end{array}$$

In addition to the notation simplifications they allow, these series of functions present convenient mathematical properties. Their values at the origin are

$$\begin{array}{}\text{(7)}& E_n\left(\mathrm{0}\right)=\left\{\begin{array}{ll}+\mathrm{\infty }& \text{for}n=\mathrm{1},\\ \frac{\mathrm{1}}{n-\mathrm{1}}& \text{for}n>\mathrm{1}.\end{array}\right.\end{array}$$

They also follow simple recurrence relations:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}{E}_n\left(x\right)=\left\{\begin{array}{ll}-e^{-x}/x& \text{for}n=\mathrm{1},\\ -E_{n-\mathrm{1}}\left(x\right)& \text{for}n>\mathrm{1}.\end{array}\right.\end{array}$$

This formulation is commonly used to compute fluxes (Stamnes et al., 2017; Feigelson et al., 1991; Ridgway et al., 1991). The exponential integral functions can be evaluated numerically even though the computational time is significant in the context of a line-by-line flux calculation. To overcome this limitation, 4A-Flux uses a tabulation of the E_n functions that is pre-calculated using the GNU Scientific Library (https://www.gnu.org/software/gsl/, last access: 12 June 2022).

In order to compute the downwelling flux for the discretized model of the atmosphere used in the 4A/OP model, we decompose the vertical integral into sums over every layer of the atmosphere. The atmospheric column is decomposed into N+1 levels using subscripts denoted k ranging from 0 to N or N layers using subscripts denoted $k-\frac{\mathrm{1}}{\mathrm{2}}$ with k ranging from 1 to N. We denote the value of the optical depth between the TOA and any level k as σ_k. At the TOA σ₀=0, and at the surface σ_N=σ_BOA. At the TOA, the downwelling flux is supposed null, $F^{\downarrow }\left({\mathit{\sigma }}_{\mathrm{0}}\right)=\mathrm{0}$, and for every level k>0, we have

$$\begin{array}{}\text{(9)}& F^{\downarrow }\left({\mathit{\sigma }}_k\right)=\mathrm{2}\mathit{\pi }\sum _{i=\mathrm{1}}^k\underset{{\mathit{\sigma }}_{i-\mathrm{1}}}{\overset{{\mathit{\sigma }}_i}{\int }}B\left[T\right({\mathit{\sigma }}^{\prime }\left)\right]E_{\mathrm{2}}({\mathit{\sigma }}_k-{\mathit{\sigma }}^{\prime })\mathrm{d}{\mathit{\sigma }}^{\prime }.\end{array}$$

To compute the atmospheric contributions of each layer to the fluxes (the integrand of Eq. 9) an assumption has to be made about the variation of the Planck function at sub-layer scale. The simplest assumption considers the Planck function to be uniform at the sub-layer scale (isothermal layer) and evaluated at the average temperature of the layer $B\left[T\right({\mathit{\sigma }}^{\prime }\left)\right]=B\left[{\stackrel{\mathrm{‾}}{T}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}\right]$, where ${\stackrel{\mathrm{‾}}{T}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}=\left[T\right({\mathit{\sigma }}_i)+T({\mathit{\sigma }}_{i-\mathrm{1}}\left)\right]/\mathrm{2}$. This assumption is a rather coarse approximation of the effective emission of the layer, especially when the optical depth is important. In such a case, the effective temperature is closer to the lower (upper) boundary of the layer for the downwelling (upwelling flux). A better assumption should take into account that the effective emission temperature is related to the optical depth of the layer. Thus, in our approach we have considered a sub-layer variation of the Planck function as linear in optical depth (also referred to as B(σ) linear). According to Ridgway et al. (1991), this approach dates back to Schuster (1905). Under such an assumption, the sub-layer variation of the Planck function in the layer $i-\frac{\mathrm{1}}{\mathrm{2}}$ and for the downwelling flux at level i is approximated by

$$\begin{array}{}\text{(10)}& B\left[T\right({\mathit{\sigma }}^{\prime }\left)\right]=B\left[T\right({\mathit{\sigma }}_i\left)\right]+\left[{\displaystyle \frac{\mathrm{2}\left(B\right[{\stackrel{\mathrm{‾}}{T}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}]-B[T\left({\mathit{\sigma }}_i\right)\left]\right)}{{\mathit{\sigma }}_i-{\mathit{\sigma }}_{i-\mathrm{1}}}}\right]({\mathit{\sigma }}_i-{\mathit{\sigma }}^{\prime }).\end{array}$$

The assumption of an isothermal layer would lead to a systematic cold (warm) bias in estimated downwelling (upwelling) fluxes (Ridgway et al., 1991), whereas the B(σ) linear assumption substantially reduces this bias (the remaining bias is discussed by Wiscombe, 1976). To better understand this approximation, one can verify that Eq. (10) is equivalent to the isothermal layer assumption when the layer optical depth tends to zero and is equivalent to the Planck function taken at the lower (upper) layer boundary temperature for the downwelling (upwelling) flux when the layer optical depth tends to infinity. For a more detailed explanation the reader can refer to Clough et al. (1992). Thanks to the B(σ) linear approximation, the downwelling flux can be integrated. To simplify the notations, we introduce the optical depth between two levels i and j as ${\mathrm{\Delta }}_{i,j}={\mathit{\sigma }}_i-{\mathit{\sigma }}_j$ (for $\mathrm{0}\le i,j\le N$). The downwelling flux, for every level k>0, can finally be expressed as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{F}^{\downarrow }\left({\mathit{\sigma }}_k\right)=& \mathrm{2}\mathit{\pi }\sum _{i=\mathrm{1}}^{k}B\left[T\right({\mathit{\sigma }}_i\left)\right]\left[E_{\mathrm{3}}\right({\mathrm{\Delta }}_{k,i})-E_{\mathrm{3}}({\mathrm{\Delta }}_{k,i-\mathrm{1}}\left)\right]\\ & +\mathrm{2}\mathit{\pi }\sum _{i=\mathrm{1}}^k\mathrm{2}\left(B\right[{\stackrel{\mathrm{‾}}{T}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}]-B[T\left({\mathit{\sigma }}_i\right)\left]\right)\\ & \times \left[{\displaystyle \frac{E_{\mathrm{4}}\left({\mathrm{\Delta }}_{k,i}\right)-E_{\mathrm{4}}\left({\mathrm{\Delta }}_{k,i-\mathrm{1}}\right)}{{\mathrm{\Delta }}_{i,i-\mathrm{1}}}}-E_{\mathrm{3}}\left({\mathrm{\Delta }}_{k,i-\mathrm{1}}\right)\right].\end{array}\end{array}$$

The first sum in Eq. (11) can be viewed as the cumulative downwelling flux emitted at the temperature of the lower boundary of every layer, whereas the second sum brings the correction needed to account for the fact that the equivalent temperature of emission is somewhere between lower boundary temperature T(σ_i) and the average layer temperature ${\stackrel{\mathrm{‾}}{T}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}$ depending on the value of the optical depth of the layer ${\mathrm{\Delta }}_{i,i-\mathrm{1}}$.

2.2.2 Upwelling flux

Now that the downwelling flux has been treated, we will focus on the monochromatic upwelling flux. Unlike for the downwelling flux of Eq. (4) where the downwelling radiance only accounts for the atmospheric component, the upwelling radiance must also account for the surface emission and the reflection of the downwelling flux on the surface. Thus, the upwelling flux is written as

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{F}^{\uparrow }\left(\mathit{\sigma }\right)=& \stackrel{\text{atmospheric\hspace{0.17em} component}}{\overbrace{\mathrm{2}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\underset{\mathit{\sigma }}{\overset{{\mathit{\sigma }}_{\text{BOA}}}{\int }}B\left[T\right({\mathit{\sigma }}^{\prime }\left)\right]\exp [-({\mathit{\sigma }}^{\prime }-\mathit{\sigma })/\mathit{\mu }]\mathrm{d}{\mathit{\sigma }}^{\prime }\mathrm{d}\mathit{\mu }}}\\ & +\stackrel{\text{surface\hspace{0.17em} component}}{\overbrace{\mathrm{2}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathit{\mu }\mathit{ϵ}B\left(T_{\mathrm{s}}\right)\exp [-({\mathit{\sigma }}_{\text{BOA}}-\mathit{\sigma })/\mathit{\mu }]\mathrm{d}\mathit{\mu }}}\\ & +\stackrel{\text{reflected\hspace{0.17em} component}}{\overbrace{\mathrm{2}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathit{\mu }\left\{\begin{array}{c}(\mathrm{1}-\mathit{ϵ})I^{\downarrow }({\mathit{\sigma }}_{\text{BOA}},\mathit{\mu })\\ \times \exp [-({\mathit{\sigma }}_{\text{BOA}}-\mathit{\sigma })/\mathit{\mu }]\end{array}\right\}\mathrm{d}\mathit{\mu }}}.\end{array}\end{array}$$

The first term of this equation is the contribution of the atmospheric emission to the upwelling flux, which is analogous to the downwelling flux of Eq. (4), replacing the boundaries of the integral and the argument of the exponential function. The second term is the contribution of the surface emission at temperature T_s considering a surface spectral diffuse emittance ϵ [unit: dimensionless] at wavenumber ν (omitted in the notation). The third term represents the contribution of the downwelling flux that is reflected on the surface and consequently contributes to the upwelling flux. The presence of a reflected contribution to the upwelling flux is the consequence of the fact that the surface does not behave exactly as a black body. The flux emerging from the surface must take the spectral emittance of the surface into account. Under the assumption of specular reflection we have $I^{\uparrow }({\mathit{\sigma }}_{\text{BOA}},\mathit{\mu })=I^{\downarrow }({\mathit{\sigma }}_{\text{BOA}},\mathit{\mu })$.

For the atmospheric component at level k, denoted $F_{\text{atmo}}^{\uparrow }\left({\mathit{\sigma }}_k\right)$, we follow the same rationale as we previously have for the downwelling flux. However, we have to rewrite Eq. (10), replacing B[T(σ_i)] by B[T(σ_i−1)] and $({\mathit{\sigma }}_i-{\mathit{\sigma }}^{\prime })$ by $({\mathit{\sigma }}^{\prime }-{\mathit{\sigma }}_{i-\mathrm{1}})$, to account for the change in direction of propagation. And, we also have to integrate the radiance from the BOA to level σ_k instead of integrating from the TOA to level σ_k. We have $F_{\text{atmo}}^{\uparrow }\left({\mathit{\sigma }}_N\right)=\mathrm{0}$, and for every level k<N we finally obtain

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{F}_{\text{atmo}}^{\uparrow }\left({\mathit{\sigma }}_k\right)=& \mathrm{2}\mathit{\pi }\sum _{i=k+\mathrm{1}}^{N}B\left[T\right({\mathit{\sigma }}_{i-\mathrm{1}}\left)\right]\left[E_{\mathrm{3}}\right({\mathrm{\Delta }}_{i-\mathrm{1},k})-E_{\mathrm{3}}({\mathrm{\Delta }}_{i,k}\left)\right]\\ & +\mathrm{2}\mathit{\pi }\sum _{i=k+\mathrm{1}}^N\mathrm{2}\left(B\left[{\stackrel{\mathrm{‾}}{T}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}\right]-B\left[T\right({\mathit{\sigma }}_{i-\mathrm{1}}\left)\right]\right)\\ & \times \left[{\displaystyle \frac{E_{\mathrm{4}}\left({\mathrm{\Delta }}_{i-\mathrm{1},k}\right)-E_{\mathrm{4}}\left({\mathrm{\Delta }}_{i,k}\right)}{{\mathrm{\Delta }}_{i,i-\mathrm{1}}}}-E_{\mathrm{3}}\left({\mathrm{\Delta }}_{i,k}\right)\right].\end{array}\end{array}$$

The surface component at level k, denoted $F_{\text{surf}}^{\uparrow }\left({\mathit{\sigma }}_k\right)$, can be rewritten using the E₃ function as

$$\begin{array}{}\text{(14)}& F_{\text{surf}}^{\uparrow }\left({\mathit{\sigma }}_k\right)=\mathrm{2}\mathit{\pi }\mathit{ϵ}B\left(T_{\mathrm{s}}\right)E_{\mathrm{3}}\left({\mathrm{\Delta }}_{N,k}\right).\end{array}$$

Table 1 CPU time of the complete radiative transfer including calculation of radiative quantities for a single atmospheric situation.

Computations are performed sequentially on a quad-pro AMD Opteron 6378 16-core with 2.4 GHz and 256 GB of RAM. TIGR situations are described in Sect. 4.1, and RFMIP situations are described in Sect. 3.

For the reflected component at level k, denoted $F_{\text{refl}}^{\uparrow }\left({\mathit{\sigma }}_k\right)$, the integration of the third term of Eq. (12) leads to

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{F}_{\text{refl}}^{\uparrow }\left({\mathit{\sigma }}_k\right)=& \mathrm{2}\mathit{\pi }(\mathrm{1}-\mathit{ϵ})\sum _{i=\mathrm{1}}^{N}B\left[T\right({\mathit{\sigma }}_i\left)\right]\left[E_{\mathrm{3}}\right({\mathrm{\Delta }}_{N,k}+{\mathrm{\Delta }}_{N,i})\\ & -E_{\mathrm{3}}({\mathrm{\Delta }}_{N,k}+{\mathrm{\Delta }}_{N,i-\mathrm{1}})]\\ & +\mathrm{2}\mathit{\pi }(\mathrm{1}-\mathit{ϵ})\sum _{i=\mathrm{1}}^N\mathrm{2}\left(B\left[{\stackrel{\mathrm{‾}}{T}}_{i-\frac{\mathrm{1}}{\mathrm{2}}}\right]-B\left[T\right({\mathit{\sigma }}_i\left)\right]\right)\\ & \times [{\displaystyle \frac{E_{\mathrm{4}}({\mathrm{\Delta }}_{N,k}+{\mathrm{\Delta }}_{N,i})-E_{\mathrm{4}}({\mathrm{\Delta }}_{N,k}+{\mathrm{\Delta }}_{N,i-\mathrm{1}})}{{\mathrm{\Delta }}_{i,i-\mathrm{1}}}}\\ & -E_{\mathrm{3}}({\mathrm{\Delta }}_{N,k}+{\mathrm{\Delta }}_{N,i-\mathrm{1}})].\end{array}\end{array}$$

This equation can be understood by comparison to the equation for the downwelling flux (Eq. 11). The first difference is the (1−ϵ) factors before the sums that correspond to the fraction of downwelling radiation that is reflected by the surface. The second difference to Eq. (11) is that we sum all the layers from 1 to N because the reflected flux is a fraction of the downwelling flux that reaches the BOA. Finally, the third difference is the argument of the E_n functions. These arguments are now incremented by Δ_N,k, which corresponds to the optical depth from the surface to the level k (upwelling propagation after the reflection).

2.2.3 Radiative quantities

In 4A-Flux, the total upwelling flux is finally obtained by adding the three contributions described above:

$$\begin{array}{}\text{(16)}& F^{\uparrow }\left({\mathit{\sigma }}_k\right)=F_{\text{atmo}}^{\uparrow }\left({\mathit{\sigma }}_k\right)+F_{\text{surf}}^{\uparrow }\left({\mathit{\sigma }}_k\right)+F_{\text{refl}}^{\uparrow }\left({\mathit{\sigma }}_k\right).\end{array}$$

Once the monochromatic hemispherical fluxes are computed using Eqs. (11) and (16), 4A-Flux computes the net flux denoted F(σ_k) using the following convention:

$$\begin{array}{}\text{(17)}& F\left({\mathit{\sigma }}_k\right)=F^{\uparrow }\left({\mathit{\sigma }}_k\right)-F^{\downarrow }\left({\mathit{\sigma }}_k\right).\end{array}$$

The monochromatic OLR is simply defined as the net flux at the TOA F(σ₀).

The monochromatic vertical heating rate, denoted H (omitted spectral dependence), is defined as the divergence of the net flux and is conventionally expressed in units of Kelvin per day [unit: K d⁻¹]. It is the radiatively induced time rate of change in temperature T within a layer of the atmosphere. In the plane-parallel geometry, it is defined as

$$\begin{array}{}\text{(18)}& H={\displaystyle \frac{g}{c_p}}{\displaystyle \frac{\partial F}{\partial p}}={\displaystyle \frac{\partial T}{\partial t}},\end{array}$$

where g is the acceleration due to gravity and c_p is the specific heat of dry air at constant pressure. The heating rate is usually negative for the longwave, and thus the expression cooling rate is sometimes preferred to the expression heating rate that we are using here. This quantity represents the rate of change in temperature caused by radiative properties of the atmosphere. In 4A-Flux, the computation of the heating rate at every layer $k-\frac{\mathrm{1}}{\mathrm{2}}$ is performed from the net flux of adjacent levels k and k−1 by calculating the difference quotient:

$$\begin{array}{}\text{(19)}& H\left({\mathit{\sigma }}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\right)={\displaystyle \frac{g}{c_p}}{\displaystyle \frac{F\left({\mathit{\sigma }}_{k-\mathrm{1}}\right)-F\left({\mathit{\sigma }}_k\right)}{p\left({\mathit{\sigma }}_{k-\mathrm{1}}\right)-p\left({\mathit{\sigma }}_k\right)}}.\end{array}$$

We then convert the heating rate from SI units to the commonly used K d⁻¹ unit.

The computation time of the total radiative transfer including the calculation of the hemispherical fluxes, the net flux, and the heating rate profile is indicated in Table 1. The longwave fluxes are computed at the spectral resolution of the 4A/OP model (5 × 10⁻⁴ cm⁻¹) from 10 to 3250 cm⁻¹. The computation in the “OLR-only” option is a lot faster than the complete vertical profile computation. Indeed, to calculate the OLR the only requirement is to add up the contributions of the emission of all layers i that are transmitted to the TOA level. Therefore, the computational time is approximately proportional to the number of levels N. On top of the previous computation of linear computational complexity, the complete vertical profile requires the addition of the contributions of the emission of every layer i transmitted to every level k above layer i (below) for the upwelling (downwelling) flux. Therefore, the computational time is proportional to $N(N+\mathrm{1})/\mathrm{2}$, and the computational complexity is quadratic. The additional terms required for the complete vertical profile are not required for the OLR calculation. This can be deduced from Eqs. (11), (13), and (15) in which it can be seen that the factors $E_{\mathrm{3}}\left({\mathrm{\Delta }}_{i,j}\right)-E_{\mathrm{3}}\left({\mathrm{\Delta }}_{i^{\prime },j^{\prime }}\right)$ not only depend on the emission layers, but also on the level at which the fluxes are computed.

Cheruy et al. (1995)Scott and Chedin (1981)Buehler et al. (2018)Dudhia (2017)Edwards and Slingo (1996)Walters et al. (2019)Clough et al. (2005)

Table 2 List of the six benchmark models that have contributed to RFMIP-IRF and that are compared in this study.

IPSL: Institut Pierre Simon Laplace; UHH: Universität Hamburg; NOAA-GFDL: National Oceanic and Atmospheric Administration, Geophysical Fluid Dynamics Laboratory; MOHC: Met Office Hadley Centre; NERC: Natural Environment Research Council; AER: Research and Climate Group, Atmospheric and Environmental Research.

3 Evaluation of 4A-Flux over the RFMIP database

The Radiative Forcing Model Intercomparison Project (RFMIP⁴; protocol paper: Pincus et al., 2016) is one of the 23 intercomparison projects supported by the World Climate Research Program (WRCP) in the context of the 6th phase of the Coupled Model Intercomparison Project (CMIP6, see Eyring et al., 2016). The objective of the experiment RFMIP-IRF (RFMIP Instantaneous Radiative Forcing from greenhouse gases) is to characterize the accuracy of the parameterization of the IRF used in climate models under clear-sky and aerosol-free conditions. This characterization is performed in comparison to line-by-line radiative transfer models that are recognized for their very high accuracy. Thanks to the new implementation of the 4A-Flux module, the 4A/OP radiative transfer model contributed to RFMIP-IRF by providing longwave downwelling and upwelling fluxes. Calculations of radiative forcing by greenhouse gases with the six benchmark models that are participating in RFMIP-IRF, 4A/OP included, are listed in Table 2. These models and the radiative forcing evaluation results are presented in Pincus et al. (2020).

RFMIP provides a sample of 100 atmospheric situations from reanalysis of present-day conditions containing profiles of pressure, temperature, greenhouse gas concentration, and surface properties. When averaged using the provided weights, the fluxes estimated from individual profiles can be used to estimate the time-averaged global mean fluxes with very small sampling errors (relative to disagreements between models). Along with the 100 present-day atmospheric situations, 17 perturbations around the present day are also provided, leading to a total sample size of 1800 atmospheric situations (perturbations are listed in Appendix A). All six benchmark models have provided longwave integrated upwelling and downwelling flux profiles for the whole set of atmospheric profiles.

Participating models have usually provided several sets of spectrally integrated longwave flux profiles for slightly different sets of greenhouse gases, called forcing variants as described in Meinshausen et al. (2017). For the first forcing variant (denoted f1), models are free to use any subset of the 43 greenhouse gases specified in Meinshausen et al. (2017). For forcing variant f1, 4A/OP takes into account 16 gases out of the 43 gases included into the RFMIP database⁵. Both the second and third forcing variant (respectively denoted f2 and f3) use CO₂, CH₄, and N₂O. In addition to those three gases, forcing variant f2 also uses CFC-12 concentration and a modified concentration of CFC-11 to account for all omitted gases, whereas forcing variant f3 uses modified concentrations of CFC-11 and HFC-134a to represent all omitted gases.

Some participating models have also provided several sets of results based on slightly different model configurations, called physics variants. As described in Pincus et al. (2020), the model ARTS 2.3 only includes CO₂ line mixing for physics variant p2 and not for physics variant p1. The model HadGEM3-GC3.1 physics variant p1 uses a lower spectral resolution than physics variant p2 that corresponds to the high-resolution configuration. The physics variant p3 is the high-resolution configuration with MT_CKD 3.2 treatment of the water vapor continuum instead of the CAVIAR continuum.

Using the RFMIP database, we have evaluated the downwelling and upwelling fluxes computed with the 4A-Flux module of the 4A/OP model against calculations of other participating radiative transfer models. To perform this comparison, we have averaged downwelling and upwelling fluxes for the 1800 atmospheric situations using the provided weights to average the 100 atmospheric situations and uniformly averaging the 18 perturbations. For any given pair of model and physics variant, we have also averaged the downwelling and upwelling fluxes for every provided forcing variant. Forcing variants are averaged for clarity purposes as for every model, variability between forcing variants is smaller than the variability between models and physics variants (not shown here). For HadGEM3-GC3.1, we have not retained the physics variant p1 for the comparison as it uses a lower spectral resolution and thus would not usefully serve as a reference. Figure 2 shows the distance $\mathrm{\Delta }{F}_i^{\uparrow \downarrow }\left(\mathrm{p}\right)$ from the mean downwelling and upwelling flux profiles of each model i to the mean of all profiles for the eight retained model–physics variant pairs. This distance is computed as

$$\begin{array}{}\text{(20)}& \mathrm{\Delta }{F}_i^{\uparrow \downarrow }\left(\mathrm{p}\right)=F_i^{\uparrow \downarrow }\left(\mathrm{p}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{8}}}\sum _i{F}_i^{\uparrow \downarrow }\left(\mathrm{p}\right).\end{array}$$

Figure 2 Comparison of net flux profiles between line-by-line models participating in RFMIP-IRF. Panels (a) and  (b) respectively show the distance to the mean of the mean downwelling and upwelling flux profiles for every pair of model and physics variant. The averaging is first performed on the 100 atmospheric profiles using provided weights and then on the 18 perturbations (uniformly weighted). Finally, for every pair of model and physics variant, all available forcing variants are averaged. The ± 1 standard deviation range (±σ) is represented by the gray shaded area.

The ± 1 standard deviation range between models is also displayed in Fig. 2 (gray shaded area). We first notice the strong agreement between all compared models. The standard deviation between models is smaller than 0.61 W m⁻² for the downwelling flux and 0.75 W m⁻² for the upwelling flux. 4A/OP model outputs are included in the ± 1 standard deviation range at almost every level with the exceptions of the 993–976 hPa pressure range for the downwelling flux and of the 300–401 hPa pressure range for the upwelling flux; they nonetheless stay close to the ± 1 standard deviation range. The extremum difference for the 4A/OP model ranges from −0.24 to 0.49 W m⁻² for the downwelling flux profile and from −0.43 to 0.20 W m⁻² for the upwelling flux profile.

Table 3 Means and standard deviations of OLR differences (model minus 4A/OP) in W m⁻².

The OLR estimations are compared using the same averaging methodology. However, we have not averaged forcing variants to give more detailed results. The statistics for all models are calculated by first averaging forcing variants for every model–physics variant pair, and then we average these pairs together. Table 3 presents these statistics (model minus 4A/OP). The mean difference between 4A/OP and other models is −0.148 W m⁻² and the mean standard deviation is 0.218 W m⁻². This shows very good agreement between 4A/OP and other models. Note that 4A/OP simulates OLR values that are very close to the LBLRTM simulations, especially for forcing index f2 and f3 for which the exact same gas list has been simulated. In the forcing index f3, the mean difference between 4A/OP and LBLRTM is 0.007 W m⁻² with a difference standard deviation of no more than 0.098 W m⁻². The difference to all models is always smaller than 1 W m⁻² except for ARTS physics index p1 that does not include CO₂ line mixing. The greater observed difference in standard deviation noticeable for forcing variant f1 compared to f2 and f3 variants is due to the gases that are not yet implemented in 4A/OP and thus not simulated in f1, whereas all gases of f2 and f3 are identical between models. This is further explained in Appendix A where experiments, forcing, and physics variants are compared without any averaging.

The results of this comparison show very good agreement of 4A/OP with the five other participating models. The 4A/OP simulation yields results to LBLRTM, especially for the longwave upwelling flux and the OLR. This good very close agreement, while the two models differ completely in both numerical methods and the spectroscopic databases used (GEISA for 4A/OP, HITRAN for LBLRTM), is an important factor in the confidence that we can have in their results. Furthermore, these deviations between models are well below the sensitivities to typical perturbations of surface and atmospheric parameters that we will describe next.

4 Sensitivity of radiative fluxes and vertical heating rate to temperature and composition

In this section, we present a first application of the 4A-Flux module newly implemented in the radiative transfer code 4A/OP. Following Clough and Iacono (1995), we perform a sensitivity study of the OLR and the vertical heating rate to various atmospheric and surface parameters. This sensitivity is based on the Thermodynamic Initial Guess Retrieval (TIGR) atmospheric database that is developed at the Laboratoire de Météorologie Dynamique. We will first present the TIGR database and then describe the sensitivity study that we have performed on the OLR and on the longwave vertical heating rate.

4.1 TIGR atmospheric database

The Thermodynamic Initial Guess Retrieval (TIGR) atmospheric database is a climatological library containing 2311 atmospheric situations that have been carefully selected by statistical methods from radiosonde reports (Chédin et al., 1985; Achard, 1991; Chevallier et al., 1998). The atmospheric samples are described by their vertical profiles of temperature, water vapor, and ozone concentration on a pressure grid of 43 levels. The pressure levels range from the surface at 1013 up to 0.0026 hPa and the density of points along the vertical increases, while pressure decreases to correctly account for the upper atmospheric layers. The 2311 atmospheres are classified into five air mass types. The Tropical class contains situations in the tropics, the Midlat 1 class contains midlatitude situations, the Midlat 2 class contains cold midlatitude and summer polar situations, the Polar 1 class contains Northern Hemisphere very cold polar situations, and the Polar 2 class contains winter polar situations for both hemispheres. This classification has been performed by a hierarchical ascending classification depending on their virtual temperature profiles (Achard, 1991; Chédin et al., 1994).

Figure 3 Average atmospheric profiles of temperature (a), water vapor concentration (b), and ozone concentration (c) of the mean atmosphere of each of the five classes of TIGR. The temperature profile is presented on a linear scale, whereas the concentration profiles are presented on a logarithmic scale in units of grams of gas per gram of dry air. The upper panels present the profiles on a logarithmic pressure scale, and the lower panels present them on a linear pressure scale.

Figure 4 Spectral OLR at a resolution of 1 cm⁻¹ computed for the mean atmosphere of each of the five classes of TIGR. Spectrally integrated OLR values in the 10–2500 cm⁻¹ range are displayed in the legend.

Figure 3 represents the mean profiles of temperature, water vapor, and ozone concentration computed for each of the five classes of the TIGR database. OLR value calculated by 4A-Flux using these five averaged atmospheres are represented in Fig. 4. Except for temperature, water vapor, and ozone, other simulated variables are identical for all five atmospheres. Table 4 presents the concentration of well-mixed greenhouse gases (WMGHGs) considered as a reference for our simulations. For all the simulations we perform here, the surface spectral diffuse emittance is set at a constant value of 0.98 on the entire longwave spectrum. To be more concise, the averaged atmospheres of the five classes will be called by the name of their class (e.g., Tropical atmosphere or Tropical situation instead of mean atmosphere of the Tropical class).

Table 4 Concentration of WMGHGs and minor gases used for the simulation. The concentrations of all these gases are assumed to be vertically constant.

The Tropical temperature at the surface and in the troposphere is significantly higher than all other classes, especially both Polar situations that exhibit lower temperatures. Higher temperatures lead to a higher OLR, which is especially noticeable in the atmospheric IR window. This difference is lower where the absorption is important such as on the carbon dioxide ν₂ absorption band. At the tropopause, the temperature profile of the Tropical situation is the lowest. Above the tropopause, the temperature profiles are closer together.

Mainly present in the troposphere, water vapor is 15 times more concentrated at the surface of the Tropical atmosphere than at the surface of both Polar atmospheres. The Tropical atmosphere has the highest water vapor content compared to all other classes, especially compared to the dry atmospheres of the Polar classes. Consequently, we observe less of an impact of water vapor absorption on the OLR spectra of the drier atmospheres.

In all atmospheres, the ozone concentration is high in the stratosphere and in the mesosphere and relatively low in the troposphere. However, it is less spread and peaks higher in the Tropical class compared to other classes. Thus, the OLR of the mean Tropical atmosphere has more of an impact on the ozone absorption band than all other classes.

As they represent very different and opposite situations, we will only focus next on two specific classes: the hot and humid Tropical class and the cold and dry Polar 2 class. We will now present the results of a sensitivity study of the OLR and longwave heating rate profile to several perturbations of surface and atmospheric parameters using TIGR Tropical and Polar 2 atmospheric situations as input of the 4A/OP model.

4.2 OLR sensitivity

Here, we seek to quantify the effects of perturbations applied to atmospheric and surface parameters on the OLR spectrum and on the spectrally integrated OLR. To perform this analysis, we have calculated the OLR differences between a reference atmosphere and a modified atmosphere affected by a series of perturbations. The reference atmospheres used here are computed with the TIGR mean Tropical and mean Polar 2 atmospheric situations that are presented in Sect. 4.1. We first present the reference OLR spectra in both Tropical and Polar 2 situations, and then we present the effects of several thermodynamic and composition perturbations. We will describe the effects of an increase in surface temperature of 1 K above the reference and of a vertically uniform increase in the atmospheric temperature of 1 K. And, we will analyze the effects doubling the concentrations of the three main WMGHGs: carbon dioxide, methane, and nitrous oxide. Each perturbation is studied independently of one another and applied uniformly on the vertical profile. And then, we will examine the effects of a +1 % vertically uniform increase in water vapor concentration and ozone concentration.

Figure 5 Spectral OLR at a spectral resolution of 1 cm⁻¹ computed for the TIGR mean Tropical atmosphere (a) and TIGR mean Polar 2 atmosphere (b). Corresponding OLR values spectrally integrated between 10 and 2500 cm⁻¹ are provided in the legends. Spectral OLR sensitivity to temperature and composition change for the TIGR mean Tropical (c) and mean Polar 2 (d) atmospheres. Results are expressed as [OLR_perturb.−OLR_ref.] as a function of the wavenumber. The spectrally integrated sensitivities are specified in the legends and in Table 5.

Table 5 Spectrally integrated OLR sensitivity to the variations of thermodynamic and composition variables.

Figure 5 presents the OLR spectra at resolution 1 cm⁻¹ simulated by the 4A-Flux module in the Tropical situation in Fig. 5a and in the Polar 2 situation in Fig. 5b. Figure 5 also shows the effects of the studied perturbations on OLR as a function of the wavenumber in the Tropical situation in Fig. 5c and in the Polar 2 situation in Fig. 5d. Table 5 summarizes the spectrally integrated OLR variation caused by these perturbations.

The Tropical OLR spectrum and longwave integrated OLR (panel a) are greater than the Polar 2 OLR spectrum (panel b). This is mainly due to the substantially higher surface and atmospheric temperatures of the Tropical situation compared to the Polar 2 situation (almost 40 K difference between the surface temperatures). The well-known signatures of the main longwave active gases are also clearly visible in Fig. 5: carbon dioxide (580–750, 2200–2400 cm⁻¹), ozone (980–1080 cm⁻¹), water vapor (especially the far infrared and large band around 1600 cm⁻¹), methane, and nitrous oxide (1250–1350 cm⁻¹).

According to Eq. (14), the increase in the surface temperature is expected to lead to an increase in the surface Planck emission for all wavenumbers at the BOA. At the TOA, the OLR is expected to increase only in the atmospheric IR window. This is well verified here in panels (c) and (d) (pink curve). With the Tropical atmosphere being optically thicker (higher humidity), the OLR sensitivity is spectrally confined to the infrared atmospheric window compared to the Polar 2 atmosphere. Because the atmospheric window spectral range is relatively free from absorbers other than water vapor, a lower concentration of water vapor leads to a higher OLR sensitivity to the surface temperature in this spectral range.

The vertically uniform increase in the atmospheric temperature profile is expected to lead to an increase in the Planck emission of every layer of the atmosphere (Eq. 13), leading to an increase in the OLR on bands where there is gas absorption (assuming the absorption lines of the gases are marginally impacted by the temperature perturbation). This is also verified in Fig. 5. In the Tropical case, the increase extends to the whole longwave spectrum. However, for the dry Polar 2 atmosphere, there is very low absorption and emission on the infrared window.

The increase in the WMGHG concentration leads to a decrease in the surface contribution to the OLR (Eq. 14) as the total surface-to-space optical depth increases. For the contribution of the atmospheric layers to the OLR, however (Eq. 13), there is competition between the increase in the thermal emission of atmospheric layers (increasing the OLR) and the increase in absorption by gases (decreasing the OLR). The ν₂ absorption band of CO₂ shows a remarkable pattern. The increase in CO₂ causes an OLR decrease in the band wings, where the absorption dominates; however, it causes an OLR increase in the band center, where the emission dominates. Integrated over the entire longwave, the absorption dominates, causing an OLR decrease in the longwave. Jeevanjee et al. (2021) clearly explain that in the absence of H₂O, the CO₂ forcing can be considered a swap of surface emission in band wings for the stratospheric emission in the band center and that, in the presence of H₂O, the surface emission in band wings is being replaced by the emission of a colder atmosphere. In Fig. 5, in the tropical atmosphere (high RH) we observe on the CO₂ band wings a high sensitivity shared between H₂O and CO₂, whereas in the Polar 2 atmosphere (low RH) the high sensitivity is shared between the surface temperature and CO₂. The sensitivity of N₂O (550–640, 1140–1320 cm⁻¹) is dominated by the absorption, and we do not observe any local increase in OLR. The CH₄ sensitivity (1200–1370 cm⁻¹) is also dominated by the absorption.

Unlike for the WMGHGs, water vapor absorbs on the entire longwave region. The sensitivity to water vapor concentration is especially important in the Tropical situation and less important in the drier Polar 2 situation. In the Tropical case, the OLR decreases almost everywhere on the longwave spectrum, except where the ν₂ carbon dioxide band dominates the absorption. In the Polar 2 case, the sensitivity to water vapor almost completely disappears in the atmospheric window. A reduction in water vapor sensitivity is also notable on spectral bands where other gases absorb and compete with water vapor absorption (O₃ at 1000–1080 cm⁻¹ as well as both CH₄ and N₂O at 1250–1350 cm⁻¹). With water vapor being more concentrated in the lower troposphere, the WMGHGs play a relatively greater part at the TOA.

An increase in ozone concentration leads to an OLR reduction (600–800, 980–1080 cm⁻¹). With ozone being more concentrated in the stratosphere around 100 hPa, an increase in its concentration leads to a decrease in the surface and the tropospheric and lower stratospheric contributions to the OLR. The increase in the emission to space is also present but remains small as it can be seen in the 600–700 cm⁻¹ range wherein there is a local increase in the spectral OLR.

Figure 6 Spectral OLR variation due to a complete removal of atmospheric HFC134a (a), HCFC22 (b), and CFC113 (c) calculated from the TIGR mean Tropical atmosphere. The total longwave OLR variations are reported in the legends. Note that the y axis ranges from 0 to 2.5 × 10⁻³ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}({\mathrm{cm}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$.

Figure 6 shows the OLR sensitivity to three anthropogenic trace gases: HFC134a, HCFC22, and CFC113. Even if their individual contributions to the total OLR are small, these three gases together represent 0.19 W m⁻² for the tropical atmosphere, which is approximately the standard deviation between 4A/OP and other RFMIP benchmark models. However, as demonstrated and discussed by Pincus et al. (2020), the standard deviation between the six RFMIP benchmark models in terms of radiative forcing is less than 0.025 W m⁻² – way below the standard deviation obtained from fluxes because of the cancelations when computing the radiative forcing as a difference of hemispheric fluxes (Pincus et al., 2020; Mlynczak et al., 2016). Accounting for more minor gases in line-by-line radiative transfer models could also improve our understanding of model intercomparison discrepancies.

These applications show that the 4A-Flux module of 4A/OP can be used to study OLR sensitivity to thermodynamic variables and composition of the atmosphere. Thanks to its very fine spectral resolution (down to 5 × 10⁻⁴ cm⁻¹), 4A/OP can help evaluate the subtle impacts of minor longwave active gases such as HFC134a, HCFC22, and CFC113.

4.3 Vertical heating rate sensitivity

In this section, we focus on the sensitivity of the vertical heating rate profile caused by a vertically uniform increase of 1 % in the water vapor concentration and by a doubling in carbon dioxide concentration using the TIGR mean Tropical atmospheric situation and the TIGR Polar 2 atmospheric situation as reference atmospheres. First, we describe the longwave heating rate vertical and spectral distribution for the two reference atmospheric situations, and then we analyze the sensitivities that are calculated as the difference between the perturbed atmospheric state and the reference state: H_perturb.(p)−H_ref.(p).

Figure 7 presents the vertical heating rate as computed by 4A-Flux for the Tropical (Fig. 7a1 to c1) and Polar 2 (Fig. 7d1 to f1) atmospheric situations. Figure 7a1 and d1 show the vertical heating rate spectrum as a function of the wavenumber (horizontal axis) and the pressure level (vertical axis). Figure 7b1 and e1 show the heating rate profile spectrally integrated between 10 and 2500 cm⁻¹. In Fig. 7c1 and f1, the range of the x axis has been modified to focus on the tropospheric variations. In all panels, the vertical representation has been divided into two separate pressure scales. The lower sub-panels, wherein the pressure scale is linear, highlight the tropospheric heating rate, and the upper sub-panels, wherein the pressure scale is logarithmic, give better insight into the heating rate above the tropopause. The color scale has been chosen to emphasize the main radiative signature of the different gases at multiple scales. Cool colors represent a net cooling, while warm colors represent a net heating, except magenta, which represents extreme cooling values.

Figure 7 Top row (a1–f1): vertical heating rate computed by 4A-Flux using TIGR mean Tropical (Polar 2) atmospheric situation in the left (right) panels. Panels (a1) and (d1) present the spectral heating rate profile as a function of the wavenumber and the pressure level. The heating rate profile integrated on the 10–2500 cm⁻¹ spectral range is represented in panels (b1) and (e1), and a focus on the tropospheric variations is displayed in panels (c1) and (f1). Middle row (a2–f2): vertical heating rate sensitivity to a 1 % increase in H₂O concentration. Panels (a2) and (d2) present the spectral heating rate sensitivity profiles [H_perturb.(p)−H_ref.(p)] as a function of the wavenumber (horizontal axis) and the pressure level (vertical axis). The heating rate sensitivity profile integrated on the 10–2500 cm⁻¹ spectral range is represented in panels (b2) and (e2), and a focus on the stratospheric variations is displayed in panels (c2) and (f2). Bottom row (a3–f3): analog to the middle row for a doubling of CO₂ concentration (instead of H₂O). Panels (c3) and (f3) show a focus on the tropospheric variations.

Tropical and Polar 2 longwave integrated heating rate profiles show similarities (Fig. 7b1, c1, e1 and f1). The heating rate profiles are negative almost everywhere on the vertical and show a characteristic shape. In both cases, we notice that the profiles are almost constant in the troposphere before increasing around a certain pressure level, referred to as the kink, which is recognized to be associated with several atmospheric processes (Jeevanjee and Fueglistaler, 2020, and references therein). Above the kink, the heating rate increases and reaches its maximum value around 100 hPa where the heating rate is close to 0 K d⁻¹. Above 100 hPa, the heating rate profile decreases again until 1 hPa, then increases until 0.1 hPa to finally decrease at the TOA. We also notice differences of heating rate profiles between the two atmospheric situations. First, the values of the heating rate until the kink are not the same (around −2 K d⁻¹ in the Tropical case and −1 K d⁻¹ in the Polar 2 case). Second, the kink occurs around 250 hPa in the Tropical atmosphere and around the 400 hPa in the Polar 2 atmosphere, which approximately corroborates the T_kink = 220 K specified in Jeevanjee and Fueglistaler (2020). Above the kink, we notice a greater slope in the Tropical atmosphere than in the Polar 2 atmosphere. As demonstrated in Jeevanjee and Fueglistaler (2020), the upper-tropospheric kink in the heating rate profile originates from the distribution of absorption coefficients in the water vapor rotational band. More precisely, in the H₂O rotational band, the occurrence of more weakly absorbing wavenumbers is relatively constant, whereas the occurrence of more strongly absorbing wavenumbers sharply declines above a certain threshold corresponding to the kink observed in the heating rate profile. Being linked to the water vapor absorption, this also explains why the kink appears at different altitudes for the different scenarios.

The characteristic variations of the heating rate profile previously described can be further analyzed using the spectral and vertical distribution of absorption by longwave active gases. The representation in Fig. 7a1 and d1 give insight into these distributions. In both atmospheric situations, the spectra of heating rate profile display similarities. We observe the signatures of the greenhouse gases accounted for in the model. The water vapor effects dominate the Earth longwave cooling as it absorbs on a wide portion of the longwave spectrum. This is primarily due to the strong absorption of the pure rotational band peaking in the 100–200 cm⁻¹ range but extending from the 10–1000 cm⁻¹ range in the Tropical case. The water vapor effects are also caused by the relatively lower absorption by the rotational–vibrational band centered on 1600 cm⁻¹ and by the continuum of absorption. The radiative cooling induced by the ν₂ absorption band of carbon dioxide is dominant in the 580–750 cm⁻¹ range. Due to the very high absorption of carbon dioxide in the band center, most of the radiative cooling occurs in the stratosphere. The cooling is very small under the tropopause temperature inversion, whereby a net heating occurs. Above this tropopause heating, an intense cooling occurs in the stratosphere and mesosphere until the TOA. On the band wings of the carbon dioxide ν₂ absorption band, we observe a net cooling in the whole profile except at the tropopause in the Tropical case in which we observe a net heating. We observe the effects of the ozone absorption in its main absorption region: the 980–1080 cm⁻¹ range. A strong heating occurs in the upper troposphere, tropopause, lower stratosphere, and high in the mesosphere, with an intense cooling in the upper stratosphere and stratopause. The same pattern, with a lower amplitude, appears in the 1250–1350 cm⁻¹ range wherein methane and nitrous oxide absorb. The main differences of the heating rate profile spectrum between the Tropical and the Polar 2 atmospheres is the amplitude of cooling and the pressure range at which the cooling occurs. In the troposphere of the Tropical atmosphere, the cooling is more intense due to the association of a steeper temperature gradient and a higher water vapor concentration. We also notice that, due to the more vertically spread ozone concentration profile in the Polar 2 atmosphere compared to the Tropical atmosphere, the ozone cooling pressure range is wider.

The vertical heating rate sensitivity to a 1 % increase in water vapor concentration is presented in Fig. 7a2, b2, and c2 for the Tropical atmosphere and Fig. 7d2, e2, and f2 for the Polar 2 atmosphere. For the Tropical atmosphere, the spectrally integrated heating rate sensitivity profile shows a decrease in the heating rate at almost every level except for the first levels near the surface. The decrease is especially important in the 800–950 hPa range. Figure 7a2 informs us about the spectral distribution of this decrease. On both the pure rotational and rotational–vibrational transition bands, where water vapor absorption is important, we observe a heating rate increase at the level of minimum heating rate (maximum cooling rate) and a heating rate decrease of the same magnitude above the level of minimum heating rate. This behavior corresponds to an elevation of the altitude of the minimum heating rate. The decrease in heating rate in the ozone absorption band (980–1080 cm⁻¹ range) in the upper troposphere and the stratosphere illustrates the competition between water vapor absorption and ozone absorption in this region of the longwave spectrum. Compared to the Tropical atmosphere, the impact of a 1 % increase in water vapor concentration on the relatively dry Polar 2 atmosphere is less important. As represented in Fig. 7e2 and f2, the vertical heating rate sensitivity is negative for all levels. In the troposphere, the decrease is limited to −0.07 K d⁻¹. A noticeable decrease reaching −0.2 K d⁻¹ is present in the upper atmosphere with the same amplitude as in the Tropical atmosphere. Visible in Fig. 7d2, this decrease is caused by the combination of contributions of the pure rotational and rotational–vibrational transition bands of water vapor adding up at levels between 0.1 and 1 hPa, while they compensate at higher levels. Here, we also notice the elevation of the altitude of the minimum heating rate in the troposphere.

Figure 7a3 to f3 present the vertical heating rate sensitivity to a doubling in carbon dioxide concentration. As can be seen in Fig. 7b3 and e3, the heating rate profile is slightly affected in the troposphere and highly affected in the stratosphere. In the Tropical case, we notice a small increase in the heating rate in the lower troposphere until 400 hPa. Above this level, the heating rate stays constant until the tropopause. In the Polar 2 atmosphere, we observe the same tropospheric effect with a diminished amplitude. The heating rate profile stops increasing at 400 hPa. In the stratosphere and above, Tropical and Polar 2 atmospheres show a sensitivity similar to the doubling of CO₂ concentration. In this pressure range, we observe an important decrease in the heating rate. At the TOA, the cooling to space is more intense in the Tropical atmosphere than in the Polar atmosphere.

The analysis of the sensitivity spectrum shows that the main impact of doubling the CO₂ concentration on the heating rate is located in the ν₂ absorption band. In band wings, we notice an increase in the heating rate in the lower troposphere and a small decrease in the upper troposphere, followed by an important heating rate decrease in the stratosphere and above. In the band core, however, no change occurs under the tropopause. We notice a small increase in the lower stratosphere followed by a very high decrease in the upper stratosphere. In both atmospheres, we notice a small sensitivity in the stratosphere and mesosphere in the 2200–2400 cm⁻¹ range.

5 Conclusions

Taking advantage of its pre-calculated optical depth lookup table, the fast and accurate radiative transfer model 4A/OP calculates the radiance spectra for a user-defined layered atmospheric model. In an effort to enhance the capabilities offered by 4A/OP, we have developed a module called 4A-Flux that computes the hemispheric radiative flux profiles, the net flux profile, and the heating rate profile in clear-sky conditions under the assumption of local thermodynamic equilibrium, a plane-parallel atmosphere, and specular reflection on the surface. When the only required output is the OLR, it is possible to set an option that substantially decreases computation time. The angular integration is performed using the exponential integral functions E_n. The linear variation of the sub-layer Planck function for the optical depth has been implemented to better represent the emission of layers with a high optical depth.

With its 4A-Flux module, 4A/OP has contributed to RFMIP-IRF, providing the hemispheric radiative flux profiles. The hemispheric radiative flux profiles and OLR calculated by 4A/OP have been compared to the outputs of other contributing state-of-the-art radiative transfer models, showing good agreement between 4A/OP and the other models in terms of hemispheric flux profiles, with the difference to the mean of all models always being lower than 0.49 W m⁻² and almost always included in the standard deviation of all models. In terms of OLR, the mean difference between 4A/OP and other models is −0.148 W m⁻² and the mean standard deviation is 0.218 W m⁻², also showing good agreement between 4A/OP and other models. We have also shown that 4A/OP outputs are especially close to LBLRTM.

Using 4A/OP models with the 4A-Flux module, we have computed the sensitivities of the OLR and vertical heating rate to several perturbations of atmospheric and surface parameters. Our study confirms the typical OLR sensitivities to thermodynamic and composition variables. We have also seen that the increase in the water vapor concentration leads to an elevation of the altitude of the minimum monochromatic heating rate. The spectral resolution of 5 × 10⁻⁴ cm⁻¹ offered by 4A/OP allows very fine spectral sensitivity studies such as the effects of minor trace gases (e.g., HFC134a, HCFC22, and CFC113) on the OLR and vertical heating rate.

The development the 4A-Flux module and its implementation into 4A/OP radiative transfer code offer multiple possibilities for studies on the clear-sky radiative fluxes and vertical heating rate at a very high spectral resolution and on any arbitrary atmospheric and surface description. It has already contributed to the RFMIP-IRF – along with several other reference radiative transfer models – to serve as benchmark models to characterize the accuracy of the parameterization of the IRF used in climate models. The very good agreement among these reference models, while they differ in numerical methods, coding, and spectroscopic databases, is an important factor in the confidence that we can have in their results. To go further and improve our understanding of the remaining discrepancies between them, we plan to compare 4A/OP outputs using either the GEISA or HITRAN spectroscopic databases and thus be able to identify the differences in results and discriminate between those that are due to the spectroscopic properties of gases and those that are due to the radiative transfer modeling.

Appendix A: Complementary information on the evaluation of 4A-Flux over the RFMIP database

This Appendix presents further information on the differences between 4A-Flux and the other benchmark models of RFMIP-IRF. Table A1 lists the 18 experiments used in RFMIP; more detailed information is provided in the protocol paper of RFMIP (Pincus et al., 2016). Figure A1 presents the detailed comparison of global mean OLR (weighted average of the 100 atmospheric situations) between 4A/OP and the five other benchmark models of RFMIP-IRF. Here, experiments, forcing variants, and physics variants are not averaged, allowing a thorough comparison. Conditional errors that can possibly be hidden in the averaging process are clearly visible here. The results displayed in Table 3 are simply the averaged values over the experiment indices presented in this figure.

Figure A1 Means and standard deviations of OLR differences (model minus 4A/OP; W m⁻²) as a function of the experiment index. The five models compared to 4A/OP are represented by colors. The three forcing indices f1, f2, and f3 are respectively represented in individual panels (a–c). The different physics variants are differentiated with marker shapes.

Table A1 List of atmospheric conditions and gas concentrations defining the 18 experiment variants of RFMIP.

Note: for a complete presentation, please refer to the protocol paper of RFMIP (Pincus et al., 2016); PMIP stands for Paleoclimate Modelling Intercomparison Project (Kageyama et al., 2016).

As discussed in the main discussion, the 4A/OP model is especially close to the results of LBLRTM for the three forcing indices. We notice here that in terms of OLR, the three models GFDL-GRTCODE, LBLRTM, and 4A/OP are close together. The distance between 4A/OP OLR is higher for the forcing variant f1 and experiments 4 and 17 that correspond to “future” gas concentrations. The OLR computed with 4A/OP is significantly higher due to the limited gas list that has been accounted for in the simulation (16 out of the 43 specified in the RFMIP database). These minor gases have a small impact on the pre-industrial and present-day radiative forcing but play a greater part in the future radiative forcing. The concentration of some of these gases increases by orders of magnitude in the future scenario compared to the present-day situation. The difference completely disappears for forcing variants f2 and f3 for which the same ensemble of gases has been considered for all models. Gases that are neglected in the present-day forcing can play an important part in the future forcing as their concentrations increase and should thus be accounted for in such cases. If we exclude experiments 4 and 17 as well as ARTS physics variant p1 that does not include CO₂ line mixing, the difference between 4A/OP and other models never exceeds ± 1 W m⁻².