the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An improved and extended parameterization of the CO_{2} 15 µm cooling in the middle and upper atmosphere (CO2_cool_fort1.0)
Federico Fabiano
Victor Fomichev
Bernd Funke
Daniel R. Marsh
The radiative infrared cooling of CO_{2} in the middle atmosphere, where it emits under nonlocal thermodynamic equilibrium (nonLTE) conditions, is a crucial contribution to the energy balance of this region and hence to establishing its thermal structure. The nonLTE computation is too CPU timeconsuming to be fully incorporated into climate models, and hence it is parameterized. The most used parameterization of the CO_{2} 15 µm cooling for Earth's middle and upper atmosphere was developed by Fomichev et al. (1998). The valid range of this parameterization with respect to CO_{2} volume mixing ratios (VMRs) is, however, exceeded by the CO_{2} of several scenarios considered in the Coupled Climate Model Intercomparison Projects, in particular the abrupt4×CO_{2} experiment. Therefore, an extension, as well as an update, of that parameterization is both needed and timely. In this work, we present an update of that parameterization that now covers CO_{2} volume mixing ratios in the lower atmosphere from ∼0.5 to over 10 times the CO_{2} preindustrial value of 284 ppmv (i.e. 150 to 3000 ppmv). Furthermore, it is improved by using a more contemporary CO_{2} line list and the collisional rates that affect the CO_{2} cooling rates. Overall, its accuracy is improved when tested for the reference temperature profiles as well as for measured temperature fields covering all expected conditions (latitude and season) of the middle atmosphere. The errors obtained for the reference temperature profiles are below 0.5 K d^{−1} for the presentday and lower CO_{2} VMRs. Those errors increase to ∼1–2K d^{−1} at altitudes between 110 and 120 km for CO_{2} concentrations of 2 to 3 times the preindustrial values. For very high CO_{2} concentrations (4 to 10 times the preindustrial abundances), those errors are below ∼1 K d^{−1} for most regions and conditions, except at 107–135 km, where the parameterization overestimates them by ∼1.2 %. These errors are comparable to the deviation of the nonLTE cooling rates with respect to LTE at about 70 km and below, but they are negligible (several times smaller) above that altitude. When applied to a large dataset of global (pole to pole and four seasons) temperature profiles measured by MIPAS (Michelson Interferometer for Passive Atmospheric Spectroscopy) (middle and upperatmosphere mode), the errors of the parameterization for the mean cooling rate (bias) are generally below 0.5 K d^{−1}, except between $\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ and $\mathrm{3}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (∼85–98 km), where they can reach biases of 1–2 K d^{−1}. For singletemperature profiles, the cooling rate error (estimated by the root mean square – rms – of a statistically significant sample) is about 1–2 K d^{−1} below $\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ hPa (∼85 km) and above $\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (∼102 km). In the intermediate region, however, it is between 2 and 7 K d^{−1}. For elevated stratopause events, the parameterization underestimates the mean cooling rates by 3–7 K d^{−1} (∼10 %) at altitudes of 85–95 km and the individual cooling rates show a significant rms (5–15 K d^{−1}). Further, we have also tested the parameterization for the temperature obtained by a highresolution version of the Whole Atmosphere Community Climate Model (WACCMX), which shows a large temperature variability and wave structure in the middle atmosphere. In this case, the mean (bias) error of the parameterization is very small, smaller than 0.5 K d^{−1} for most atmospheric layers, reaching only maximum values of 2 K d^{−1} near $\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (∼ 96 km). The rms has values of 1–2 K d^{−1} (∼20 %) below $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{2}}$ hPa (∼80 km) and values smaller than 4 K d^{−1} (∼2 %) above 10^{−4} hPa (∼105 km). In the intermediate region between $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ and $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (85–102 km), the rms is in the range of 5–12 K d^{−1}. While these values are significant in percentage at $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$–$\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ hPa, they are very small above $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (96 km). The routine is very fast, taking (1.5–7.5) $\times {\mathrm{10}}^{\mathrm{5}}$ s, depending on the extension of the atmospheric profile, the processor and the Fortran compiler.
 Article
(11986 KB)  Fulltext XML

Supplement
(3176 KB)  BibTeX
 EndNote
Carbon dioxide is the major infrared cooler of the atmosphere from the lower stratosphere up to the lower thermosphere, where emission by nitric oxide becomes important (LópezPuertas and Taylor, 2001). However, the CO_{2} infrared emissions in the ν_{2} bands near 15 µm that are responsible for the cooling are in nonlocal thermodynamic equilibrium (nonLTE) above around 70 km (see e.g. LópezPuertas and Taylor, 2001). Thus, in addition to the difficulty of computing the cooling rates in LTE, which requires the solution of the radiative transfer equation (RTE), i.e. a nonlocal problem, we have to calculate the nonLTE populations of the emitting levels. Thus, the calculation of the nonLTE cooling rates requires the solution of the statistical equilibrium equations for all energy levels producing a significant emission and the corresponding RTE equations for all bands originating from them (see e.g. Chap. 3 in LópezPuertas and Taylor, 2001). Therefore, the exact calculation of nonLTE cooling rates in general circulation models (GCMs) or climate models that extend in height above the stratopause is virtually impractical, and hence efficient parameterizations of the CO_{2} infrared cooling have been developed for implementation in such models.
The most used parameterization of the CO_{2} 15 µm cooling for Earth's middle and upper atmosphere was developed by Fomichev et al. (1998). That parameterization is applicable for a limited range of CO_{2} abundances, up to double the preindustrial CO_{2} concentration. Nowadays, however, with the rapid increase in the CO_{2} concentration in the atmosphere and its expected increase in the coming decades, climate model projections are being carried out in much higher CO_{2} scenarios that even quadruple the preindustrial CO_{2} abundance (van Vuuren et al., 2011; O'Neill et al., 2016). For example, such scenarios are considered in the Coupled Climate Model Intercomparison Projects. Therefore, parameterizations coping with such large CO_{2} concentrations are highly in demand, and that is precisely the purpose of our work. In addition to that of Fomichev et al. (1998), other parameterizations of the CO_{2} 15 µm cooling rates have been developed in the past. In the case of Earth's atmosphere, it worth mentioning the comprehensive review of the early works reported by Fomichev et al. (1998), the summary presented in Sect. 5.8 of LópezPuertas and Taylor (2001) and the more recent work of Feofilov and Kutepov (2012). Further, just before this work was submitted, a new parameterization was developed that uses the accelerated lambda iteration and opacity distribution function techniques (Kutepov and Feofilov, 2023). For the Mars and Venus atmospheres, where CO_{2} is the most abundant species, the problem has been tackled in several studies, e.g. LópezValverde et al. (1998), Hartogh et al. (2005), LópezValverde et al. (2008) and Gilli et al. (2017, 2021). In our case, we have the option of developing a completely new parameterization to adapt other CO_{2} parameterizations (like those cited above) or to extend and improve the parameterization of Fomichev et al. (1998). Attending mainly to practical reasons of promptness, we opted for the latter.
The paper is structured as follows. A very basic description of the parameterization is presented in Sect. 2. Section 3 describes the input atmospheric parameters used in the parameterization and required for calculating the reference cooling rates. In Sect. 4 we describe the calculations of the reference LTE and nonLTE cooling rates. A detailed description of the parameterization is presented in Sect. 5. The testing and accuracy of the parameterization against (i) the reference cooling rates, (ii) the measured temperature fields of the middle atmosphere of MIPAS (Michelson Interferometer for Passive Atmospheric Spectroscopy) and (iii) the temperature profiles obtained by a highresolution version of the Whole Atmosphere Community Climate Model (WACCMX) are discussed in Sects. 6, 7 and 8, respectively. The previous cooling rate parameterization was commonly used together with a parameterization of the CO_{2} nearinfrared heating rates (Ogibalov and Fomichev, 2003) in GCMs. As we have extended the former to higher CO_{2} volume mixing ratios (VMRs) and we do not plan to extend the latter to higher CO_{2} VMRs in the near future, we assess in Sect. 9 the validity of the current nearinfrared heating parameterization for high CO_{2} VMRs. In Sect. 10, we summarize the main conclusions of the study.
As discussed above, this parameterization is essentially based on that of Fomichev et al. (1998). To compute the CO_{2} cooling rate, the atmosphere is divided into five regions: one LTE and four different nonLTE regions. The method and approximations for computing the cooling rates in those regions are described in detail in Sect. 5. The new parameterization has a finer grid and, because it has been developed to cover a larger range of CO_{2} VMRs, the boundaries of the nonLTE regions were revised and, in general, their upper boundaries were shifted to higher altitudes. The scheme consists of 83 levels in $x=\mathrm{log}(\mathrm{1000}/p(\mathrm{hPa})$), covering x=0.125 to x=20.625 spaced by 0.25. The parameterization, however, can also be used for x>20.625, where it uses the same scheme as for the NLTE4 region (see Sect. 5.4). The relationship between pressure and geometrical altitude for the reference temperature profiles is shown in Fig. S1 in the Supplement. To a first approximation, the geometric altitude z below ∼120 km is related to x by z (km) ≈7x.
The parameterization computes cooling rates for given inputs of temperature and concentrations of CO_{2}, O(^{3}P), O_{2} and N_{2} as a function of pressure. No specific grid is required, and it can be irregular. The routine interpolates the given parameters to its internal pressure grid. Possible cooling effects caused by temperature disturbances at vertical scales smaller than the internal grid of the parameterization, ≲1.75 km (see e.g. Kutepov et al., 2013), are not taken into account. That is, we assume that the cooling induced by nonresolved gravity waves propagating with a vertical wavelength of the order of or smaller than ≲1.75 km would be taken into account in the GCMs by using an appropriate GW (gravity wave) parameterization.
Further, the collisional (de)activation of CO_{2}(v_{2}) levels by the main atmospheric molecules (N_{2} and O_{2}) and by O(^{3}P) can also be prescribed. To compute the different coefficients employed by the parameterization (see Sect. 5), reference LTE and nonLTE cooling rates are required (see Sect. 4.1 and 4.2). These are calculated for selected reference atmospheres and are described in the next section.
3.1 Temperature
We used the same six pressure–temperature reference atmospheres as in Fomichev et al. (1998) for altitudes below ∼120 km. Above this altitude, they were extended up to ∼200 km with the empirical US Naval Research Laboratory Mass Spectrometer Incoherent Scatter Radar version 2.0 (MSIS2) model (Emmert et al., 2021) for medium conditions of solar activity: F_{10.7}=103 sfu (June 2011) for all atmospheres except for MLE, which was F_{10.7}=142 sfu (September 2011). These six p–T profiles cover the envelope of the climatological zonal mean temperatures of the current middle atmosphere very well, e.g. as measured by MIPAS from 2007 to 2012 (see Sect. 7). However, they do not cover the smallscale temporal and spatial temperature variability (see Fig. B6). The performance of the parameterization for such variability is addressed in Sect. 7. Further, the range of the six temperature profiles does not cover the episodes of stratospheric warming with an elevated stratopause well. During these events, the altitude region of the typical stratopause, at about 50 km, is much colder, being even 50 K colder than during normal conditions, and the altitude of the typical mesopause, near 85–90 km, is warmer by a similar amount (see Fig. 1a, profile ES). For these conditions the temperature profile is nearly isothermal from the tropopause up to about 0.1 hPa and exhibits an inversion above, with a peak near the mesopause. We anticipate that, for these rare conditions, the error incurred by this parameterization can be significant (see Sect. 7.2).
We should also mention that the envelope of these reference atmospheres does not fully cover the predicted temperatures for the end of this century for projections with high CO_{2} emissions. In particular, WACCM simulations for this century under the RCP6.0 scenario (Marsh, 2011; Marsh et al., 2013; Garcia et al., 2017) yield zonal mean temperatures which are colder in the middle atmosphere. In order to cover such predictions, the envelope of the six p–T profiles assumed here would have to be widened by about −30 K in the upper stratosphere and by about −20 K in the mesosphere. The parameterization accuracy for such predicted temperatures has not been fully assessed in this work as we have considered only the projection of high CO_{2} VMR profiles but not the corresponding predicted temperature fields. This will be the subject of future work.
3.2 CO_{2}, O(^{3}P), O_{2} and N_{2} abundances
The valid range of the parameterization of Fomichev et al. (1998) with respect to CO_{2} volume mixing ratios is exceeded by the CO_{2} concentration of several scenarios considered in the 6th Coupled Climate Model Intercomparison Project (CMIP6), in particular for the 4×CO_{2} experiment. Several CO_{2} scenarios have been proposed for the future. Thus, van Vuuren et al. (2011) proposed the scenario RCP2.6, which reaches tropospheric CO_{2} values near 1000 ppmv by the end of the century. Likewise, Meinshausen et al. (2011) suggested the high CO_{2} scenario of RCP8.5 (CMIP5), which has CO_{2} concentrations of 2000 ppmv in the second half of the 23rd century or even higher than 2000 ppmv (e.g. SSP58.5 in CMIP6; see O'Neill et al., 2016). Here we used a wide range of tropospheric CO_{2} values ranging from about half of the preindustrial (1851) value of 285 ppmv to about 10 times this value (see Fig. 1b). The specific profiles were built from a WACCM run under the CMIP6 SSP58.5 scenario (Marsh, 2011; Marsh et al., 2013; Garcia et al., 2017). Global annual mean profiles of CO_{2} were taken from WACCM simulations for years: 1851 (preindustrial), CO_{2} profile no. 2; 2014, CO_{2} profile no. 3; 2050, CO_{2} profile no. 4 ($\sim \mathrm{2}\times $preindustrial); and 2099, CO_{2} profile no. 6 ($\sim \mathrm{4}\times $preindustrial). In addition, we set up the low CO_{2} profile (no. 1) by halving preindustrial profile no. 2, the intermediate CO_{2} profile no. 5 ($\sim \mathrm{3}\times $preindustrial) from the mean of WACCM outputs for 2050 and 2099, the high CO_{2} profile no. 7 ($\sim \mathrm{5}\times $preindustrial) by multiplying WACCM output for 2099 by a factor of 1.25, and the highest CO_{2} profile no. 8 ($\sim \mathrm{10}\times $ preindustrial) by multiplying WACCM output for 2099 by a factor of 2.7. In addition to those CO_{2} VMR profiles, we also composed the intermediate profile nos. 9, 10 and 11 for testing the parameterization (see Sect. 6.2), which are shown in Fig. 1b with dashed lines. Profile nos. 9 and 11 were obtained by multiplying WACCM outputs for the years 2050 and 2099 by factors of 0.979 and 1.8, respectively. Profile no. 10 was calculated by weighting the WACCM annual mean for the years 2050 and 2099 by 0.76875 and 0.256250, respectively. WACCM provides the CO_{2} VMR profiles up to about 130 km. Above that altitude, they were calculated by using a WACCMX run for 2008, which provides a CO_{2} VMR up to near 500 km, and scaling them, in pressure levels, by the CO_{2} value of the corresponding CO_{2} profile at a pressure of $\mathrm{5}\times {\mathrm{10}}^{\mathrm{6}}$ hPa.
As discussed above, the parameterization requires the N_{2}, O_{2} and O(^{3}P) volume mixing ratio profiles for the six p–T reference atmospheres. They were taken from the MSIS2 model (Emmert et al., 2021) and are shown in Fig. S2.
We describe in this section the nonLTE cooling rates used as a reference. To compute the coefficients of the parameterization and the boundaries of the different layers, they also require the calculations of the cooling rates in LTE, which are also described in this section. Further, we have assessed the accuracy of the LTE cooling rates by comparing them with those calculated by an independent code, the Reference Forward Model (RFM, Dudhia, 2017).
4.1 Reference LTE cooling rates
The LTE cooling rates have been computed using a modified Curtis matrix formulation (Funke et al., 2012). In that computation we used as the basis for the radiative transfer calculations (e.g. the optical depths, the transmittances and their differences) the Karlsruhe Optimised and Precise Radiative Transfer Algorithm (KOPRA, Stiller et al., 2002). This code is a welltested generalpurpose linebyline radiative transfer model that includes all the known relevant processes for performing accurate radiative transfer calculations in planetary atmospheres. We used the CO_{2} line list of HITRAN 2016 (Gordon et al., 2017) and the line shapes were modelled with a Voigt profile including the pressure and temperature dependencies of the Doppler and Lorentz halfwidths. The line mixing, although of little importance in this case because the transmittances are integrated over a wide spectral range, was also taken into account (see Stiller et al., 2002). The flux transmittances were computed using a 10point Gaussian quadrature. The wavenumber grid was 0.0005 cm^{−1}. The LTE cooling rates have been computed for the CO_{2} bands associated with the vibrational states of the ν_{1}ν_{2} mode manifold covering the spectral range from 540 to 800 cm^{−1} in intervals of 10 cm^{−1}. All bands listed in HITRAN 2016 for the six most abundant isotopes in those spectral regions were included in the calculation. For reference, the accurate cooling rates computed assuming LTE conditions for the six p–T profiles and the reference eight CO_{2} VMRs are shown in Figs. S3 and S4.
In order to ensure the accuracy of these LTE cooling rates, we have compared them with those obtained with another very welltested and widely used radiative transfer code, RFM (Dudhia, 2017). This code has been used in many studies relevant to the MIPAS instrument (Fischer et al., 2008) and for the retrieval of MIPAS level2 data obtained by the University of Oxford. It worth mentioning that RFM uses a classical Curtis matrix method (doubleflux transmittance differences), while we use the modified Curtis matrix method. Figure 2 shows the results of the comparison for the US standard temperature profile and the CO_{2} VMR of Fomichev et al. (1998). Here we used a common finealtitude grid of 0.5 km. We see that the agreement between both codes is very good, with differences at most of the altitudes smaller than 0.1–0.2 K d^{−1}. Note that some of the major differences appear to be associated with small oscillations in the RFM results.
The same formulation has been used to calculate the Curtis matrices of all the CO_{2} ν_{2} bands which are required to compute the coefficients of the parameterization in the LTE region (see Sect. 5).
4.2 Reference nonLTE cooling rates
The reference linebyline nonLTE cooling rates have been computed by using the GRANADA nonLTE code. The details of the method for solving the system of equations for CO_{2} are given in Funke et al. (2012). In addition to the solution of the statistical and radiative transfer equations described in that work for the calculation of the nonLTE populations of the CO_{2} levels, here, in order to compute accurate nonLTE cooling rates and to account for the overlap between the different CO_{2} ν_{2} bands, we included an additional final iteration computing the radiation fields in all the bands by using the lambda iteration method. This algorithm shares the radiative transfer algorithm with KOPRA (Stiller et al., 2002). Thus, the details of the radiative transfer calculation related to KOPRA, e.g. line shape, spectroscopic data, wavenumber grid, or line mixing, given in LTE Sect. 4.1 above, also apply to the nonLTE calculations described here.
For this case of nonLTE cooling rates, each rovibrational band contributes according to the nonLTE populations of their upper and lower levels. The nonLTE cooling rates calculated here comprise 16 ν_{2} vibrational bands emitting and absorbing in the 15 µm region, i.e. the fundamental ν_{2} band, three first hot ν_{2} bands, seven second ν_{2} hot bands of the major CO_{2} isotopologue and ν_{2} fundamental bands of isotopologues ^{16}O^{13}C^{16}O, ^{16}O^{12}C^{18}O, ^{16}O^{12}C^{17}O, ^{16}O^{13}C^{18}O and ^{16}O^{13}C^{17}O. The contributions of other weaker ν_{2} bands arising from higher v_{2} levels, e.g. v_{2}=4, 5 or 6, are included in the calculation but have negligible contributions for the conditions of Earth's atmosphere.
Funke et al. (2012)^{a} Rate coefficient for the forward sense of the process (cm^{3} s^{−1}). ^{b} This rate is taken as 10^{−15} cm^{3} s^{−1} for temperatures lower than 150 K (see Funke et al., 2012). T is temperature (K). i and j are different CO_{2} isotopologues. i=1–6 except as noted. v_{2} denotes equivalent 2v_{1}+v_{2} states: for example, v_{2}=2 is the triad (10002, 02201, 10001).
For the calculations of the nonLTE cooling rates, a collisional scheme and collisional rates are required. Although the collisional rates affecting the CO_{2} v_{2} levels are an input parameter for the parameterization, here we have used, for the calculations of the reference cooling rates and for testing the parameterization, the collisional rates described in Funke et al. (2012). They have been recently revised and used in the nonLTE retrieval of temperature from SABER and MIPAS measurements (GarcíaComas et al., 2008, 2023). The most relevant collisional processes concerning the populations of the levels emitting in the different ν_{2} bands described above, and their rates, are listed in Table 1 for easier reference. We should note that these rates and their temperature dependencies are different from those used in the previous parameterization of Fomichev et al. (1998). The values are in general of very similar magnitudes, except for the k_{O} rate (process 1c in Table 1) that has been considered here with its upper limit, i.e. about a factor of 2 larger than in the parameterization of Fomichev et al. (1998). This rate coefficient is not well known, with uncertainties of the order of a factor of 2 (see e.g. GarcíaComas et al., 2008). While laboratory measurements are in the range of 1.5 to $\mathrm{2}\times {\mathrm{10}}^{\mathrm{12}}$ cm^{3} s^{−1}, the values derived from atmospheric observations are close to $\mathrm{6}\times {\mathrm{10}}^{\mathrm{12}}$ cm^{3} s^{−1}. Although this rate can be chosen when using this parameterization, we have optimized it for its larger value (see Table 1), as this rate has been used in the most recent nonLTE retrievals of temperature from SABER and MIPAS measurements. The effects of using half of this value on the cooling rates are discussed in Sect. 6.2. In the comparisons shown in the next sections, we consistently used the collisional rates in Table 1 for the two parameterizations.
The cooling rates near 15 µm change very little with the illumination conditions. However, those cooling rates (or more strictly speaking, the flux divergence) of the CO_{2} ν_{2} bands computed by GRANADA under daytime conditions are affected by some emission from the relaxation and/or redistribution of the solar energy absorbed in the nearinfrared bands (see e.g. LópezPuertas et al., 1990). As this absorption or heating is already taken into account by the nearinfrared (NIR) solar heating parameterization (see Sect. 9), all the nonLTE cooling rates computed here have been performed for nighttime conditions.
The results for the accurate, linebyline nonLTE cooling rates computed for the six reference p–T profiles and the eight CO_{2} VMRs are shown in Fig. 3 from the stratosphere up to the lower thermosphere and in Fig. S5 for the upper part of the parameterization, i.e. from 80 to 200 km. We observe that the altitude distribution of the cooling rates depends very much on the temperature profile. This is the major difficulty in building the parameterization. A general common feature is the maximum near the stratopause, because at these altitudes the nonLTE cooling rates do not differ significantly from those in LTE, and these are mainly driven by the high temperature of this region. Above the stratopause, the total nonLTE cooling rates depend very much on the contributions of the different bands, e.g. the ν_{2} fundamental band of the major isotopologue (FB), the contributions of the first and hot bands (Hots) and those of the ν_{2} fundamental bands of the five minor isotopologues. These contributions are shown in Fig. 4 for the contemporary CO_{2} VMR profile (no. 3) and four p–T profiles. The nonLTE cooling rates generally decrease with altitude above the stratopause, reaching a minimum near the mesopause for several p–T profiles; see e.g. the SAS and MLE atmospheres in Figs. 3 and 4. The cooling can even be negative, e.g. heating for the very cold subarctic summer (SAS) mesopause, where heating can be several Kelvin per day (see the bottomleft panel of Fig. 3). Exceptional cases are the winter atmospheres (midlatitude winter, MLW (not shown), and subarctic winter, SAW), where the mesopause is warmer and the cooling rates are high in this region. Above the mesopause, the cooling rate rapidly increases following the enhancement of the kinetic temperature. Above about 130 km, however, the cooling rates decline (see Fig. B1). At these altitudes, cooling to space is a very good approximation to the nonLTE cooling rate, which, when expressed in Kelvin per day, is proportional to the CO_{2} VMR, to the atomic oxygen density [O] and to the temperature through $\mathrm{exp}(E/kT)$ (see e.g. Sect. 9.2 in LópezPuertas and Taylor, 2001). As altitude increases, the CO_{2} VMR decreases, and so does [O]. These two effects overcome the temperature increase, leading to a net cooling rate decrease. Note the significant contribution of the hot bands in the lower thermosphere (120–150 km; see Fig. B1), essentially due to the first hot band of the major isotopologue at these altitudes, which is about 10 % of the total cooling.
The dependence of the nonLTE cooling rates on the CO_{2} abundances is illustrated in Fig. 3. We observe that, in general, the cooling rate correlates very well with the CO_{2} abundance, although that correlation is not always linear and generally depends on altitude. This is also true for the cases where we have net heating for low CO_{2} VMR, e.g. the SAS atmosphere between about 75 km and 95 km. For the MLE and MLS atmospheres, the cooling rate near ∼90 km changes from net cooling to net heating for the largest CO_{2} VMRs. Further, the very low cooling for the tropical (TRO) p–T profile near 70 km remains very low even when the CO_{2} VMR varies by a factor of 20. In the lower thermosphere, e.g. above around ∼110 km, however, the dependency of the cooling rate on the CO_{2} is very close to being linear (see Fig. S5).
A comparison of the nonLTE and LTE cooling rates for the six p–T reference atmospheres for CO_{2} VMR nos. 3 (current VMR) and 6 (4 times the preindustrial value) is shown in Figs. 5 and S6. This comparison is necessary in order to establish the boundaries of the different atmospheric regions of the parameterization (see Sect. 5). We first observe (Fig. 5) that the altitude of the departure of the cooling rate from LTE (considered the altitude where the nonLTE–LTE difference is larger than 5 %) depends on the temperature profile and ranges from $\mathrm{5.2}\times {\mathrm{10}}^{\mathrm{2}}$ hPa (∼72.5 km) for the SAS atmosphere to $\mathrm{1.2}\times {\mathrm{10}}^{\mathrm{2}}$ hPa (∼78.7 km) for the TRO atmosphere. A similar plot but for the higher CO_{2} profile no. 6 ($\sim \mathrm{4}\times $ preindustrial) is shown in Fig. S6. An overview of the altitude or pressure level of the deviation of the cooling rate from LTE is shown in Fig. B2 for the six p–T profiles and the eight CO_{2} VMR profiles. We see that the lower altitudes (higher pressures) occur for the SAS and SAW reference atmospheres. It is also evident that this altitude increases with the CO_{2} VMR, except for the SAS and SAW cases, for which it is nearly independent of the CO_{2} VMR. That is expected as, for a more abundant CO_{2} atmosphere, the 15 µm bands become optically thicker and fewer collisions are sufficient for keeping the emitting levels in LTE. Figure B2 suggests that, for higher CO_{2} VMRs, the LTE–nonLTE transition region could be placed at higher altitudes. However, as the parameterization is intended to cover the full range of CO_{2} VMR profiles, we have to be conservative, and we placed it at the lowest altitude for any p–T or CO_{2} VMR profile. Thus, it has been taken at x=9.875 ($p=\mathrm{5.14}\times {\mathrm{10}}^{\mathrm{2}}$ hPa, z≈70 km), which, except for the SAW atmosphere with the lowest CO_{2} profile, fulfils the LTE–nonLTE transition region for all the p–T and CO_{2} VMR profiles.
For completeness, Fig. B3 shows an example of the comparison of nonLTE and LTE cooling rates, including the thermosphere for the six p–T profiles and the current CO_{2} values. This shows the enormous difference between LTE and nonLTE cooling rates (nonLTE values being much smaller) in the thermosphere.
The atomic oxygen concentration is an input to the parameterization and plays a crucial role in determining the CO_{2} infrared cooling. As a consequence, it is very important in establishing the different nonLTE regions of the parameterization (Fomichev et al., 1998). To identify the atmospheric regions where it is important, we have performed a calculation by dividing the k_{O} collisional rate by a factor of 2, which is almost equivalent to changing the O(^{3}P) concentration by the same factor. Figure 6 shows this effect for four p–T profiles and considering the current CO_{2}. Generally, it is most important above around 10^{−3} hPa (∼95 km). However, for the SAS and SAW atmospheres, it is also important down to $\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ hPa (∼85 km). The fact that its importance starts being significant at different atmospheric levels for the different p–T profiles poses an additional difficulty in the development of the parameterization.
Essentially, here we follow the parameterization developed by Fomichev et al. (1998). A brief description of the method, including the most important features and equations, is given in this section.
The atmosphere is divided into five different regions (see Fig. 7) where different approaches are used for calculating the cooling rates. These regions are qualitatively the same as those defined by Fomichev et al. (1998), but their altitude extensions (except for the LTE region) have been significantly revised, mainly as a consequence of the ample range of CO_{2} abundances for which this parameterization is developed. In fact, their upper boundaries have been moved upwards (except for LTE), resulting in the following ranges: LTE x=0–9.875 (z=0–≈70 km), NLTE1 x=9.875–12.625 (z≈70–87 km), NLTE2 x=12.625–16.375 (z≈87–109 km), NLTE3 x=16.375–19.875 (z≈109–180 km) and NLTE4 x>19.875 (z≳180 km).
The lowermost (LTE) and uppermost (NLTE4) regions are the most straightforward and also the regions where the errors are in general smaller. The most difficult parts are the transition regions from LTE to nonLTE, where (i) several bands contribute to the cooling with different source functions and their relative contributions depend very much on the actual temperature profiles (see Fig. 4), and (ii) the exchange of radiation between the layers is significant and different for the considered bands. Further, although most of the radiative excitation at a given layer is produced by the absorption of photons travelling from below, the absorption of photons travelling downwards can also contribute significantly. This is the case, for example, for the stronger fundamental band near the mesopause. Further, the cooling above around 90 km also depends on the collisions with atomic oxygen. This effect can be accurately taken into account in the upper nonLTE regions where all the bands become optically thin. However, it is very difficult to represent it properly between around the mesopause and a few tens of kilometres above, where the atomic oxygen concentration varies largely and the exchange of radiation between the layers is still important.
5.1 The LTE region
The parameterization in the LTE region is based on the Curtis matrix method. The cooling rate ${\mathit{\u03f5}}_{i}^{t}\left(\mathit{\nu}\right)$ at a given pressure level x_{i}, in a spectral region ν and for a particular temperature profile t is given by
where the indices i,j refer to pressure levels x_{i} and x_{j}, and the sum is extended over pressure levels x_{j} ranging from the lower boundary, x_{j}=0, until ${x}_{{j}_{\mathrm{CM}}}=\mathrm{13.875}$. This upper boundary of the Curtis matrix has been selected to minimize the error in the lowest nonLTE region, NLTE1 (see more details in Sect. 5.2 below). ${\mathcal{A}}_{i,j}^{t}\left(\mathit{\nu}\right)$ is the modified Curtis matrix (slightly different from its usual definition; see e.g. LópezPuertas and Taylor, 2001). The factor ${\mathit{\phi}}_{j}^{t}\left(\mathit{\nu}\right)$ in Eq. (1) represents the exponential part of the Planck function and is given by
where h is the Planck constant, k_{B} is the Boltzmann constant, and ${T}_{j}^{t}$ is the temperature of the p–T profile t at level x_{j}. Similarly, $\mathit{\phi}({T}_{\mathrm{surf}}^{t},\mathit{\nu})$ corresponds to the exponential part of the Planck function for the surface temperature ${T}_{\mathrm{surf}}^{t}$. The ${v}_{i}^{t}\left(\mathit{\nu}\right)$ term accounts for the absorption at level i times the transmission from the surface up to level i at ν, so that ${v}_{i}^{t}\left(\mathit{\nu}\right)\phantom{\rule{0.125em}{0ex}}\mathit{\phi}({T}_{\mathrm{surf}}^{t},\mathit{\nu})$ accounts for the heating rate due to the absorption of the radiation from the surface (or lower boundary). The cooling rate is calculated in the spectral range from 540 to 800 cm^{−1} divided into frequency intervals, ν, 10 cm^{−1} wide. Those LTE cooling rate profiles have been calculated for each of the six p–T reference atmospheres and the eight CO_{2} VMR profiles.
In the parameterization, the Curtis matrix is expressed with an explicit temperature dependence by
where the matrix coefficients ${\mathbf{a}}_{i,j}^{t}\left(\mathit{\nu}\right)$ and ${\mathbf{b}}_{i,j}^{t}\left(\mathit{\nu}\right)$ are given by
and S_{0}(ν), S_{1}(ν) and S_{2}(ν) are the band strengths of the fundamental, first hot and second hot bands, respectively, in the ν interval. In this way, the temperature dependence, mainly caused by the band strength of the first and second hot bands, is carried out in ${\mathit{\phi}}_{i}^{t}\left(\mathit{\nu}\right)$. Those matrix coefficients are calculated for each spectral interval ν. We obtain the coefficients for the entire spectral region of the CO_{2} 15 µm bands by summing over all the ν intervals and weighting with the ν dependency of the ${\mathit{\phi}}_{i}^{t}\left(\mathit{\nu}\right)/{\mathit{\phi}}_{i}^{t}\left({\mathit{\nu}}_{\mathrm{0}}\right)$ factor, e.g.
with ν_{0}=667.3799 cm^{−1} being the frequency of the fundamental band of the major isotopologue.
Next, we define global a_{i,j} and b_{i,j} matrix coefficients, to be used for any input temperature profile, as weighted averages of ${\mathbf{a}}_{i,j}^{t}$ and ${\mathbf{b}}_{i,j}^{t}$ for the six reference p–T profiles. We introduce a set of normalized weights ${\mathit{\xi}}_{i}^{t}$, altitudedependent, for each temperature profile so that
Analogously to the matrix coefficients ${\mathbf{a}}_{i,j}^{t}\left(\mathit{\nu}\right)$ and ${\mathbf{b}}_{i,j}^{t}\left(\mathit{\nu}\right)$, we define the corresponding vector coefficients for the surface flux, ${a}_{\mathrm{surf},i}^{t}\left(\mathit{\nu}\right)$ and ${b}_{\mathrm{surf},i}^{t}\left(\mathit{\nu}\right)$, so that
with
In this way, the cooling rate ϵ_{i}, at a pressure level x_{i}, for a given input temperature profile T_{inp} is calculated in the parameterization by
The weights ${\mathit{\xi}}_{i}^{t}$ are obtained by minimizing the cost function χ(x_{i}) at each pressure level x_{i}, given by the sum of the square of the differences of the reference LTE cooling rates, ${\mathit{\u03f5}}_{\mathrm{ref}}^{t}$, and those computed by the parameterization by using Eq. (5), ${\mathit{\u03f5}}_{\mathrm{par},i}^{t}$, for each p–T profile t, e.g.
or, in more detail, by
The normalized coefficients η^{t} were originally introduced to consider different fractions of the area of Earth ascribed to each p–T reference profile. Thus, in the previous parameterization they were taken to be equal to 0.05 for the subarctic (winter and summer) profiles, 0.1 for the midlatitude (winter and summer) profiles, 0.4 for the tropical profile and 0.3 for the midlatitude equinox p–T profile. In this study, we have explored different options including the original coefficients and a uniform weighting for the six p–T profiles, and we found a smaller χ for the latter, e.g. $\mathit{\eta}=\mathrm{1}/\mathrm{6}$ for all the profiles. Hence, that was included in this version.
In that way, we have parameterized the cooling rates as a function of temperature. The cooling rates also depend on the CO_{2} VMR profiles (see Fig. 3). The parameterization incorporates the dependence on the CO_{2} abundance by calculating a_{i,j} and b_{i,j} for a generic CO_{2} profile by assuming a linear interpolation in $\mathrm{log}[{\mathbf{a}}_{i,j}/\mathrm{VMR}({x}_{i}\left)\right]$ and $\mathrm{log}[{\mathbf{b}}_{i,j}/\mathrm{VMR}({x}_{i}\left)\right]$ from the adjacent CO_{2} VMR profiles. Thus, the a_{i,j} and b_{i,j} coefficients of Eq. (3) have been calculated (and are provided) for the eight CO_{2} VMRs shown in Fig. 1b.
5.2 The NLTE1 region: the transition from LTE to nonLTE
This region is difficult to parameterize because we have several bands contributing to the cooling (see Fig. 4), and their relative contributions depend significantly on both the temperature structure and the CO_{2} VMR profile. Note that, at certain levels, the cooling induced by the weaker hot bands is larger than that of the stronger fundamental band. We should also note that the contribution of the first hot bands at high altitudes, ∼110–150 km, is not negligible (5 %–10 %, Fig. B1). This contribution is accounted for in the parameterization by implicitly assuming that it is produced by the fundamental band of the main isotopologue (see below and Sect. 4.2).
The lower boundary of this region, i.e. the LTE–NLTE1 transition, occurs, depending on the temperature profile, at altitudes from ∼70 km up to ∼85 km (0.08 to 0.004 hPa; see Fig. B2), taking place a few kilometres lower for the subarctic summer atmosphere and the lowest CO_{2} VMR. This transition region occurs at higher altitudes for larger CO_{2} VMRs. That is, the atmosphere becomes optically thicker and fewer collisions are enough to keep the levels in LTE. However, since we also need to represent the low CO_{2} VMRs, we decided to conservatively set up this region at rather low altitudes, ${x}_{b,\mathrm{1}}=\mathrm{9.875}$ (∼70 km), the same level used in the previous parameterization.
The upper limit of this region was set up in the previous parameterization at the pressure levels where collisions with O(^{3}P) start affecting the cooling rates significantly. Again, that pressure level depends on the temperature profile (and also on the O(^{3}P) concentration), being lower at ∼0.004 hPa (x≈12.4), for the subarctic summer and winter conditions (see Fig. 6). Here, we have taken the upper boundary of ${x}_{b,\mathrm{2}}=\mathrm{12.625}$ (≈87 km), slightly higher than the 12.5 value assumed in the original parameterization.
In this region we followed, as in Fomichev et al. (1998), the matrix approach discussed in Sect. 5.1 above. Thus, Eq. (5) was used but with modified ${\mathbf{a}}_{i,j}^{\prime \phantom{\rule{0.125em}{0ex}}t}$ and ${\mathbf{b}}_{i,j}^{\prime \phantom{\rule{0.125em}{0ex}}t}$ coefficients that account for the nonLTE corrections. For each p–T profile, t, we define
and
where ${\mathit{\u03f5}}_{(\mathrm{ref},\mathrm{lte})}^{t}$ and ${\mathit{\u03f5}}_{(\mathrm{ref},\mathrm{nlte})}^{t}$ are the reference LTE and nonLTE cooling rates, respectively. Then, the general ${\mathbf{a}}_{i,j}^{\prime}$ and ${\mathbf{b}}_{i,j}^{\prime}$ coefficients were calculated by following the same procedure as for the LTE region, i.e. by weighting the p–Tspecific ${\mathbf{a}}_{i,j}^{\prime \phantom{\rule{0.125em}{0ex}}t}$ and ${\mathbf{b}}_{i,j}^{\prime \phantom{\rule{0.125em}{0ex}}t}$ coefficients with a set of altitudedependent weights ${\mathit{\xi}}_{i}^{\prime \phantom{\rule{0.125em}{0ex}}t}$ and minimizing the total cost function χ(x_{i}) (see Sect. 5.1). In this way, we obtain
and
and the cooling rates are computed by using Eq. (5) but replacing a_{i,j} and b_{i,j} by ${\mathbf{a}}_{i,j}^{\prime}$ and ${\mathbf{b}}_{i,j}^{\prime}$, respectively, i.e.
This procedure, while producing a perfect match for a single atmosphere by construction, generates irregularities for other atmospheres at some levels, e.g. when using ${\mathbf{a}}_{i,j}^{\prime \phantom{\rule{0.125em}{0ex}}t}$ for atmosphere t^{′} at points where ${\mathit{\u03f5}}_{(\mathrm{ref},\mathrm{lte}),i}^{t}$ is close to zero. We observed that the irregularities were significantly mitigated by reducing the dimensions of the Curtis matrix from 83×83 to 55×55, where i=55 corresponds to x_{CM}=13.875 ($p=\mathrm{9.422}\times {\mathrm{10}}^{\mathrm{4}}$ hPa, z≈94 km), i.e. by placing x_{CM} slightly above the boundary between the NLTE1 and NLTE2 regions. Errors induced in the LTE cooling rates by the matrix reduction are negligible (smaller than 0.05 K d^{−1} at the upper boundary).
5.3 The NLTE2 and NLTE3 regions: the recurrence formula with and without correction
The parameterization in the NLTE2, NLTE3 and NLTE4 regions is based on the recurrence formula proposed by Kutepov and Fomichev (1993). This approach is valid when the cooling rate is dominated by the fundamental band, and the absorption of radiation coming from the layers above the layer at work can also be neglected (Kutepov and Fomichev, 1993; Fomichev et al., 1998). The boundaries of these regions are then adapted to the applicability of that approach. The NLTE3 boundaries were chosen to embrace the region where those conditions are fulfilled to a large degree. In the layers below, i.e. in the NLTE2 region, that formula is not accurate and requires a correction term to account for the absorption of radiation coming from the layers above and for the cooling of bands other than the fundamental one of the main isotopologue. The recurrence formula is also the basis for the calculation of the cooling rate in the NLTE4 region (see Sect. 5.4), but it is simplified because the exchange of photons within the layers of this region can be neglected. Further, we should emphasize that the dependence of the cooling rate on the CO_{2} VMR in these regions is mainly twofold: on the one hand, its direct dependence (see Eq. 7 below) and, on the other hand, the escape function which depends on the CO_{2} column above a given layer (see Figs. S7 and S8). We discuss below the boundaries of the NLTE2 and NLTE3 layers and the expressions used for the cooling rates, i.e. the recurrence formulation.
The lower boundary of the NLTE2 region is set up in the layer where the cooling rate obtained by the corrected recurrence formula is more accurate than that given by the nonLTEcorrected Curtis matrix approach (used in NLTE1). This has been set up at ${x}_{b,\mathrm{2}}=\mathrm{12.625}$ (≈87 km), which is very similar to the value in the original parameterization of 12.5 (≈85 km). The upper boundary of the NLTE2 region is set up at the pressure level where the recurrence formula does not need to be corrected to yield an accurate estimation of the cooling rate. In this work, it has been set up at ${x}_{b,\mathrm{3}}=\mathrm{16.375}$ (≈109 km), which is significantly higher than the value of 14 (≈93 km) set up in the previous parameterization. The main reason for this is that, for the higher CO_{2} VMRs used here, the atmosphere becomes optically thicker; hence, the absorption of radiation of the layer above needs to be taken into account, also at lower pressures.
Thus, the cooling rates in the NLTE2 and NLTE3 regions are calculated by
where ${\mathit{\kappa}}_{F}=\mathrm{2.55520997}\times {\mathrm{10}}^{\mathrm{11}}$ is a constant that depends on the Einstein coefficient of the fundamental band (A), on ν_{0} and on the units of ϵ (Fomichev et al., 1998)^{1}. VMR(x) is the CO_{2} VMR; M(x) is the mean molecular weight, $\mathit{\lambda}\left({x}_{i}\right)=A\phantom{\rule{0.125em}{0ex}}/\phantom{\rule{0.125em}{0ex}}[A+{l}_{t}({x}_{i}\left)\right]$, where ${l}_{t}\left({x}_{i}\right)={k}_{\mathrm{N}\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\left[{\mathrm{N}}_{\mathrm{2}}\right]+{k}_{\mathrm{O}\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\left[{\mathrm{O}}_{\mathrm{2}}\right]+{k}_{\mathrm{O}}\phantom{\rule{0.125em}{0ex}}\left[\mathrm{O}\right]$, ${k}_{{\mathrm{N}}_{\mathrm{2}}}$, ${k}_{{\mathrm{O}}_{\mathrm{2}}}$ and k_{O} are the collisional rate constants with N_{2}, O_{2} and O(^{3}P) (see Table 1); and [N_{2}], [O_{2}] and [O(^{3}P)] are the concentrations of the respective species. Note that the collisional rates depend on x_{i} through their temperature dependencies.
The $\stackrel{\mathrm{\u0303}}{\mathit{\u03f5}}$ at level x_{i}, $\stackrel{\mathrm{\u0303}}{\mathit{\u03f5}}\left({x}_{i}\right)$, is obtained by the recurrence formula
starting from the lower boundary at x_{i}=x_{b2}, where, using Eq. (7),
and ϵ(x_{b2}) is obtained by Eq. (6). The D_{i} coefficients above are given by
where
L(u) is the escape function which mainly depends on the CO_{2} column, u, above a given level x_{i}. The temperature of those layers affects this function as it influences the line shape of the CO_{2} lines and hence the probability of photons escaping into space. This is reflected in Fig. S7a, which shows L(u) as a function of the CO_{2} column for the six p–T reference atmospheres and a single CO_{2} VMR profile (no. 3, current VMR). The dependence of L(u) with the CO_{2} VMR profiles is shown in Fig. S7. In our calculations, we have used for L(u) the average of this function for the six p–T reference atmospheres.
The α(x_{i},u) parameter entering Eq. (11) and needed in the NLTE2 region has been computed by minimizing the following cost function at each point x_{i}:
After performing some sensitivity tests, we used uniform weighting for the different reference atmospheres (${\mathit{\eta}}^{t}=\mathrm{1}/\mathrm{6}$ for all atmospheres) rather than the area weighting used in the previous parameterization. Other tests were performed to determine the optimal upper boundary for the α correction: extending the region upwards reduces the error in the x=16–19 region but results in a spurious jump at the uppermost boundary, which is avoided when using a lower x_{b3} of 16.375. It is worth noting that α above x=14.5 takes values below unity, thus decreasing the escape in the region. For the fit of the optimal α, the parameterized value of ϵ(x_{b2}) is considered a starting point rather than the reference value.
In the NLTE3 region, we used the same method as in region NLTE2, except that no correction for the L(u) function is applied, i.e. $\mathit{\alpha}({x}_{i},u)=\mathrm{1}$.
5.4 The NLTE4 region
The recurrence formula described above is also valid in the uppermost NLTE4 region, but, as the CO_{2} bands are so optically thin here, the exchange of radiation within the layers of this region can be neglected, and the recurrence formula is reduced to a simpler expression, i.e. a coolingtospace term and an additional term that accounts for the absorption of the radiation emitted by the layers below its boundary. Thus, the cooling rate for this region is computed by using Eq. (7) but with a simple expression for $\stackrel{\mathrm{\u0303}}{\mathit{\u03f5}}\left({x}_{i}\right)$, $\stackrel{\mathrm{\u0303}}{\mathit{\u03f5}}\left({x}_{i}\right)=\mathrm{\Phi}\left({x}_{b\mathrm{3}}\right)\mathit{\phi}\left({x}_{i}\right)$, that gives a smooth transition to the coolingtospace approximation, e.g.
where Φ(x_{b3}) is obtained from the boundary condition
and uses the recurrence formula in Eq. (8).
The parameterization has been tested against the reference cooling rates calculated for the reference atmospheres (the six p–T profiles and the eight CO_{2} VMR profiles) (see the next section) and for intermediate CO_{2} VMRs and the k_{O} collisional rate (Sect. 6.2). Further, it has been verified for measured temperature profiles that exhibit a large variability (Sect. 7) and for the temperature profiles obtained by a highresolution version of WACCMX (Sect. 8).
6.1 Accuracy of the parameterization for the reference atmospheres
In this section, we discuss the accuracy of the current parameterization for the assumed reference atmospheres. The nonLTE models used in the original (Fomichev et al., 1998) and current parameterizations are different. Hence we expect some differences caused not just by the parameterization itself, but possibly also by the different nonLTE models. Figure 8 shows the cooling rates of this parameterization compared to those of the previous parameterization and the reference ones for a contemporary CO_{2} VMR profile (no. 3) and the six p–T profiles. The comparisons for the lower CO_{2} profiles are shown in Figs. S9 and S10 and in Figs. S11–S15 for high CO_{2} VMRs. We should clarify that, to make the comparison meaningful, the three sets of cooling rates shown here include the same updated collisional rates (Table 1). Note that these rates are different from those used in the previous parameterization (Fomichev et al., 1998). The new parameterization also supports the previous collisional rates, but it has been optimized for the new ones in Table 1. As expected, larger differences are obtained in the region between 10^{−2} hPa (∼80 km) and $\mathrm{2}\times {\mathrm{10}}^{\mathrm{5}}$ hPa (∼120 km) and are more marked for the SAS and SAW atmospheres.
The differences are more clearly illustrated in Figs. 9 and S16, where we show the mean and standard deviation of the differences for the four lowest and four highest CO_{2} VMR profiles, respectively. The improvement of the new parameterization is noticeable (compare the blue and red lines). In general, the cooling rates of the current parameterization are more accurate than in the previous one for most of the regions and temperature structures. We observe that the errors (e.g. the differences with respect to the reference nonLTE cooling rates) of the new parameterization (red curves) are very small overall. They are below ∼0.5 K d^{−1} for the current and lower CO_{2} abundances (see Fig. 9). For higher CO_{2} concentrations, between about 2 and 3 times the preindustrial values, the largest errors are ∼1–2 K d^{−1} and are located near 110–120 km (see Fig. 9 and the topleft panel in Fig. S16). The quoted values refer to the mean of the differences, although they are larger for the individual p–T atmospheres. The spread of these values is larger in the region of 10^{−2} hPa (∼80 km) to 10^{−4} hPa (∼105 km), where the root mean square (rms) reaches values between −2 and +2 K d^{−1} (Fig. 9).
For the very high CO_{2} concentrations (4, 5 and 10 times the preindustrial abundances), the errors are also very small, below ∼1 K d^{−1} for most regions and conditions, except in the 107–135 km region, where we found maximum positive biases of ∼4, ∼5 and ∼16 K d^{−1} for 4×, 5× and 10× the preindustrial CO_{2} VMR profiles (see Fig. S16). Those maximum errors are however comparable when expressed in relative terms, all about 1.2 %. The significant rms in the region of ∼80–120 km is also notable; clearly, this region is more difficult to parameterize, particularly for such a large range of CO_{2} abundances.
That increase in the differences of the new parameterization with respect to the reference calculations for the very high CO_{2} VMRs near 110 km seems to be related to the transition region from NLTE2 to NLTE3 (see Fig. 7). It looks like the cooling in the lower part of the NLTE3 region also requires correction by the α factor for high CO_{2} VMRs. This suggests that, for higher CO_{2} VMRs, the parameterization would be more accurate if this transition altitude has risen. Such a rise, however, would worsen the cooling below this boundary. This manifests the difficulty in obtaining very accurate cooling rates for a large range of CO_{2} VMRs with this method.
6.2 Assessment of the cooling rates for intermediate CO_{2} VMRs and for the k_{O} collisional rate
The aim of the parameterization is to be used for any CO_{2} VMR input profile in the range of profile nos. 1 and 8 of Fig. 1b and any plausible value for the k_{O} rates discussed in Sect. 4.2. In this section we demonstrate that the parameterization is also very accurate for CO_{2} VMRs that fall between the reference profiles used for its development and also when using different k_{O} values. In particular, we show results for the intermediate CO_{2} VMR profile nos. 9, 10 and 11 (see Fig. 1b) and the k_{O} collisional rate used in the reference calculations divided by a factor of 2.
Figure B4 shows the results of the calculation for the intermediate CO_{2} VMR profile no. 9, which is between the current CO_{2} VMR value and that projected for 2050 (2 times the preindustrial value). We can observe similar features in the calculations for the contemporary CO_{2} VMR profile (no. 3) (see Fig. 8), although the differences are slightly larger because the CO_{2} VMR profile is larger. The distinctions are more clearly seen in Fig. 10, where we show the mean and standard deviation of the differences for the six p–T profiles. The patterns in the differences, as well as their values and spreads, are very similar to those described above in Sect. 6 for the CO_{2} reference profiles. The major differences appear between 105 and 135 km, reaching maximum values of 1, 2 and 9 K d^{−1} for VMR profile nos. 9, 10 and 11, respectively. Again, we observe that the new parameterization is more accurate at practically all altitude levels. Further, the maximum values of the standard deviations of the differences for the various p–T profiles also resemble very much those discussed before, reaching maximum values of about 2 K d^{−1}, 3 K d^{−1} and 15 K d^{−1} for the respective CO_{2} VMR profiles. Note that, although these values are larger for higher CO_{2} VMRs, they are very similar when expressed in percentage.
As the CO_{2}(v_{2})–O(^{3}P) collisional rate, k_{O}, is still uncertain nowadays by about a factor of 2 (see e.g. GarcíaComas et al., 2012) and we intend this parameterization to also be used for rates different from the nominal value, we have tested its accuracy for its lowest likely value. Figure B5 shows the results of decreasing the collisional rate k_{O} by a factor of 2 for CO_{2} VMR profile no. 3 (current value) and the six p–T profiles. The errors incurred when using this rate are slightly larger than for the nominal rate. We see that, for the reduced rate, the differences are generally below 1 K d^{−1} but can have values of up to 2 K d^{−1} near 90 km for the midlatitude summer and midlatitude winter atmospheres and between ∼85 and 100 km for the SAS conditions. The improvement with respect to the previous parameterization is not that large for this case (see Fig. 11), only below 70 km and near 90 km, mainly caused by the significant difference incurred by the previous parameterization for the SAW atmosphere (see the bottomright panel in Fig. B5). The smaller differences between both versions of the parameterizations for the reduced k_{O} are likely caused because the previous one was optimized for this reduced rate.
6.3 Performance of the parameterization
In Table 2 we list some examples of the time taken to execute the parameterization for two processors, two compilers and three atmospheric intervals. Better performance is noticeable (up to a factor of 5 faster) when using the Intel Fortran Compiler Classic (ifort
) with respect to gfortran
. We did not test the ifort
compiler in the second processor, which suggests that the times obtained with ifort
and processor no. 1 could be improved by a factor of 2.7. We did not try other more modern Fortran compilers like ifx
, which could run the parameterization even faster. It is also worth mentioning that, if the atmosphere is cut in the first 50 km, the execution is about 1.76 times faster. This is because, in the lower region, where the cooling occurs under LTE conditions, the calculation involves the Curtis matrix and operations of the coefficient matrices a and b are required. Reducing this region, e.g. starting near 50 km where LTE still prevails, makes the parameterization significantly faster. Reducing the atmosphere in the upper layers, however, hardly decreases the computation time. Thus, when using a processor of type 2 with an ifort
compiler for an atmosphere in the range of 50–270 km, the calculation of a cooling rate profile could be as low as 0.015 ms and still have a margin of improvement when using modern Fortran compilers like ifx
.
7.1 Solstice and equinox conditions
We have compared the cooling rates estimated by the parameterization with the reference ones for realistic, e.g. measured, temperature profiles that present a large variability and very variable vertical structure (see e.g. Fig. B6). Specifically, we compared them for the p–T profiles measured by the MIPAS instrument (GarcíaComas et al., 2023) for 5 full days of measurements (about 2500 profiles) with global latitude coverage and covering 2 d for solstice conditions (14 January and 13 June) and 2 d for equinox conditions (25 March and 21 September) for 2010. Further, we compare the results for the temperatures of 15 February 2009 when strong stratospheric warming followed by an elevated stratopause event occurred in the northern polar hemisphere (see Sect. 7.2). The comparison is carried out for the MIPAS measurements taken only under nighttime conditions, as the MIPAS nonLTE cooling rates for daytime, obtained simultaneously with the temperature inversion, also include the fraction of the 15 µm cooling which is produced by the relaxation of the solar energy absorbed by CO_{2} NIR bands, which is not accounted for in this parameterization (see Sect. 9). The zonal means of the temperatures, CO_{2} VMRs and O(^{3}P) abundances for those conditions are shown in Figs. B7, S19 and S20, respectively.
The results are presented in Fig. 12 for the zonal mean of the differences for 1 d of solstice and 1 d of equinox conditions and in Fig. 13 as the global mean difference for all latitudes for each of the 4 individual days. In general, the new parameterization is slightly more accurate. For example, the deviations of the cooling rates from the reference calculations in the altitude range of 105–115 km are larger in the old parameterization (about 2 K d^{−1}) than in the new one and are negligible in this region. Also, the differences with respect to the reference calculations are larger in the altitude range of 80–95 km for solstice conditions and at altitudes of 80–100 km for equinox conditions (see Fig. 13).
Overall, the errors in the mean profiles of the cooling rates of the new parameterization for 1 d of measurements are below 0.5 K d^{−1}, except in the region between $\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ and $\mathrm{3}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (∼85–95 km), where they can reach values of 1–2 K d^{−1}. This region is the most difficult to parameterize because several bands contribute to the cooling rate, and they are very sensitive to the temperature structure of the middle atmosphere (e.g. even outside this region). Note also that this is precisely the region where the rms of the differences of the cooling rate with respect to the reference ones are largest, reaching values of up to 6 K d^{−1} (see Fig. 13).
7.2 Elevated stratopause conditions
The comparison of the cooling rates estimated by the old and new parameterizations with respect to the reference calculations for 15 February 2009, a day with a pronounced and unusual elevated stratopause event (see the zonal mean temperatures in Fig. B9), is shown in Fig. 14. Similar features to the other conditions shown above can be appreciated, except in the polar winter region. The mean of the differences and the standard deviations for all the profiles at latitudes north of 50° N are shown in Fig. 15. The differences are significantly larger than for other latitudes in the 80–95 km altitude region. Both parameterizations underestimate the cooling in that atmospheric region. The new parameterization has, however, better performance above about 80 km, but in the elevated stratopause region (80–100 km), it still underestimates the cooling by 3–7 K d^{−1} (∼10 %).
It seems clear that part of this underestimation is caused by the fact that such atypical temperature profiles (see Sect. 3.1) were not considered in the parameterization. However, its inclusion would not solve the problem, as in the calculations of the coefficients a tradeoff of the weighting of the different p–T reference atmospheres has to be chosen (see Sect. 5.1 and 5.2). Thus, it might ameliorate the inaccuracy for these elevated stratopause events but would worsen the accuracy for other general situations. This manifests the difficulty or limitation of this method in providing accurate nonLTE cooling rates for all the temperature structures (gradients) that we might find in the real atmosphere. Nevertheless, we have to keep in mind that these situations are sporadic and limited to high polar regions. Hence they should not impact significantly the accuracy of the cooling rates of this parameterization in global multiyear GCM simulations.
In addition to the tests above, we have also tested the parameterizations for the temperature structure obtained with a highresolution version of the WACCMX model (Liu et al., 2024). This version of WACCMX has a fine grid of 0.25° × 0.25° in latitude × longitude and a vertical resolution of 0.1 scale heights (∼0.5 km) in most of the middle and upper atmosphere, transitioning to 0.25 scale heights in the top three scale heights^{2}. With such a fine grid, the model itself can internally generate gravity waves, thus providing temperature profiles with a vertical structure very similar to that measured by highverticalresolution lidars of the mesosphere and lower thermosphere. Some examples of p–T profiles of the model exhibiting those vertical features and the latitudinal and longitudinal variabilities are shown in Fig. S23 of the Supplement. Further, the model spans from the surface up to nearly 600 km ($\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{10}}$ hPa), which is ideal for testing the parameterization. In addition to the pressure–temperature profiles of the model, their O(^{3}P), O_{2} and N_{2} VMR profiles have been used. A contemporary CO_{2} VMR (profile no. 3) was included in these calculations.
The parameterization has been tested for a total of 225 temperature profiles. They have been selected from the model output for January conditions at four latitudes, 20, 40, 60 and 70° N of the northern winter hemisphere, and two additional latitudes, 60 and 70° S of the southern summer hemisphere. For each latitude, 36 profiles corresponding to longitudes from 0 to 360° every 10° were selected. A few p–T profiles are shown in the left column of Fig. 16, and all the profiles for latitudes 20° N, 60° N, 70° N and 70° S are shown in Fig. S23 in the Supplement. We should note that those temperature structures very much resemble those measured by lidar instruments.
The results for a few representative p–T profiles are shown in Fig. 16. A few more examples are shown in Figs. S24–S26 in the Supplement. In general, the results are very similar to those obtained for the MIPAS temperatures. The parameterization works very well below $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ hPa (∼85 km), with differences generally smaller than 1–2 K d^{−1}. In the upper regions, above $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (∼105 km, not fully shown in Fig. 16 because of the small scale chosen to highlight the differences in the region of larger differences), it also works very well. In this region, the cooling rate differences are slightly larger than near $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ hPa (∼85 km), generally below 5 K d^{−1}, but they are much smaller in relative terms since the cooling rates at high altitudes are very large (of the order of 100–300 K d^{−1}). In the intermediate region, between $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ hPa and $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ hPa, the algorithm still reproduces the reference calculations rather well but not as well as in the other regions. The cooling rate differences can reach up to 10 K d^{−1} at a few isolated levels of some individual profiles which can represent up to about 20 %. It is noticeable though that, while the differences can be significant, the profile shapes of the reference and parameterized calculations are very similar (see the middle column of Fig. 16).
To have a global perspective, we have plotted in Fig. 17 the mean of the differences for all the p–T profiles together with their rms. We see that the mean (bias) of the parameterization is very small, practically below 1.5 K d^{−1} anywhere and below 0.5 K d^{−1} for most of the atmospheric layers. The rms, a representative error of individual profiles, is also small at levels below $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{2}}$ hPa (∼80 km) with values of 1–2 K d^{−1} (∼20 %) and above 10^{−4} hPa (∼105 km) with values smaller than 4 K d^{−1} (∼2 %). In the intermediate region, between $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ and $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ hPa, the rms is however significant, with most values in the range of 5–12 K d^{−1}. While these values are significant in percentage at $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$–$\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ hPa, they are very small above $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ hPa.
Some of the GCM models use the parameterization of the CO_{2} 15 µm cooling together with that of the CO_{2} NIR heating of Ogibalov and Fomichev (2003). Hence, as we are updating the former for larger CO_{2} abundances and as the update of the NIR heating parameterization for the large CO_{2} abundances is beyond the scope of this work, we investigate whether the latter is still valid for the large CO_{2} abundances. For that purpose, we compute CO_{2} NIR heating rates with the parameterization of Ogibalov and Fomichev (2003) and with the GRANADA model for the large CO_{2} concentrations for the six p–T reference atmospheres. We should note that the nonLTE models used in both parameterizations are different, and hence we expect some difference caused not just by the parameterization itself but also by the differences between the underlying nonLTE model in the NIR heating parameterization and GRANADA. The CO_{2} NIR heating rates of GRANADA were calculated with the rate coefficients and photolysis rates described in Funke et al. (2012) but updated with those described in JuradoNavarro et al. (2015, 2016) and also with those described below. In particular, the ${J}_{{\mathrm{O}}_{\mathrm{3}}}$ rate used in these calculations is ∼10 % smaller than in JuradoNavarro et al. (2015) below 100 km and thus leads to an [O(^{1}D)] of ∼10 % smaller below 90 km but that is very similar near 100 km. Above ∼100 km, the ${J}_{{\mathrm{O}}_{\mathrm{2}}}$ coefficient used in the present calculations is about 40 % smaller than in JuradoNavarro et al. (2015), leading to a similar reduction in [O(^{1}D)]. Further, we updated the following collisional rates. The rate coefficient of N_{2} + O(^{1}D) → N_{2}(1) + O has been increased by a factor of 1.08, and the collisional deactivation of N_{2}(1) with atomic oxygen (which has an important role in the heating rates; see e.g. LópezPuertas et al., 1990) has been updated from values of $\mathrm{4.5}\times {\mathrm{10}}^{\mathrm{15}}$ $(T/\mathrm{300}{)}^{\mathrm{1.5}}$ to $\mathrm{4.3}\times {\mathrm{10}}^{\mathrm{15}}$ $(T/\mathrm{300}{)}^{\mathrm{2.9}}$ cm^{3} s^{−1}.
The results of the comparison are shown in Fig. 18 for the tropical atmosphere and an intermediate solar zenith angle (SZA) of 44.5°. The region of most importance for the CO_{2} NIR heating rates is that comprised between 0.1 and 0.01 hPa (Fomichev et al., 2004). In this region, the differences between the algorithm of Ogibalov and Fomichev (2003) and GRANADA are in the range of +0.2 to −0.5 K d^{−1} for CO_{2} VMRs up to 5 times the preindustrial CO_{2} profile, e.g. about 10 % to 15 %. Hence, given that they have been computed with very different nonLTE models and the significant effect that parameters like the CO_{2} VMR above ∼90 km, the collisional rate between N_{2}(1) and O(^{3}P), the O(^{3}P) concentration itself and the rate of exchange of CO_{2} v_{3} quanta with N_{2}(1) have on this solar heating (see e.g. LópezPuertas et al., 1990), these differences are reasonable. Hence the new CO_{2} cooling rate parameterization reported here can be safely used together with the CO_{2} solar NIR heating parameterization of Ogibalov and Fomichev (2003) for CO_{2} VMRs up to 5 times the preindustrial CO_{2} profile.
An improved and extended parameterization of the CO_{2} 15 µm cooling rates of Earth's middle and upper atmosphere has been developed. It essentially follows the same method of the parameterization of Fomichev et al. (1998). The major novelty is its extended range of CO_{2} abundances, ranging from CO_{2} profiles with tropospheric values close to half of the preindustrial value to 10 times that value. This extension of CO_{2} profiles can still be safely applied to the parameterization of the CO_{2} nearinfrared heating of Ogibalov and Fomichev (2003) up to at least 5 times the preindustrial CO_{2} values and is normally combined with this cooling rate parameterization.
Other improvements or updates are as follows. They have an extended and finer vertical grid, increasing the number of levels from 8 to 83. The CO_{2} line list has been updated, from HITRAN 1992 to HITRAN 2016. Although the collisional rate coefficients affecting the CO_{2} v_{1} and v_{2} levels are input parameters for the parameterization, in this version we have used more contemporary values, e.g. as currently used in the nonLTE retrieval of temperature from CO_{2} 15 µm emissions of SABER and MIPAS measurements (GarcíaComas et al., 2008, 2023). The rate coefficients are in general of a very similar magnitude, except for the collisional deactivation of CO_{2}(v_{1},v_{2}) levels by atomic oxygen, which is now larger by approximately a factor of 2, e.g. close to its accepted upper limit. As a consequence of the larger range of CO_{2} VMR profiles, the different NLTE layers for computing the cooling rates have been significantly revised. For example, it is worth mentioning that the lowermost altitude of the coolingtospace approximation (the uppermost NLTE layer) has risen from ∼110 km to 160–170 km.
The new parameterization has been thoroughly tested against linebyline LTE and nonLTE cooling rates for (i) the six p–T reference atmospheres; (ii) the two most important input parameters (besides temperature), the CO_{2} VMR profiles and the collisional rate of CO_{2}(v_{1},v_{2}) by atomic oxygen; (iii) realistic measured temperature fields of the middle atmosphere (about 2500 profiles), including an episode of strong stratospheric warming with a very elevated stratopause; and (iv) the temperature profiles (225 profiles) obtained by a highresolution version of WACCMX capable of generating internally gravity waves and hence with temperatures showing a large variability and pronounced vertical wave structures. Further, to illustrate the improvements, the comparisons of points (i) to (iii) have also been performed for the previous parameterization.
For the reference temperature profiles, the errors of the new parameterization (mean of the differences in the cooling rates with respect to the reference calculations for the six p–T atmospheres) are below 0.5 K d^{−1} for the current and lower CO_{2} VMRs. For higher CO_{2} concentrations, between about 2 and 3 times the preindustrial values, the largest errors are ∼1–2 K d^{−1} and are located near 110–120 km. For the very high CO_{2} concentrations (from 4 to 10 times the preindustrial abundances) the errors are also very small, below ∼1 K d^{−1}, for most regions and conditions, except in the 107–135 km region, where the parameterization overestimates them in a few Kelvin per day (∼1.2 %). For these reference atmospheres, the new parameterization has a better performance for most of the atmospheric layers and temperature structures.
From the testing of the parameterization for realistic current temperature fields of the middle atmosphere as measured by MIPAS, we found that, in general, the new parameterization is slightly more accurate. In particular, in the 105–115 km range, the previous parameterization overestimates the cooling rate by 1.5 K d^{−1}, while the new one is very accurate. However, in the other height regions the difference is not so important. The new parameterization has a better performance in the 80–95 km altitude region. Overall, the errors in the mean profiles (bias) of the cooling rates of the new parameterization, calculated for four different atmospheric conditions with about 500 profiles in each of them, are below 0.5 K d^{−1}, except between $\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ and $\mathrm{3}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (∼85–95 km), where they can reach biases of 1–2 K d^{−1}. That region is the most challenging to parameterize because several CO_{2} 15 µm bands contribute to the cooling rate, and they depend very heavily on the temperature structure of the whole middle atmosphere (e.g. even outside this region). For singletemperature profiles, the cooling rate error (characterized by the rms of the difference between the reference and the parameterized cooling rates) is about 1–2 K d^{−1} below $\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ hPa (∼85 km) and above $\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (∼100 km). In the intermediate region, however, it is significant, between 2 and 7 K d^{−1}. We have further tested the parameterization against very rare and demanding situations, such as the temperature structures of stratospheric warming events with an elevated stratopause. In these situations, however, the parameterization underestimates the cooling rates by 3–7 K d^{−1} (∼10 %) at altitudes of 80–100 km, and the individual cooling rates show a significant rms (5–15 K d^{−1}).
In addition, we have tested the parameterization for the temperature structure obtained by a highresolution version of WACCMX, with the temperatures showing a large variability and pronounced vertical wave structure. The mean (bias) error of the parameterization is very small, smaller than 0.5 K d^{−1} for most atmospheric layers, and below 1.5 K d^{−1} for almost any altitude from the surface up to 200 km. The rms of the differences in the cooling rates from the parameterization and the reference model is similar to that obtained for MIPAS temperatures, with values of 1–2 K d^{−1} (∼20 %) below $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{2}}$ hPa (∼80 km) and smaller than 4 K d^{−1} (∼2 %) above 10^{−4} hPa (∼105 km). In the intermediate region, between $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ and $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ hPa, they are slightly larger than for MIPAS, with values in the range of 5–12 K d^{−1}. These values, while they are significant in relative terms at $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$–$\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ hPa, are very small in percentage above $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ hPa.
As has been shown, this parameterization has some limitations (see Sects. 6, 7.2 and 8). In order to be able to apply specific approximations for the cooling rates, it has been designed for fixed atmospheric regions where specific radiative transfer regimes prevail. Thus, its extension to a very large range of CO_{2} abundances inevitably causes a loss of accuracy for extreme cases in specific atmospheric layers. A possible solution for future updates could be to use different extensions of the nonLTE regions (i.e. Fig. 7) for different abundances of CO_{2}. Likewise, this parameterization (like the original one) was devised for use in GCMs, i.e. to produce accurate cooling rates globally, e.g. when considering all expected temperature profiles covering the different latitudinal and seasonal conditions. Thus, the ability of the parameterization to compute accurate cooling rates for individual temperature profiles with large temperature gradients in the $\mathrm{5}\times {\mathrm{10}}^{\mathrm{3}}$ hPa (∼85 km) to $\mathrm{3}\times {\mathrm{10}}^{\mathrm{4}}$ hPa (∼95 km) region is limited. On the contrary, it is extremely fast. The routine takes only 15 µs of CPU time to calculate a profile in the range of 50 to 270 km on a machine with an Intel Core i7 4.2 GHz processor when compiled with ifort
. This is more than 6600 times faster than the best option of the NLTE15μmCoolE v1.0
routine recently reported by Kutepov and Feofilov (2023). To conclude, parameterizations overcoming those limitations but retaining that speed are highly desirable for development in the future.
The routine source code is written in Fortran 90 and is available at https://doi.org/10.5281/zenodo.10849970 (LópezPuertas et al., 2024). It has been devised for implementation in general circulation models, although it can also be used for other purposes, e.g. to compute the CO_{2} 15 µm cooling rate for a given reference atmosphere.
The code is organized in a library (in the directory source/modules/) that can be included in a more complex GCM model. The subroutine to be called is CO2_NLTE_COOL
inside module file co2cool.f90
.
The following inputs are required (in order) by CO2_NLTE_COOL
.

Atmospheric profiles as a function of pressure for temperature and four VMRs of CO_{2}, O, O_{2} and N_{2}

lev0
: the index of the given pressures so that p(lev0)
is the maximum pressure level (lower boundary) to be considered for calculating the heating rate. Heating rates will be calculated from that pressure up to the minimum pressure specified in the pressure array. For example, if p is given in the range of 10^{3} to 10^{−6} hPa (or 10^{−6} to 10^{3} hPa) and p(lev0)
= 1 hPa, the heating rate will be calculated in the range of 1 to 10^{−6} hPa. 
surf_temp
: surface temperature – if set to a negative value, the temperature of the maximum pressure level will be used. 
hr
: heating rate. This is an input–output array with the same dimension of pressure. It will only be calculated at pressures in the range ofp(lev0)
(the maximum pressure considered) to the minimum specified pressure (minimum(pressure)). Note that throughout this paper we have used the term “cooling rates”, e.g. thehr
values with changed signs. 
The values are the temperature (K), pressure (hPa), VMRs (mol mol^{−1}, not ppm) and heating rate (K d^{−1}).

Input profiles can run either from the ground to the top of the atmosphere (decreasing pressures) or reverse (top to ground with increasing pressures). The pressure grid can be irregular.

Important notes:
calculations in the LTE region.
Pressure levels should include the surface pressure (near 10^{3} hPa), even if the 15 µm cooling is to be calculated only at lower pressure levels (higher altitudes), i.e.
p(lev0)
≪ 10^{3} hPa. 
If 15 µm cooling is only calculated in the nonLTE regime, it is recommended to set up the lower boundary,
p(lev0)
, close to the limit of the LTE–nonLTE transition, e.g. near 1 hPa. In this way, more timeconsuming calculations in the LTE region will be avoided.

The output is expressed as the heating rate (K d^{−1}) on the given input grid in the range of p(lev0)
to the minimum specified pressure.
To compile the routine, follow these steps.

Edit
Makefile
and change the Fortran compiler to your preferred choice (e.g.gfortran
orifort
). 
From this folder, run
make
.
The compilation produces a test program run_cool
(see below) and a module library file lib/libco2_cool.a
.
A test program, source/main.f90
, is also provided to test the parameterization on individual profiles. Its input file input.dat
has a fixed format. Do not change the number of commented lines.

First input at line 9:
n_lev
,lev0
andT_surf

Start from line 12.

Six atmospheric profiles are read (
n_lev
rows are expected).
The output heating rates are written in the output.dat
file. To test main.f90
, two input–output files are provided: input_test.dat
and input_test2.dat
with their corresponding output_test.dat
and output_test2.dat
output files. The first computes the heating in the full pressure range provided, the second only at pressures smaller than ∼1 hPa. To test the routine, follow these steps.

cp input_test.dat input.dat

./run_cool

Check that the results in
output.dat
are consistent withoutput_test.dat
. 
The same procedure can be done for test 2.
The routine is supplied with the collisional rates described in this paper (see Table 1). Nevertheless, they can be changed by the user. They are prescribed in the constants.f90
module.

The rates are defined in the form $z=a\sqrt{T}+b\phantom{\rule{0.125em}{0ex}}\mathrm{exp}(g\phantom{\rule{0.125em}{0ex}}{T}^{\mathrm{1}/\mathrm{3}})$.
The coefficients are specified as follows.

for CO_{2}–O:
a_zo, b_zo, g_zo
(default: $\mathrm{3.5}\times {\mathrm{10}}^{\mathrm{13}}$, $\mathrm{2.32}\times {\mathrm{10}}^{\mathrm{9}}$, 76.75) 
for CO_{2}–O_{2}:
a_zo2, b_zo2, g_zo2
(default: $\mathrm{7.0}\times {\mathrm{10}}^{\mathrm{17}}$, $\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{9}}$, 83.8) 
for CO_{2}–N_{2}:
a_zn2, b_zn2, g_zn2
(default: $\mathrm{7.0}\times {\mathrm{10}}^{\mathrm{17}}$, $\mathrm{6.7}\times {\mathrm{10}}^{\mathrm{10}}$, 83.8)
The most likely rate to be changed is k_{O}, probably by using smaller values. We tested a collisional k_{O} rate 2 times smaller than used in the development of the parameterization and found that its accuracy did not change significantly (see Sect. 6.2).
Although the parameterization is specifically developed for the CO_{2} 15 µm nonLTE region, it also works for the LTE region, but the user should be cautious that other important cooling rates in the LTE region, such as those of O_{3} and H_{2}O, are not included. We recommend that GCM users utilize their radiation scheme in the LTE region and this parameterization in the nonLTE region (e.g. above ∼50 or 60 km).
Boundaries of the parameterization. About the lower boundary (maximum pressure), see the notes above. About the upper boundary (minimum pressure), there is in principle no limitation, but we recommend setting it as high as the upper lid of your model. There is a large number of GC and CC models with an upper lid at $\sim {\mathrm{10}}^{\mathrm{2}}$ hPa (or ∼80 km). This parameterization can be used for such models to compute the CO_{2} nonLTE cooling rates between ∼50 and ∼80 km. We note that, under these circumstances, the cooling rates near the upper lid might not be accurate, as the contribution of the layers above the upper lid is not considered. This, however, is not a limitation of the parameterization itself but an intrinsic limitation of this kind of model. There is no restriction either on the upper limit of the upper boundary, provided it is physically meaningful. That is, it can be placed at altitudes as high as 500 km or higher.
The code is available at https://doi.org/10.5281/zenodo.10849970 (LópezPuertas et al., 2024). The parameterization is also available as a Python routine for calculating cooling rates for specific purposes at https://doi.org/10.5281/zenodo.10567258 (Fabiano et al., 2024). Note that the Python version is much slower than the Fortran version, and it is not recommended for use in GCMs.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1744012024supplement.
MLP performed the LTE and nonLTE reference calculations, participated in the adaptation of the original parameterization, wrote the manuscript and had the final editorial responsibility for this paper. FF led (together with BF) the adaptation of the original parameterization of Fomichev et al. (1998), wrote the code of the new parameterization, performed all the tests of the parameterization and calculated the cooling rates of the previous and new parameterizations presented here. VF made a critical contribution to the adaptation of its original algorithm and to its update. BF coled the adaptation of the original parameterization and designed the accuracy tests. DRM provided the CO_{2} abundance and advice on the development of the parameterization for its use and implementation in climate models. All the authors contributed to the discussions and provided text and comments.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The IAA team acknowledges financial support from the Agencia Estatal de Investigación, MCIN/AEI/10.13039/501100011033, through grant nos. PID2019110689RBI00, PID2022141216NBI00 and CEX2021001131S. We thank two anonymous referees for their very valuable suggestions leading to an improvement of this work.
This research has been supported by the Agencia Estatal de Investigación (grant nos. PID2019110689RBI00, PID2022141216NBI00 and CEX2021001131S).
The article processing charges for this openaccess publication were covered by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).
This paper was edited by Tatiana Egorova and reviewed by two anonymous referees.
Dudhia, A.: The Reference Forward Model (RFM), J. Quant. Spectrosc. Ra., 186, 243–253, https://doi.org/10.1016/j.jqsrt.2016.06.018, 2017. a, b
Emmert, J. T., Drob, D. P., Picone, J. M., Siskind, D. E., Jones, M., Mlynczak, M. G., Bernath, P. F., Chu, X., Doornbos, E., Funke, B., Goncharenko, L. P., Hervig, M. E., Schwartz, M. J., Sheese, P. E., Vargas, F., Williams, B. P., and Yuan, T.: NRLMSIS 2.0: A WholeAtmosphere Empirical Model of Temperature and Neutral Species Densities, Earth Space Sci., 8, e2020EA001321, https://doi.org/10.1029/2020EA001321, 2021. a, b
Fabiano, F., LópezPuertas, M., and Bernd, F.: CO_{2} cool – v1.1 (v1.1), Zenodo [code], https://doi.org/10.5281/zenodo.10567258, 2024. a
Feofilov, A. and Kutepov, A.: Infrared Radiation in the Mesosphere and Lower Thermosphere: Energetic Effects and Remote Sensing, Surv. Geophys., 33, 1231–1280, https://doi.org/10.1007/s1071201292040, 2012. a
Fischer, H., Birk, M., Blom, C., Carli, B., Carlotti, M., von Clarmann, T., Delbouille, L., Dudhia, A., Ehhalt, D., Endemann, M., Flaud, J. M., Gessner, R., Kleinert, A., Koopman, R., Langen, J., LópezPuertas, M., Mosner, P., Nett, H., Oelhaf, H., Perron, G., Remedios, J., Ridolfi, M., Stiller, G., and Zander, R.: MIPAS: an instrument for atmospheric and climate research, Atmos. Chem. Phys., 8, 2151–2188, https://doi.org/10.5194/acp821512008, 2008. a
Fomichev, V. I., Ogibalov, V. P., and Beagley, S. R.: Solar Heating by the NearIR CO2 Bands in the Mesosphere, Geophys. Res. Lett., 31, L21102, https://doi.org/10.1029/2004GL020324, 2004. a
Fomichev, V. L., Blanchet, J.P., and Turner, D. S.: Matrix parameterization of the 15 µm CO_{2} band cooling in the middle and upper atmosphere for variable CO_{2} concentration, J. Geophys. Res., 103, 11505–11528, 1998. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w
Funke, B., LópezPuertas, M., GarcíaComas, M., Kaufmann, M., Höpfner, M., and Stiller, G. P.: GRANADA: a Generic RAdiative traNsfer AnD nonLTE population Algorithm, J. Quant. Spectrosc. Ra., 113, 1771–1817, https://doi.org/10.1016/j.jqsrt.2012.05.001, 2012. a, b, c, d, e, f, g
Garcia, R., Smith, A., Kinnison, D., de la Camara, A., and Murphy, D.: Modification of the gravity wave parameterization in the Whole Atmosphere Community Climate Model: Motivation and results, J. Atmos. Sci., 74, 275–291, https://doi.org/10.1175/JASD160104.1, 2017. a, b
GarcíaComas, M., LópezPuertas, M., Marshall, B., Wintersteiner, P. P., Funke, B., BermejoPantaléon, D., Mertens, C. J., Remsberg, E. E., Gordley, L. L., Mlynczak, M., and Russell, J.: Errors in SABER kinetic temperature caused by nonLTE model parameters, J. Geophys. Res., 113, D24106, https://doi.org/10.1029/2008JD010105, 2008. a, b, c
GarcíaComas, M., Funke, B., LópezPuertas, M., BermejoPantaleón, D., Glatthor, N., von Clarmann, T., Stiller, G., Grabowski, U., Boone, C. D., French, W. J. R., Leblanc, T., LópezGonzález, M. J., and Schwartz, M. J.: On the quality of MIPAS kinetic temperature in the middle atmosphere, Atmos. Chem. Phys., 12, 6009–6039, https://doi.org/10.5194/acp1260092012, 2012. a
GarcíaComas, M., Funke, B., LópezPuertas, M., Glatthor, N., Grabowski, U., Kellmann, S., Kiefer, M., Linden, A., MartínezMondéjar, B., Stiller, G. P., and von Clarmann, T.: Version 8 IMK–IAA MIPAS temperatures from 12–15 µm spectra: Middle and Upper Atmosphere modes, Atmos. Meas. Tech., 16, 5357–5386, https://doi.org/10.5194/amt1653572023, 2023. a, b, c
Gilli, G., Lebonnois, S., GonzálezGalindo, F., LópezValverde, M. A., Stolzenbach, A., Lefèvre, F., Chaufray, J. Y., and Lott, F.: Thermal Structure of the Upper Atmosphere of Venus Simulated by a GroundtoThermosphere GCM, Icarus, 281, 55–72, https://doi.org/10.1016/j.icarus.2016.09.016, 2017. a
Gilli, G., Navarro, T., Lebonnois, S., Quirino, D., Silva, V., Stolzenbach, A., Lefèvre, F., and Schubert, G.: Venus Upper Atmosphere Revealed by a GCM: II. Model Validation with Temperature and Density Measurements, Icarus, 366, 114432, https://doi.org/10.1016/j.icarus.2021.114432, 2021. a
Gordon, I. E., Rothman, L. S., Hill, C., Kochanov, R. V., Tan, Y., Bernath, P. F., Birk, M., Boudon, V., Campargue, A., Chance, K. V., Drouin, B. J., Flaud, J. M., Gamache, R. R., Hodges, J. T., Jacquemart, D., Perevalov, V. I., Perrin, A., Shine, K. P., Smith, M. A. H., Tennyson, J., Toon, G. C., Tran, H., Tyuterev, V. G., Barbe, A., Császár, A. G., Devi, V. M., Furtenbacher, T., Harrison, J. J., Hartmann, J. M., Jolly, A., Johnson, T. J., Karman, T., Kleiner, I., Kyuberis, A. A., Loos, J., Lyulin, O. M., Massie, S. T., Mikhailenko, S. N., MoazzenAhmadi, N., Müller, H. S. P., Naumenko, O. V., Nikitin, A. V., Polyansky, O. L., Rey, M., Rotger, M., Sharpe, S. W., Sung, K., Starikova, E., Tashkun, S. A., Auwera, J. V., Wagner, G., Wilzewski, J., Wcisło, P., Yu, S., and Zak, E. J.: The HITRAN2016 molecular spectroscopic database, Satellite Remote Sensing and Spectroscopy: Joint ACEOdin Meeting, October 2015, 203, 3–69, 2017. a
Hartogh, P., Medvedev, A. S., Kuroda, T., Saito, R., Villanueva, G., Feofilov, A. G., Kutepov, A. A., and Berger, U.: Description and Climatology of a New General Circulation Model of the Martian Atmosphere, J. Geophys. Res.Planets, 110, E11008, https://doi.org/10.1029/2005JE002498, 2005. a
JuradoNavarro, A. A., LópezPuertas, M., Funke, B., GarcíaComas, M., Gardini, A., Stiller, G. P., and von Clarmann, T.: Vibrationvibration and vibrationthermal energy transfers of CO_{2} with N_{2} from MIPAS high resolution limb spectra, J. Geophys. Res., 120, 8002–8022, https://doi.org/10.1002/2015JD023429, 2015. a, b, c
JuradoNavarro, Á. A., LópezPuertas, M., Funke, B., GarcíaComas, M., Gardini, A., GonzálezGalindo, F., Stiller, G. P., Clarmann, T. V., Grabowski, U., and Linden, A.: Global distributions of CO2 volume mixing ratio in the middle and upper atmosphere from daytime MIPAS highresolution spectra, Atmos. Meas. Tech., 9, 6081–6100, https://doi.org/10.5194/amt960812016, 2016. a
Kutepov, A. and Feofilov, A.: New Routine NLTE15µmCoolE v1.0 for Calculating the nonLTE CO_{2} 15 µm Cooling in GCMs of Earth's atmosphere, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd2023115, in review, 2023. a, b
Kutepov, A. A. and Fomichev, V. I.: Application of the SecondOrder Escape Probability Approximation to the Solution of the NLTE VibrationRotational Band Radiative Transfer Problem, J. Atmos. Terr. Phys., 55, 1–6, 1993. a, b
Kutepov, A. A., Feofilov, A. G., Medvedev, A. S., Berger, U., Kaufmann, M., and Pauldrach, A. W. A.: InfraRed Radiative Cooling/Heating of the Mesosphere and Lower Thermosphere Due to the SmallScale Temperature Fluctuations Associated with Gravity Waves, in: Climate and Weather of the SunEarth System (CAWSES), edited by: Lübken, F.J., 429–442, Springer Netherlands, Dordrecht, ISBN 9789400743472, 9789400743489, https://doi.org/10.1007/9789400743489_23, 2013. a
Liu, H.L., Lauritzen, P. H., and Vitt, F.: Impacts of Gravity Waves on the Thermospheric Circulation and Composition, Geophys. Res. Lett., 51, e2023GL107453, https://doi.org/10.1029/2023GL107453, 2024. a
LópezPuertas, M. and Taylor, F. W.: NonLTE radiative transfer in the Atmosphere, World Scientific Pub., Singapore, 2001. a, b, c, d, e, f
LópezPuertas, M., LópezValverde, M. A., and Taylor, F. W.: Studies of Solar Heating by CO2 in the Upper Atmosphere Using a NonLTE Model and Satellite Data, J. Atmos. Sci., 47, 809–822, https://doi.org/10.1175/15200469(1990)047<0809:SOSHBC>2.0.CO;2, 1990. a, b, c
LópezPuertas, M., Fabiano, F., Fomichev, V., Funke, B., and Marsh, D. R.: CO_{2} cool (fortran version), Zenodo [code], https://doi.org/10.5281/zenodo.10849970, 2024. a, b
LópezValverde, M. A., Edwards, D. P., LópezPuertas, M., and Roldán, C.: NonLocal Thermodynamic Equilibrium in General Circulation Models of the Martian Atmosphere 1. Effects of the Local Thermodynamic Equilibrium Approximation on Thermal Cooling and Solar Heating, J. Geophys. Res., 103, 16799–16812, https://doi.org/10.1029/98JE01601, 1998. a
LópezValverde, M. A., LópezPuertas, M., and GonzálezGalindo, F.: New Parameterization of CO_{2} Cooling Rates at 15 µm for the EMGCM, ESA Rep. ESA Rep., ESA, 2008. a
Marsh, D. R.: ChemicalDynamical Coupling in the Mesosphere and Lower Thermosphere, in: Aeronomy of the Earth's Atmosphere and Ionosphere, 2, 3–17, Springer, Dordrecht, iaga special sopron book edn., 2011. a, b
Marsh, D. R., Mills, M. J., Kinnison, D. E., Lamarque, J.F., Calvo, N., and Polvani, L. M.: Climate Change from 1850 to 2005 Simulated in CESM1(WACCM), J. Climate, 26, 7372–7391, https://doi.org/10.1175/JCLID1200558.1, 2013. a, b
Meinshausen, M., Smith, S. J., Calvin, K., Daniel, J. S., Kainuma, M. L. T., Lamarque, J.F., Matsumoto, K., Montzka, S. A., Raper, S. C. B., Riahi, K., Thomson, A., Velders, G. J. M., and van Vuuren, D. P.: The RCP Greenhouse Gas Concentrations and Their Extensions from 1765 to 2300, Clim. Change, 109, 213, https://doi.org/10.1007/s105840110156z, 2011. a
Ogibalov, V. P. and Fomichev, V. I.: Parameterization of Solar Heating by the near IR CO2 Bands in the Mesosphere, Adv. Space Res., 32, 759–764, https://doi.org/10.1016/S02731177(03)800698, 2003. a, b, c, d, e, f, g
O'Neill, B. C., Tebaldi, C., van Vuuren, D. P., Eyring, V., Friedlingstein, P., Hurtt, G., Knutti, R., Kriegler, E., Lamarque, J.F., Lowe, J., Meehl, G. A., Moss, R., Riahi, K., and Sanderson, B. M.: The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6, Geosci. Model Dev., 9, 3461–3482, https://doi.org/10.5194/gmd934612016, 2016. a, b
Stiller, G. P., von Clarmann, T., Funke, B., Glatthor, N., Hase, F., Höpfner, M., and Linden, A.: Sensitivity of trace gas abundances retrievals from infrared limb emission spectra to simplifying approximations in radiative transfer modelling, J. Quant. Spectrosc. Ra., 72, 249–280, 2002. a, b, c
van Vuuren, D. P., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., Hurtt, G. C., Kram, T., Krey, V., Lamarque, J.F., Masui, T., Meinshausen, M., Nakicenovic, N., Smith, S. J., and Rose, S. K.: The Representative Concentration Pathways: An Overview, Clim. Change, 109, 5–31, https://doi.org/10.1007/s105840110148z, 2011. a, b
Note that this constant has been changed from its value of 2.63187×10^{11} in Fomichev et al. (1998) to the value used here of 2.55520997×10^{11}.
We recall that the Δx grid of the parameterization is 0.25.
 Abstract
 Introduction
 Framework of the parameterization
 The reference atmospheres
 Cooling rates for the reference atmospheres
 The parameterization
 Testing the parameterization for the reference atmospheres
 Testing the parameterization for the MIPASmeasured temperatures
 Testing the parameterization for WACCMX highresolution temperatures
 Discussion: the use of this parameterization with a previous CO_{2} solar NIR heating rate parameterization
 Summary and conclusions
 Appendix A: Notes and recommendations for using the parameterization
 Appendix B: Additional figures
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Framework of the parameterization
 The reference atmospheres
 Cooling rates for the reference atmospheres
 The parameterization
 Testing the parameterization for the reference atmospheres
 Testing the parameterization for the MIPASmeasured temperatures
 Testing the parameterization for WACCMX highresolution temperatures
 Discussion: the use of this parameterization with a previous CO_{2} solar NIR heating rate parameterization
 Summary and conclusions
 Appendix A: Notes and recommendations for using the parameterization
 Appendix B: Additional figures
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement