A mass-weighted isentropic coordinate for mapping chemical tracers and computing atmospheric inventories

We introduce a transformed isentropic coordinate $M_{{\mathit{\theta }}_{\mathrm{e}}}$, defined as the dry air mass under a given equivalent potential temperature surface (θ_e) within a hemisphere. Like θ_e, the coordinate $M_{{\mathit{\theta }}_{\mathrm{e}}}$ follows the synoptic distortions of the atmosphere but, unlike θ_e, has a nearly fixed relationship with latitude and altitude over the seasonal cycle. Calculation of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ is straightforward from meteorological fields. Using observations from the recent HIAPER Pole-to-Pole Observations (HIPPO) and Atmospheric Tomography Mission (ATom) airborne campaigns, we map the CO₂ seasonal cycle as a function of pressure and $M_{{\mathit{\theta }}_{\mathrm{e}}}$, where $M_{{\mathit{\theta }}_{\mathrm{e}}}$ is thereby effectively used as an alternative to latitude. We show that the CO₂ seasonal cycles are more constant as a function of pressure using $M_{{\mathit{\theta }}_{\mathrm{e}}}$ as the horizontal coordinate compared to latitude. Furthermore, short-term variability in CO₂ relative to the mean seasonal cycle is also smaller when the data are organized by $M_{{\mathit{\theta }}_{\mathrm{e}}}$ and pressure than when organized by latitude and pressure. We also present a method using $M_{{\mathit{\theta }}_{\mathrm{e}}}$ to compute mass-weighted averages of CO₂ on a hemispheric scale. Using this method with the same airborne data and applying corrections for limited coverage, we resolve the average CO₂ seasonal cycle in the Northern Hemisphere (mass-weighted tropospheric climatological average for 2009–2018), yielding an amplitude of 7.8 ± 0.14 ppm and a downward zero-crossing on Julian day 173 ± 6.1 (i.e., late June). $M_{{\mathit{\theta }}_{\mathrm{e}}}$ may be similarly useful for mapping the distribution and computing inventories of any long-lived chemical tracer.

1 Introduction

The spatial and temporal distribution of long-lived chemical tracers like CO₂, CH₄, and O₂ ∕ N₂ typically includes regular seasonal cycles and gradients with latitude and pressure (Conway and Tans, 1999; Ehhalt, 1978; Randerson et al., 1997; Rasmussen and Khalil, 1981; Tohjima et al., 2012). These patterns are evident in climatological averages but are potentially distorted on short timescales by synoptic weather disturbances, especially at middle to high latitudes (i.e., poleward of 30^∘ N and S) (Parazoo et al., 2008; Wang et al., 2007). With a temporally dense dataset such as from satellite remote sensing or tower in situ measurements, climatological averages can be created by averaging over this variability. For temporally sparse datasets such as from airborne campaigns, it may be necessary to correct for synoptic distortion.

A common approach to correct synoptic distortion is to use transformed coordinates rather than geographic coordinates (i.e., pressure–latitude), to take into account atmospheric dynamics and transport barriers. Such coordinate transformation has been used, for example, to reduce dynamically induced variability in the stratosphere using equivalent latitude rather than latitude as the horizontal coordinate (Butchart and Remsberg, 1986); to diagnose the tropopause profile using a tropopause-based rather than surface-based vertical coordinate (Birner et al., 2002); to study the transport regime in the Arctic using a horizontal coordinate based on the polar dome (Bozem et al., 2019); and to study UTLS (upper troposphere–lower stratosphere) tracer data by using tropopause-based, jet-based, and equivalent latitude coordinates (Petropavlovskikh et al., 2019). In the troposphere, a transformed coordinate, the isentropic coordinate (θ), has been widely applied to evaluate the distribution of tracer data (Miyazaki et al., 2008; Parazoo et al., 2011, 2012). As air parcels move with synoptic disturbances, θ and the tracer tend to be similarly displaced so that the θ–tracer relationship is relatively conserved (Keppel-Aleks et al., 2011). Furthermore, vertical mixing tends to be rapid on θ surfaces, so θ and tracer contours are often nearly parallel (Barnes et al., 2016). However, θ varies greatly with latitude and altitude over seasons due to changes in heating and cooling with solar insolation, which complicates the interpretation of θ–tracer relationships on seasonal timescales.

During analysis of airborne data from the HIAPER Pole-to-Pole Observations (HIPPO) (Wofsy, 2011) and the Atmospheric Tomography Mission (ATom) (Prather et al., 2018) airborne campaigns, we have found it useful to transform potential temperature into a mass-based unit, M_θ, which we define as the total mass of dry air under a given isentropic surface in the hemisphere. In contrast to θ, which has large seasonal variation, M_θ has a more stable relationship to latitude and altitude, while varying in parallel with θ on synoptic scales. Also, for a tracer which is well-mixed on θ, a plot of this tracer versus M_θ can be directly integrated to yield the atmospheric inventory of the tracer, because M_θ directly corresponds to the mass of air. We note that a similar concept to $M_{{\mathit{\theta }}_{\mathrm{e}}}$ has been introduced in the stratosphere by Linz et al. (2016), in which M(θ) is defined as the mass above the θ surface, to study the relationship between age of air and diabatic circulation of the stratosphere.

Several choices need to be made in the definition of M_θ, including defining boundary conditions (e.g., in altitude and latitude) for mass integration and whether to use potential temperature θ or equivalent potential temperature θ_e. Here, for boundaries, we use the dynamical tropopause (based on the potential vorticity unit, PVU) and the Equator, thus integrating the dry air mass of the troposphere in each hemisphere. We also focus on M_θ defined using equivalent potential temperature (θ_e) to conserve moist static energy in the presence of latent heating during vertical motion, which improves alignment between mass transport and mixing especially within storm tracks in mid-latitudes (Parazoo et al., 2011; Pauluis et al., 2008, 2010). We call this tracer $M_{{\mathit{\theta }}_{\mathrm{e}}}$.

In this paper we describe the method for calculating $M_{{\mathit{\theta }}_{\mathrm{e}}}$ and discuss its variability on synoptic to seasonal scales. We also discuss the time variation in the θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ relationship within each hemisphere and explore the stability of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ and the θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ relationship using different reanalysis products. To illustrate the application of $M_{{\mathit{\theta }}_{\mathrm{e}}}$, we map CO₂ data from two recent airborne campaigns (HIPPO and ATom) on $M_{{\mathit{\theta }}_{\mathrm{e}}}$. Further, we show how $M_{{\mathit{\theta }}_{\mathrm{e}}}$ can be used to accurately compute the average CO₂ concentration over the entire troposphere of the Northern Hemisphere using measurements from the same airborne campaigns. We examine the accuracy of this method and propose an appropriate way to sample the atmosphere with aircraft to compute the average of a chemical tracer within a large zonal domain.

2 Methods

2.1 Meteorological reanalysis products

The calculation of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ requires the distribution of dry air mass and θ_e. For these quantities, we alternately use three reanalysis products: ERA-Interim (Dee et al., 2011), NCEP2 (Kanamitsu et al., 2002), and Modern-Era Retrospective analysis for Research and Applications Version 2 (MERRA-2) (Gelaro et al., 2017). All products have a 2.5^∘ horizontal resolution. NCEP2 has a daily resolution, and we average 6-hourly ERA-Interim fields and 3-hourly MERRA-2 fields to yield daily fields. ERA-Interim has 32 vertical levels from 1000 to 1 mbar, with approximately 20 to 27 levels in the troposphere. NCEP2 has 17 vertical levels from 1000 to 10 mbar, with approximately 8 to 12 levels in the troposphere. MERRA-2 has 42 vertical levels from 985 to 0.01 mbar, with approximately 21 to 25 levels in the troposphere.

2.2 Equivalent potential temperature (θ_e) and dry air mass (M) of the atmospheric fields

We compute θ_e (K) using the following expression:

$$\begin{array}{}\text{(1)}& {\mathit{\theta }}_{\mathrm{e}}=\left(T+{\displaystyle \frac{L_{\mathrm{v}}\left(T\right)}{C_{\mathrm{pd}}}}\cdot w\right)\cdot {\left({\displaystyle \frac{P_{\mathrm{0}}}{P}}\right)}^{\frac{R_{\mathrm{d}}}{C_{\mathrm{pd}}}}\end{array}$$

from Stull (2012). T (K) is the temperature of air; w (kg of water vapor per kg of air mass) is the water vapor mixing ratio; R_d (287.04 J kg⁻¹ K⁻¹) is the gas constant for air; C_pd (1005.7 J kg⁻¹ K⁻¹) is the specific heat of dry air at constant pressure; P₀ (1013.25 mbar) is the reference pressure at the surface, and L_v(T) is the latent heat of evaporation at temperature T. L_v(T) is defined as 2406 kJ kg⁻¹ at 40 ^∘C and 2501 kJ kg⁻¹ at 0 ^∘C and scales linearly with temperature.

Following Bolton (1980), we compute the water vapor mixing ratio (w) from relative humidity (RH; kg kg⁻¹) provided by the reanalysis products and the formula for the saturation mixing ratio of water vapor (P_s,v; mbar) modified by Wexler (1976).

$$\begin{array}{}\text{(2)}& {\displaystyle }{P}_{\mathrm{s},\mathrm{v}}=\mathrm{0.06122}\cdot e^{\frac{\mathrm{17.67}\cdot T}{T+\mathrm{243.5}}}\text{(3)}& {\displaystyle }w=\mathrm{RH}\cdot \mathrm{0.622}\cdot {\displaystyle \frac{P_{\mathrm{s},\mathrm{v}}}{P-P_{\mathrm{s},\mathrm{v}}}}\end{array}$$

We compute the total air mass of each grid cell x at time t, M_x(t), shown in Eq. (4), from the product of the pressure range and surface area and divided by a latitude- and height-dependent gravity constant provided by Arora et al. (2011). The surface area is computed by using latitude (Φ), longitude (λ), and the radius of the Earth (R, 6371 km). The total air mass of each grid cell is computed from

$$\begin{array}{}\text{(4)}& M_x={\displaystyle \frac{\mathrm{\Delta }P}{g}}\cdot \left|\mathrm{\Delta }\sin \left(\mathrm{\Phi }\right)\cdot \mathrm{\Delta }\mathit{\lambda }\right|\cdot R^{\mathrm{2}},\end{array}$$

where Δ represents the difference between two boundaries of each grid cell.

The gravity constant (g; kg m⁻²) is computed following Arora et al. (2011) as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}g(\mathrm{\Phi },h)& =g_{\mathrm{0}}\cdot \left(\mathrm{1}+\mathrm{0.0053}\cdot {\sin }^{\mathrm{2}}\left(\mathrm{\Phi }\right)\right.\\ & \left.-\mathrm{0.000006}\cdot {\sin }^{\mathrm{2}}(\mathrm{2}\cdot \mathrm{\Phi })\right)-\mathrm{0.000003086}\cdot h\end{array},\end{array}$$

where the reference gravity constant (g₀) is assumed to be 9.78046 m s⁻² and the height (h) in units of meters is computed from

$$\begin{array}{}\text{(6)}& P=P_{\mathrm{0}}\cdot e^{-\frac{h}{H}},\end{array}$$

where H is the scale height of the atmosphere and assumed to be 8400 m.

The dry air mass is then computed by subtracting the water mass, computed from relative humidity, the saturation water vapor mass mixing ratio, and the total air mass of the grid cell (Eq. 3). Since this study focuses on tracer distributions in the troposphere, we compute $M_{{\mathit{\theta }}_{\mathrm{e}}}$ with an upper boundary at the dynamical tropopause defined as the 2 PVU (potential vorticity units, 10⁻⁶ K kg⁻¹ m² s⁻¹) surface.

ERA-Interim and NCEP2 include hypothetical levels below the true land or sea surface, for example, the 850 hPa level over the Himalaya, which we exclude in the calculation of $M_{{\mathit{\theta }}_{\mathrm{e}}}$.

2.3 Determination of $M_{{\mathit{\theta }}_{\mathrm{e}}}$

We show a schematic of the conceptual basis for the calculation of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ in Fig. 1. To compute $M_{{\mathit{\theta }}_{\mathrm{e}}}$, we sort all tropospheric grid cells in the hemisphere by increasing θ_e and sum the dry air mass over grid cells following

$$\begin{array}{}\text{(7)}& M_{{\mathit{\theta }}_{\mathrm{e}}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)=\sum {\left.M_x\left(t\right)\right|}_{{\mathit{\theta }}_{{\mathrm{e}}_x}<{\mathit{\theta }}_{\mathrm{e}}},\end{array}$$

where M_x(t) is the dry air mass of each grid cell x at time t and ${\mathit{\theta }}_{{\mathrm{e}}_x}$ is the equivalent potential temperature of the grid cell. The sum is over all grid cells with ${\mathit{\theta }}_{{\mathrm{e}}_x}$ less than θ_e.

Figure 1 Schematic of the conceptual basis to calculate $M_{{\mathit{\theta }}_{\mathrm{e}}}$. $M_{{\mathit{\theta }}_{\mathrm{e}}}$ of a given θ_e surface is computed by summing all dry air mass with a low equivalent potential temperature in the troposphere of the hemisphere. This calculation yields a unique θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ relation at a given time point.

This calculation yields a unique value of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ for each value of θ_e. We refer to the relationship between θ_e and $M_{{\mathit{\theta }}_{\mathrm{e}}}$ as the “θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ look-up table”, which we generate at daily resolution. We provide this look-up table for each hemisphere computed from ERA-Interim from 1980 to 2018 with a daily resolution and from the lowest to the highest θ_e surface in the troposphere with 1 K intervals (see data availability).

3 Characteristics of $M_{{\mathit{\theta }}_{\mathrm{e}}}$

3.1 Spatial and temporal distribution of $M_{{\mathit{\theta }}_{\mathrm{e}}}$

Figure 2 shows snapshots of the distribution of zonal average θ_e and $M_{{\mathit{\theta }}_{\mathrm{e}}}$ with latitude and pressure at two arbitrary time slices (1 January 2009, 1 July 2009). $M_{{\mathit{\theta }}_{\mathrm{e}}}$ is not continuous across the Equator because it is defined separately in each hemisphere. By definition, each $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surface is exactly aligned with a corresponding θ_e surface, and $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces have the same characteristics as θ_e surfaces, which decrease with latitude and generally increase with altitude. Whereas the zonal average θ_e surfaces vary by up to 20^∘ in latitude over seasons, the meridional displacement of zonal average $M_{{\mathit{\theta }}_{\mathrm{e}}}$ is much smaller, with less than 5^∘ in latitude poleward of 30^∘ N and S, as expected, because the zonal average displacement of atmospheric mass over seasons is small. This small seasonal displacement is closely associated with the seasonality of vertical sloping of θ_e surfaces (Fig. 2). As the mass under each $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surface is always constant, the change in tilt must cause the meridional displacement. In the summer, the tilt is steeper (due to increased deep convection), so $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces move poleward in the lower troposphere but move equatorward in the upper troposphere.

Figure 2 Snapshot of the distribution of (a) zonal average θ_e surfaces on 1 January 2009 (solid lines) and 1 July 2009 (dashed lines) and (b) zonal average $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces on 1 January 2009 (solid lines) and 1 July 2009 (dashed lines). The zonal average tropopause is also shown here for 1 January 2009 (solid black line) and 1 July 2009 (dashed black line). θ_e, $M_{{\mathit{\theta }}_{\mathrm{e}}}$, and the tropopause are computed from ERA-Interim.

Figure 3 $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces as January 2009 average (solid lines) and July 2009 average (dashed lines) for (a) 180^∘ E (mostly over the Pacific Ocean), and (b) 100^∘ E (mostly over the Eurasia land in the Northern Hemisphere). $M_{{\mathit{\theta }}_{\mathrm{e}}}$ and the tropopause are computed from ERA-Interim.

$M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces at given meridians (Fig. 3) in the Northern Hemisphere show clear zonal asymmetry, with larger and more complex displacements compared to the zonal averages, associated with differential heating by land and ocean and orographic stationary Rossby waves (Hoskins and Karoly, 1981; Wills and Schneider, 2018). For example, over the Northern Hemisphere ocean at 180^∘ E (Fig. 3a) and from the summer to winter, $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces move poleward in the middle to high latitudes (e.g., poleward of 45^∘ N) but move equatorward in the mid- to low-latitude lower troposphere (e.g., equatorward of 45^∘ N, 900–700 mbar), with the magnitude smaller than 10^∘ latitude in both. In comparison, over the Northern Hemisphere land at 100^∘ E (Fig. 3b) and from the summer to winter, $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces move equatorward by up to 30^∘ latitude, except in the high-latitude middle troposphere (e.g., poleward of 70^∘ N, ∼ 500 mbar), where the flat $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces lead to slightly poleward displacements. In the Southern Hemisphere, in contrast, the summer-to-winter displacements of the 180 and 100^∘ E sections are similar to the zonal average.

At lower latitudes, the zonal averages of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ and θ_e both exhibit strong secondary maxima near the surface associated with the Hadley circulation (equatorward of 30^∘ N and S) and in the summer, driven by high water vapor. From the contours in Fig. 2, this surface branch of high $M_{{\mathit{\theta }}_{\mathrm{e}}}$ and θ_e appears disconnected from the upper tropospheric branch. In fact, these two branches are connected through air columns undergoing deep convection, which are not resolved in the zonal means shown in Fig. 2 but are resolved in some meridians (e.g., Fig. 3a). We also note that, over the land at 100^∘ E (Fig. 3b), the two disconnected $M_{{\mathit{\theta }}_{\mathrm{e}}}$ and θ_e branches in the Northern Hemisphere summer are displaced poleward compared to the zonal average, consistent with a northward shift of the intertropical convergence zone (ITCZ) over southern Asia. The existence of these two branches may limit some applications of $M_{{\mathit{\theta }}_{\mathrm{e}}}$, as discussed in Sect. 4.

Figure 4 Time series of meridional displacement of selected zonal average θ_e (K) surfaces over a year at (a) 500 mbar, (b) 700 mbar, and (c) 925 mbar. Meridional displacement of selected zonal average $M_{{\mathit{\theta }}_{\mathrm{e}}}$ (10¹⁶ kg) surfaces over a year at (d) 500 mbar, (e) 700 mbar, and (f) 925 mbar. The value of each surface is labeled. θ_e and $M_{{\mathit{\theta }}_{\mathrm{e}}}$ are computed from ERA-Interim. Results shown are for the year 2009.

Figure 4 shows the zonal average meridional displacement of θ_e and $M_{{\mathit{\theta }}_{\mathrm{e}}}$ with a daily resolution. In summer, $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces displace poleward in the lower troposphere but equatorward in the upper troposphere. The displacements in the lower troposphere (925 mbar) are greater in the Northern Hemisphere, where the $M_{{\mathit{\theta }}_{\mathrm{e}}}=\mathrm{140}\times {\mathrm{10}}^{\mathrm{16}}$ kg surface, for example, displaces poleward by 10^∘ in latitude between winter and summer (Fig. 4b). Beside the seasonal variability, Fig. 4 also shows evident synoptic-scale variability.

Since the tilting of θ_e surfaces has an impact on the seasonal displacement of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces, the contribution of different pressure levels to the mass of a given $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin must also vary with season. In Fig. 5, we show these contributions as two daily snapshots on 1 January and 1 July 2009. Low $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins consist of air masses mostly below 500 mbar near the pole. As $M_{{\mathit{\theta }}_{\mathrm{e}}}$ increases, the contribution from the upper troposphere gradually increases while the contribution from the surface to 800 mbar decreases to its minimum at around 100 × 10¹⁶ to 120 × 10¹⁶ kg. The contribution from the surface to 800 mbar increases as $M_{{\mathit{\theta }}_{\mathrm{e}}}$ increases above 120 × 10¹⁶ kg. The mass fraction shows only small variations with season, with the lower troposphere (surface to 800 mbar) contributing slightly less in the low-$M_{{\mathit{\theta }}_{\mathrm{e}}}$ bands and slightly more in the high-$M_{{\mathit{\theta }}_{\mathrm{e}}}$ bands in the summer, which is closely related to the seasonal tilting of corresponding θ_e surfaces.

Figure 5 Snapshots (1 January 2009 and 1 July 2009) of the mass distribution of different $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins from three pressure bins (surface to 800 mbar, 800 to 500 mbar, and 500 mbar to tropopause). $M_{{\mathit{\theta }}_{\mathrm{e}}}$ is computed from ERA-Interim. Low $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins are seen to have larger contributions from the air near the surface, and high $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins have larger contributions from air aloft. Comparing the top and the bottom panels shows that the seasonal differences in pressure contributions are small except for the highest $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins (160 × 10¹⁶–180 × 10¹⁶ kg) and the lowest $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin in the Northern Hemisphere (0–20 × 10¹⁶ kg).

3.2 θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ relationship

Figure 6 compares the temporal variation in $M_{{\mathit{\theta }}_{\mathrm{e}}}$ of several given θ_e surfaces (i.e., θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ look-up table) computed from different reanalysis products for 2009. The deviations are indistinguishable between ERA-Interim and MERRA-2, except near θ_e=340 K, where MERRA-2 is systematically lower than ERA-Interim by 1.5 × 10¹⁶ to 6.5 × 10¹⁶ kg. NCEP2 shows slightly larger deviations from ERA-Interim but by less than 8.5 × 10¹⁶ kg. The products are highly consistent in seasonal variability, and they also show agreement on synoptic timescales. The small difference between products is expected because of different resolutions and methods (Mooney et al., 2011). We expect these differences would be negligible for most applications of $M_{{\mathit{\theta }}_{\mathrm{e}}}$.

Figure 6 Variability in $M_{{\mathit{\theta }}_{\mathrm{e}}}$ of given θ_e surfaces (i.e., θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ look-up table) over a year with a daily resolution in the Northern and Southern Hemisphere. Data from ERA-Interim are shown as a solid line; MERRA-2 data are shown as a dashed line, and NCEP2 data are shown as a dotted line. Results shown are for the year 2009.

Figure 6 shows that, in both hemispheres, $M_{{\mathit{\theta }}_{\mathrm{e}}}$ reaches its minimum in summer and maximum in winter for a given θ_e surface, with the largest seasonality at the lowest θ_e (or $M_{{\mathit{\theta }}_{\mathrm{e}}}$) values. The seasonality decreases as θ_e increases, following the reduction in the seasonality of shortwave absorption at lower latitudes (Li and Leighton, 1993). The seasonality is smaller in the Southern Hemisphere, consistent with the larger ocean area and hence greater heat capacity and transport (Fasullo and Trenberth, 2008; Foltz and McPhaden, 2006). Figure 6 also shows that $M_{{\mathit{\theta }}_{\mathrm{e}}}$ has significant synoptic-scale variability although smaller than the seasonal variability. Synoptic variability is typically larger in winter than summer, as discussed below.

3.3 Relationship to diabatic heating and mass fluxes

A key step of the application of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ for interpreting tracer data is the generation of the look-up table that relates θ_e and $M_{{\mathit{\theta }}_{\mathrm{e}}}$. In this section, we address a tangential question of what controls the temporal variation in the look-up table, which is not necessary for the application but may be of fundamental meteorological interest.

As shown in Appendix A, the temporal variation in the look-up table, ${\dot{M}}_{{\mathit{\theta }}_{\mathrm{e}}}=\frac{\partial }{\partial t}{M}_{{\mathit{\theta }}_{\mathrm{e}}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)$, can be related to underlying mass and heat fluxes according to

$$\begin{array}{}\text{(8)}& {\dot{M}}_{{\mathit{\theta }}_{\mathrm{e}}}=-{\displaystyle \frac{\mathrm{1}}{C_{\mathrm{pd}}}}{\displaystyle \frac{\partial Q_{\mathrm{dia}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}}}+m_{\mathrm{T}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)+m_{\mathrm{E}}\left({\mathit{\theta }}_{\mathrm{e}},t\right),\end{array}$$

where $\frac{\partial Q_{\mathrm{dia}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}}$ (J s⁻¹ K⁻¹) is the effective diabatic heating, integrated over the full θ_e surface per unit width in θ_e; m_T(θ_et) (kg s⁻¹) is the net mass flux across the tropopause; and m_E(θ_et) (kg s⁻¹) is the net mass flux across the Equator, including all air with equivalent potential temperature of less than θ_e. Q_dia has contributions from internal heating without ice formation ($Q_{\mathrm{int}}^{\prime }$), heating from ice formation (Q_ice), sensible heating from the surface (Q_sen), surface evaporation (Q_evap), turbulent diffusion of heat (Q_diff), and turbulent transport of water vapor ($Q_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$) following

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{Q}_{\mathrm{dia}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)& =Q_{\mathrm{int}}^{\prime }\left({\mathit{\theta }}_{\mathrm{e}},t\right)+Q_{\mathrm{ice}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)+Q_{\mathrm{sen}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)\\ & +Q_{\mathrm{evap}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)+Q_{\mathrm{diff}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)+Q_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left({\mathit{\theta }}_{\mathrm{e}},t\right).\end{array}\end{array}$$

The terms Q_evap and $Q_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ are expressed as heating rates by multiplying the underlying water fluxes by L_v(T)∕C_pd. In order to quantify the dominant processes contributing to temporal variation in $M_{{\mathit{\theta }}_{\mathrm{e}}}$, the terms in Eqs. (8) and (9) must be linked to diagnostic variables available in the reanalysis or model products. Although there was no perfect match with any of the three reanalysis products, MERRA-2 provides temperature tendencies for individual processes, which can be converted to heating rates per Eq. 9 following

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial Q_i\left({\mathit{\theta }}_{\mathrm{e}},t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}}}={\displaystyle \frac{C_{\mathrm{pd}}}{\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{e}}}}\sum _x{\left({\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}\right)}_{x,i}{M}_x,\end{array}$$

where i refers to a specific process ($Q_{\mathrm{int}}^{\prime }$, Q_ice, etc.), ${\left(\frac{\mathrm{d}T}{\mathrm{d}t}\right)}_x$ (K s⁻¹) is the temperature tendency of grid cell x, M_x (kg) is the mass of grid cell x, and Δθ_e is the width of the θ_e surface.

Table 1 Correspondence of heating variables between our derivation (Eq. 9) and MERRA-2.

There are five heating terms provided in the MERRA-2 product, which we can approximately relate to terms in Eq. (9), as shown in Table 1. The first three terms (Q_rad, Q_dyn, and Q_ana) can be summed to yield $Q_{\mathrm{int}}^{\prime }$; the fourth (Q_trb) is equal to the sum of Q_diff and Q_sen; and the fifth (Q_mst) approximates the sum of Q_ice and Q_evap. MERRA-2 does not provide terms corresponding to $Q_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ or Q_evap, but Q_mst represents heating due to moist processes, which includes Q_ice plus water vapor evaporation and condensation within the atmosphere. This water vapor evaporation and condensation should be approximately equal to Q_evap with a small time lag when integrated over a θ_e surface because mixing is preferentially along θ_e surfaces and water vapor released into a θ_e surface by surface evaporation will tend to transport and precipitate from the same θ_e surface within a short time period (Bailey et al., 2019). Thus, the MERRA-2 term for heating by moist processes (Q_mst) should approximate Q_ice+Q_evap.

Figure 7 (a) Temporal variation in $M_{{\mathit{\theta }}_{\mathrm{e}}}$ in the Northern Hemisphere at θ_e=300 K computed by integrating air mass (blue line) and estimated from the sum of five heating terms (Table 1) in MERRA-2 (black line). (b) The heating variables are decomposed into five contributions as indicated (see Table 1). Results shown are for the year 2009.

Figure 7a compares the temporal variation in ${\dot{M}}_{{\mathit{\theta }}_{\mathrm{e}}}$ computed by integrating the dry air mass (i.e., θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ look-up table) with $M_{{\mathit{\theta }}_{\mathrm{e}}}$ computed from the sum of the diabatic heating terms from MERRA-2 (via Eqs. 8 to 10). The comparison focuses on the θ_e=300 K surface, which does not intersect with the Equator or tropopause, so the two mass flux terms (m_T, m_E) vanish. These two methods have a high correlation at 0.71. We do not expect perfect agreement because ${\dot{M}}_{{\mathit{\theta }}_{\mathrm{e}}}$ computed by the sum of heating neglects turbulent water vapor transport ($Q_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$), and only approximates Q_evap as discussed above. This relatively good agreement nevertheless demonstrates that the formulation based on MERRA-2 heating terms includes the dominant processes that drive temporal variations in the look-up table. Figure 7a shows poorer agreement from late August to October, which we also find in other years (Figs. S1 and S2 in the Supplement) and on lower (e.g., θ_e=290 K, Fig. S3) but not higher (e.g., θ_e=310 K, Fig. S4) surfaces, where the two methods agree better. The poor agreement may reflect a partial breakdown of the assumption that Q_mst approximates the sum of Q_ice and Q_evap, but further analysis is beyond the scope of this study.

Figure 7b further breaks down the sum of the heating terms in Eqs. (8) and (10) from MERRA-2 into individual components. Each term clearly displays variability on synoptic to seasonal scales. To quantify the contribution of different terms on the different timescales, we separate each term into a seasonal and synoptic component, where the seasonal component is derived by a two-harmonic fit with a constant offset and the synoptic component is the residual. We estimate the fractional contribution of each heating term on seasonal and synoptic timescales separately in Table 2, using the method in Sect. S1 in the Supplement. On the seasonal timescale, the variance is dominated by radiative heating and cooling of the atmosphere and moist processes (including both ice formation and extra water vapor from surface evaporation) together, with prominent counteraction between them. On the synoptic timescale, dissipation of the kinetic energy of turbulence dominates the variance.

Table 2 Fractional contribution of the individual heating terms in Fig. 7b to their sum for θ_e=300 K. The analysis is done separately on synoptic and seasonal components. The seasonal component is based on a two-harmonic fit, and the synoptic component is defined as the residual. The fractional contributions sum to 1, while a positive contribution means in phase and negative contribution means anti-phase. A contribution of an absolute value that is bigger than 1 illustrates that the variability in the heating term is larger than the variability in the sum on the corresponding timescale.

Similar analyses on different θ_e surfaces and in different years (Figs. S1 to S4) all show that a combination of radiative heating and moist processes dominates the temporal variation in $M_{{\mathit{\theta }}_{\mathrm{e}}}$ on the seasonal timescale, while dissipation of the kinetic energy of turbulence dominates on the synoptic timescale.

4 Applications of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ as an atmospheric coordinate

To illustrate the potential application of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ for interpreting sparse data, we focus on the seasonal cycle of CO₂ in the Northern Hemisphere as resolved by two series of global airborne campaigns, HIPPO and ATom. HIPPO consisted of five campaigns between 2009 and 2011, and ATom consisted of four campaigns between 2016 and 2018. Each campaign covered from ∼ 150 to ∼ 14 000 m and from nearly pole to pole, along both northbound and southbound transects. On HIPPO, both transects were over the Pacific Ocean, while on ATom, southbound transects were over the Pacific Ocean and northbound transects were over the Atlantic Ocean. The flight tracks are shown in Fig. 8a. We aggregate data from each campaign into northbound and southbound transects within each hemisphere but only use data from the Northern Hemisphere. We only consider tropospheric observations by excluding measurements from the stratosphere, which is defined by observed water vapor of less than 50 ppm and either O₃ greater than 150 ppb or detrended N₂O to the reference year of 2009 of less than 319 ppb. Water vapor and O₃ were measured by the NOAA UCATS (Unmanned Aerial Systems Chromatograph for Atmospheric Trace Species; Hurst, 2011) instrument and were interpolated to a 10 s resolution. N₂O was measured by the Harvard QCLS (Quantum Cascade Laser System; Santoni et al., 2014) instrument. Furthermore, we exclude all near-surface observations within ∼ 100 s of takeoffs and within ∼ 600 s of landings as well as missed approaches, which usually show high CO₂ variability due to strong local influences. In situ measurements of CO₂ were made by three different instruments on both HIPPO and Atom. Of these, we use the CO₂ measurements made by the NCAR Airborne Oxygen Instrument (AO2) with a 2.5 s measurement interval (Stephens et al., 2020), for consistency with planned future applications of APO (atmospheric potential oxygen) computed from AO2. The differences between instruments are small for our application (Santoni et al., 2014). The data used in this study are averaged to a 10 s resolution, and we show the detrended CO₂ values along each airborne campaign transect for the Northern Hemisphere in Fig. 8b. Since we focus on the seasonal cycle of CO₂, all airborne observations are detrended by subtracting an interannual trend fitted to CO₂ measured at the Mauna Loa Observatory (MLO) by the Scripps CO₂ Program. This trend is computed by a stiff cubic spline function plus four-harmonic terms with linear gain to the MLO record. $M_{{\mathit{\theta }}_{\mathrm{e}}}$ is computed from ERA-Interim in this section.

Figure 8 (a) HIPPO and ATom horizontal flight tracks colored by campaigns. (b) Latitude and pressure cross section of detrended CO₂ of each airborne campaign transect. CO₂ is detrended by subtracting the MLO stiff cubic spline trend, which is computed by a stiff cubic spline function plus four-harmonic functions with linear gain to the MLO record.

4.1 Mapping Northern Hemisphere CO₂

A conventional method to display seasonal variations in CO₂ from airborne data is to plot time series of the data at a given location or latitude and different pressure levels (Graven et al., 2013; Sweeney et al., 2015). In Fig. 9, we compare this method using HIPPO and ATom airborne data, binning and averaging the data from each airborne campaign transect by pressure and latitude bins, with our new method, binning the data by pressure and $M_{{\mathit{\theta }}_{\mathrm{e}}}$. For each latitude bin, we choose a corresponding $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin which has approximately the same meridional coverage in the lower troposphere. We remind the reader that $M_{{\mathit{\theta }}_{\mathrm{e}}}$ decreases poleward while also generally increasing with altitude (Figs. 2 to 4).

Figure 9 Seasonal cycles of airborne Northern Hemisphere CO₂ data sorted by (a) $M_{{\mathit{\theta }}_{\mathrm{e}}}$–pressure bins and (b) latitude–pressure bins. $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins (10¹⁶ kg) and latitude bins are shown at the top of each panel. Pressure bins are colored. The latitude bounds are chosen to approximate the meridional coverage of each corresponding $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin in the lower troposphere. The seasonal cycle at MLO from 2009 to 2018 is shown in the 90–110 $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin panel, which spans the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ of the station. Airborne observations are first grouped into $M_{{\mathit{\theta }}_{\mathrm{e}}}$–pressure or latitude–pressure bins and then averaged for each airborne campaign transect, shown as points. We filter out the points averaged from fewer than twenty 10 s observations. The seasonal cycle of airborne data and MLO (2009–2018) are computed by a two-harmonic fit to the detrended time series. The 1σ variability about the seasonal cycle fits for each $M_{{\mathit{\theta }}_{\mathrm{e}}}$–pressure or latitude–pressure bin is labeled at the top of each panel. These 1σ values are based on the distribution of all binned observations (not shown), rather than on the distribution of average CO₂ of each bin and airborne campaign transect (shown).

As shown in Fig. 9, the transect averages of detrended CO₂ (shown as points) from both binning methods resolve well-defined seasonal cycles (based on two-harmonic fit) in all bins, with higher amplitudes near the surface (low pressure) and at high latitudes (low $M_{{\mathit{\theta }}_{\mathrm{e}}}$). However, binning by $M_{{\mathit{\theta }}_{\mathrm{e}}}$ leads to much smaller variations in the mean seasonal cycle (shown as solid curves) with pressure, as expected, because moist isentropes are preferential surfaces for mixing. Also, within individual pressure bins, the short-term variability relative to the mean cycles based on the distribution of all detrended observations (not shown as points but denoted as 1σ values in Fig. 9) is smaller when binning by $M_{{\mathit{\theta }}_{\mathrm{e}}}$ (F test, p < 0.01), except in the lower troposphere of the highest $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin (90 × 10¹⁶–110 × 10¹⁶ kg). The smaller short-term variability is expected because $M_{{\mathit{\theta }}_{\mathrm{e}}}$ tracks the synoptic variability in the atmosphere. When binning by latitude, the smallest short-term variability is found at the lowest bin (surface–800 mbar) and the largest short-term variability is found in the highest bin (500 mbar tropopause), except the highest latitude bin (45–55^∘ N). When binning by $M_{{\mathit{\theta }}_{\mathrm{e}}}$, in contrast, the short-term variability in the middle pressure bin is always smaller than the higher and lower pressure bins (F test, p < 0.01), except for the 50-to-70 $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin, where the difference between the lowest and middle pressure bins is not significant (based on 1σ levels). The lower variability in the middle troposphere may reflect the suppression of variability from synoptic disturbances, leaving a clearer signal of the influence of surface fluxes of CO₂ and stratosphere–troposphere exchanges. We compare the variance of detrended airborne observations within each $M_{{\mathit{\theta }}_{\mathrm{e}}}$–pressure bin with its fitted value. The fitted seasonal cycle of each bin explains 63.2 % to 90.5 % of the variability for different bins, with higher fractions in the middle troposphere.

Figure 9 also shows the CO₂ seasonal cycle at MLO, which falls within a single $M_{{\mathit{\theta }}_{\mathrm{e}}}$–pressure bin (90 × 10¹⁶–110 × 10¹⁶ kg, 500–800 mbar) at all seasons. Although the airborne data in this bin span a wide range of latitudes (∼ 10–75^∘ N), the seasonal cycle averaged over this bin is very similar to the cycle at MLO (airborne cycle leads by ∼ 10 d with 1.0 % lower amplitude). This small difference is within the 1σ uncertainty in our estimation from airborne observation, and some difference is expected, since we choose a $M_{{\mathit{\theta }}_{\mathrm{e}}}$–pressure bin wider than the seasonal variation in $M_{{\mathit{\theta }}_{\mathrm{e}}}$ and pressure at MLO.

Figure 10 CO₂ seasonal cycles of multiple surface stations (2009–2018) compared to seasonal cycles of airborne observations averaged over corresponding $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin. The choice of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin is to approximate the range of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ at each corresponding surface station and is shown at the top of each panel. Daily $M_{{\mathit{\theta }}_{\mathrm{e}}}$ of the station is computed from ERA-Interim, based on its location. We detrend station and airborne observations by subtracting the MLO stiff cubic spline trend. We compute average detrended CO₂ for each airborne campaign transect and each $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bin, shown as black points. The seasonal cycles are computed from a two-harmonic fit, with the seasonal amplitude (Amp) shown in the upper right of each panel.

It is also of interest to examine how CO₂ data from surface stations fit into the framework based on $M_{{\mathit{\theta }}_{\mathrm{e}}}$. Figure 10 compares the CO₂ seasonal cycle of five NOAA surface stations (Cooperative Global Atmospheric Data Integration Project, 2019) with the cycle from the airborne observations binned into selected $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins. These surface stations are chosen to be representative of different $M_{{\mathit{\theta }}_{\mathrm{e}}}$ ranges. For the comparison, we chose $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins that span the seasonal maximum and minimum $M_{{\mathit{\theta }}_{\mathrm{e}}}$ value of the station. These bins are narrower than the bins used in Fig. 9, in order to sharply focus on the latitude of the station. To maximize sampling coverage, we bin the airborne data only by $M_{{\mathit{\theta }}_{\mathrm{e}}}$ without pressure sub-bins. For mid- and high-latitude surface stations (right three panels), the seasonal amplitude of station CO₂ and corresponding airborne CO₂ are close (within 4 %–5 %), while airborne cycles lag by 2–3 weeks. The lag presumably represents the slow mixing from the mid-latitude surface to the high-latitude mid-troposphere (Jacob, 1999). In contrast, for low-latitude stations (left two panels) which generally sample trade winds, the seasonal cycles differ significantly, indicating that the air sampled at these stations is not rapidly mixed along surfaces of constant $M_{{\mathit{\theta }}_{\mathrm{e}}}$ or θ_e with air aloft. As mentioned above (Sect. 3.1), surfaces of high $M_{{\mathit{\theta }}_{\mathrm{e}}}$ within the Hadley circulation have two branches, one near the surface and one aloft. A timescale of several months for transport from the lower to the upper branch can be estimated from the known overturning flows based on air mass flux stream functions (Dima and Wallace, 2003). This delay, plus strong mixing and diabatic effects (Miyazaki et al., 2008), ensures that the lower and upper branches are not well connected on seasonal timescales. Our results nevertheless demonstrate that the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ framework combining airborne and surface data could help understanding of details of atmospheric transport both along and across θ_e surfaces.

4.2 Computing the hemispheric mass-weighted average CO₂ mole fraction

We next illustrate the use of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ for computing the mass-weighted average of a long-lived chemical tracer by performing this exercise for CO₂ in the Northern Hemisphere. We calculate the Northern Hemisphere tropospheric mass-weighted average CO₂ from each airborne transect using a method that assumes that CO₂ is uniformly mixed on θ_e surfaces throughout the hemisphere (Barnes et al., 2016; Parazoo et al., 2011, 2012). We exclude airborne observation from HIPPO-1 Northbound due to the lack of data north of 40^∘ N. We use the θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ look-up table of the corresponding date to assign a value of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ to each observation based on its θ_e. The observations for each transect are then sorted by $M_{{\mathit{\theta }}_{\mathrm{e}}}$. The hemispheric average CO₂ is calculated by trapezoidal integration of CO₂ as a function of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ and divided by the total dry air mass as computed from the corresponding range of $M_{{\mathit{\theta }}_{\mathrm{e}}}$.

To illustrate the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method, we choose HIPPO-1 Southbound and show CO₂ measurements and ΔCO₂ atmospheric inventory (Pg) as a function of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ in Fig. 11. The Northern Hemisphere tropospheric average detrended ΔCO₂ is computed by integrating the area under the curve (subtracting negative contributions) and dividing by the maximum value of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ within the hemisphere (here 195.13 × 10¹⁶ kg). This yields a mass-weighted average detrended ΔCO₂ of 1.13 ppm for the full troposphere of the Northern Hemisphere. The trapezoidal integration has a high accuracy because the data are dense over $M_{{\mathit{\theta }}_{\mathrm{e}}}$. The ΔCO₂ atmospheric inventory is dominated by the domain $M_{{\mathit{\theta }}_{\mathrm{e}}}$ < 120 ×10¹⁶ kg (mid-latitude to high latitude), which has a large CO₂ seasonal cycle driven by a temperate and boreal ecosystem, with less than 4.1 % contributed by the additional ∼ 38.8 % of the air mass outside this domain in the low latitudes or upper troposphere (Fig. 11b), where ΔCO₂ differs less from the subtracted baseline.

Figure 11 (a) Detrended CO₂ measurements from HIPPO-1 Southbound (from 12 to 17 January 2009) plotted as a function of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ in the Northern Hemisphere. The data are detrended by subtracting the MLO stiff cubic spline trend. Individual points are connected by straight line segments, and the area under the resulting curve is shaded. We note that the area under the curve has units of parts per million × kilograms, and dividing this by the total dry air mass (i.e., the range of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ of the integral) gives the parts per million unit because the mass of dry air is proportional to the moles of dry air. The Northern Hemisphere average of 1.13 ppm is indicated by the dashed line. (b) Integral of the data in panel (a), rescaled from parts per million to petagrams, integrating from $M_{{\mathit{\theta }}_{\mathrm{e}}}=\mathrm{0}$ to a given $M_{{\mathit{\theta }}_{\mathrm{e}}}$ value.

We compute Northern Hemisphere mass-weighted average detrended ΔCO₂ for each airborne campaign transect and fit the time series to a two-harmonic fit to estimate the seasonal cycle (Fig. 12). We find that the cycle has a seasonal amplitude of 7.9 ppm and a downward zero-crossing on Julian day 179, where the latter is defined as the date when the detrended seasonal cycle changes from positive to negative.

Figure 12 Comparison between the CO₂ seasonal cycle of Northern Hemisphere tropospheric average computed from airborne observation and the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method (black points and line) and the mean cycle at MLO measured by Scripps CO₂ Program from 2009 to 2018 (red line). Both are detrended by subtracting a stiff cubic spline trend at MLO. We then compute mass-weighted average detrended CO₂ for each airborne campaign transect, shown as black points, with campaigns and transects presented as different shapes. The seasonal cycle of both are computed by a two-harmonic fit to the detrended time series. The 1σ variability in the detrended average CO₂ values about the fit line is shown on the lower right. The first half year is repeated for clarity.

To address the error in our estimation of the Northern Hemisphere mass-weighted average CO₂ seasonal cycle from HIPPO and Atom airborne observation, we consider two main sources: (1) irreproducibility in the CO₂ measurements and (2) limited coverage in space and time. For the first contribution, we compute the difference between mass-weighted average CO₂ from AO2 and mean mass-weighted average CO₂ from Harvard QCLS, Harvard OMS, and NOAA Picarro for each airborne campaign transect, while masking values that are missing in any of these datasets. We compute the standard deviation of these differences (±0.15 ppm) for the mass-weighted average CO₂ of each airborne campaign transect as the 1σ level of uncertainty. We further compute the uncertainties for the seasonal amplitude of ±0.11 ppm and for the downward zero-crossing of ±0.83 d, which are calculated from 1000 iterations of the two-harmonic fit, allowing for random Gaussian uncertainty (σ $=\pm$0.15 ppm) for each transect.

Table 3 RMSE, seasonal amplitude, and day of year of the downward zero-crossing of each simulation based on the Jena CO₂ inversion. The true value (daily average CO₂) is computed by integrating over all tropospheric grid cells of the Jena CO₂ inversion, while troposphere is defined by PVU < 2 from ERA-Interim. Seasonal amplitude and downward zero-crossing of true average and each simulation is computed from two-harmonic fit to the detrended value, which is detrended by subtracting the MLO cubic stiff spline. Subsampling with randomly retaining a certain fraction of data are conducted by randomly subsampling 1000 times, thus, the seasonal amplitude and day of year of the downward zero-crossing is computed as the mean ± standard deviation of the 1000 iterations.

^∗ Each simulation yields 17 data points of different dates over the seasonal cycle from 17 airborne campaign transects. RMSE of each simulation is computed with respect to the true value.

For the contribution to the error in the amplitude and phase from limited special and temporal coverage, we use simulated CO₂ data from the Jena CO₂ inversion Run ID s04oc v4.3 (Rödenbeck et al., 2003). This model includes full atmospheric fields from 2009 to 2018, which we detrend using the cubic spline fit to the observed MLO trend. From these detrended fields, we compute the climatological cycle of the Northern Hemisphere average by integrating over all tropospheric grid cells (cutoff at PVU = 2) to produce a daily time series of the hemispheric mean, which we take as the model “truth”. We fit a two-harmonic function to this true time series to compute a true climatological cycle over the 2009–2018 period (Table 3), which is our target for validation. We then subsample the Jena CO₂ inversion along the HIPPO and ATom flight tracks and process the data similarly to the observations, using the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method and a two-harmonic fit. The comparison shows that the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method yields an amplitude which is 1 % too large and yields a downward zero-crossing date which is 6 d too late. We view these offsets as systematic biases, which we correct from the observed amplitude and phase reported above. The uncertainties in these biases are hard to quantify, but we take ±100 % as a conservative estimate. We thus allow an additional random error of ±0.08 ppm in amplitude and ±6.0 d in downward zero-crossing for uncertainty in the bias. Combining the random and systematic error contributions leads to a corrected Northern Hemisphere tropospheric average CO₂ seasonal cycle amplitude of 7.8 ± 0.14 ppm and downward zero-crossing of 173 ± 6.1 d. This corrected cycle is an estimate of the climatological average from 2009 to 2018.

The error due to limited spatial and temporal coverage can be divided into three components: limited seasonal coverage (17 transects over the climatological year), limited interannual coverage (sampling particular years instead of all years), and limited spatial coverage (under-sampling the full hemisphere). We quantify the combined biases due to both limited seasonal and limited interannual coverage by comparing the two-harmonic fit of the full true daily time series of the hemispheric mean to a two-harmonic fit of those data subsampled on the actual mean sampling dates of the 17 flight tracks. We isolate the bias associated with limited seasonal coverage by repeating this calculation, replacing the true daily time series with the daily climatological cycle. The bias associated with limited spatial coverage is quantified as the residual. Combining these results, we estimate that the limited seasonal, interannual, and spatial coverage account for biases in the downward zero-crossing of 1.1, 1.4, and 3.5 d respectively, all in the same direction (too late). The seasonal amplitude biases due to individual components are all small (< 0.5 %).

It is of interest to compare our estimate of the Northern Hemisphere average cycle with the cycle at Mauna Loa, which is also broadly representative of the hemisphere. Our comparison in Fig. 12 shows small but significant differences in both amplitude and phase, with the MLO amplitude being ∼ 11.5 % smaller than the hemispheric average and lagging in phase by ∼ 1 month. There are also differences in the shape of the cycle, with the MLO cycle rising more slowly from October to February but more quickly from February to May. These features at least partly reflect variations in the transport of air masses to the station (Harris et al., 1992; Harris and Kahl, 1990).

In Fig. 13, we compare the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method with an alternate latitude–pressure weighted average method, with no correction for synoptic variability. For this method, we bin flight track subsampled Jena CO₂ inversion data into sin(latitude)-pressure bins with 0.01 and 25 mbar as intervals respectively, while all bins without data are filtered. We further compute weighted average CO₂ for each airborne campaign transect. The root-mean-square errors (RMSEs) to the true average of the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method are ±0.32 and ±0.27 ppm for the HIPPO and ATom campaigns respectively, which are smaller than the RMSEs of the simple latitude–pressure weighted average method at ±0.82 and ±0.53 ppm.

Figure 13 Comparison between the Northern Hemisphere average CO₂ from full integration of the simulated atmospheric fields from the Jena CO₂ inversion (cutoff at PVU = 2) and from two methods that use the same simulated data subsampled with HIPPO or ATom coverage: (1) the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method (blue) and (2) simple integration by sin(latitude)–pressure (red). We divide the comparison into HIPPO (a, c) and ATom (b, d) temporal coverage. Panels (c) and (d) shows the error for individual tracks using alternate subsampling methods.

We also evaluate the biases in the hemispheric average seasonal cycles computed with the simple latitude–pressure weighted average method. As summarized in Table 3, the latitude–pressure weighted average method yields a larger error in seasonal amplitude ($M_{{\mathit{\theta }}_{\mathrm{e}}}$ method 1.0 % too large, latitude–pressure method 20.8 % too large), while both methods show a similar phasing error (6 to 7 d late). The larger error associated with the latitude–pressure weighted average method is consistent with strong influence of synoptic variability. This synoptic variability could potentially be corrected using model simulations of the 3-dimensional CO₂ fields (Bent, 2014). The $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method appears advantageous because it accounts for synoptic variability and easily yields a hemispheric average by directly integrating over $M_{{\mathit{\theta }}_{\mathrm{e}}}$.

The relative success of the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method in yielding accurate hemispheric averages using HIPPO and ATom data is attributable partly to the extensive data coverage. To explore the coverage requirement for reliably resolving hemispheric averages, we also test the integration method when applied to simulated data with lower coverage. We start with the same coverage as for ATom and HIPPO but select only subsets of the points in four groups: poleward of 30^∘ N, Equator to 30^∘ N, surface to 600 mbar, and 600 mbar to tropopause. We also examine whether we can only utilize observation along the Pacific transect by excluding measurements along the Atlantic transects (ATom northbound). We further explore the impact of reduced sampling density by subsampling the Jena CO₂ inversion based on the spatial coverage of the Medusa sampler, which is an airborne flask sampler that collected 32 cryogenically dried air samples per flight during HIPPO and ATom (Stephens et al., 2020). We further randomly retain 10 %, 5 %, and 1 % of the full flight track subsampled data, repeating each ratio with 1000 iterations. We compute the detrended average CO₂ from these nine simulations by the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method and then compute the RMSE relative to the detrended true hemispheric average, together with the seasonal magnitude and day of year of the downward zero-crossing, as summarized in Table 3. HIPPO-1 Northbound is excluded in all these simulations. The number of data points of each simulation and number of observations of the original HIPPO and ATom datasets are summarized in Table S1. These results show that limiting sampling to either equatorward or poleward of 30^∘ N yields significant error (24.3 % smaller and 24.9 % larger seasonal amplitude respectively). Additionally, there is a ∼ 25 d lag in phase if sampling is limited to equatorward of 30^∘ N. However, restricting sampling to be exclusively above or below 600 mbar or only along the Pacific transect does not lead to significant errors. Randomly reducing the sampling by 10- to 100-fold or only keeping Medusa spatial coverage also has minimal impact. This suggests that, to compute the average CO₂ of a given region, it may be sufficient to have low sampling density provided that the measurements adequately cover the full range in θ_e (or $M_{{\mathit{\theta }}_{\mathrm{e}}}$).

5 Discussion and summary

We have presented a transformed isentropic coordinate, $M_{{\mathit{\theta }}_{\mathrm{e}}}$, which is the total dry air mass under a given θ_e surface in the troposphere of the hemisphere. $M_{{\mathit{\theta }}_{\mathrm{e}}}$ can be computed from meteorological fields by integrating dry air mass under a specific θ_e surface, and different reanalysis products show a high consistency. The θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ relationship varies seasonally due to seasonal heating and cooling of the atmosphere via radiative heating and moist processes. The seasonality in the relationship is greater at low θ_e compared to high θ_e and is greater in the Northern Hemisphere than in the Southern Hemisphere. The θ_e–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ relationship also shows synoptic-scale variability, which is mainly driven by the dissipation of the kinetic energy of turbulence. $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces show much less seasonal displacement with latitude and altitude than surfaces of constant θ_e while being parallel and exhibiting essentially identical synoptic-scale variability. As a coordinate for mapping tracer distributions, $M_{{\mathit{\theta }}_{\mathrm{e}}}$ shares with θ_e the advantages of following displacements due to synoptic disturbances and aligning with surfaces of rapid mixing. $M_{{\mathit{\theta }}_{\mathrm{e}}}$ has the additional advantages of being approximately fixed in space seasonally, which allows mapping to be done on seasonal timescales, and having units of mass, which provides a close connection with atmospheric inventories.

As a coordinate, $M_{{\mathit{\theta }}_{\mathrm{e}}}$ is probably better viewed as an alternative to latitude, due to its nearly fixed relationship with latitude over season, rather than as an alternative to altitude (or pressure), as typically done for potential temperature (Miyazaki et al., 2008; Miyazaki and Iwasaki, 2005; Parazoo et al., 2011; Tung, 1982; Yang et al., 2016). Even though the contours of constant $M_{{\mathit{\theta }}_{\mathrm{e}}}$ extend over a wide range of latitudes (from low latitudes at the Earth's surface to high latitudes aloft), a close association with latitude is provided by the point of contact with the Earth's surface. Also, $M_{{\mathit{\theta }}_{\mathrm{e}}}$ is nearly always monotonic with latitude (increasing equatorward), while it is not necessarily monotonic with altitude in the lower troposphere (Figs. 2 and 3).

As a first application, we have illustrated using $M_{{\mathit{\theta }}_{\mathrm{e}}}$ the seasonal variation in CO₂ in the Northern Hemisphere, with data from the HIPPO and ATom airborne campaigns. This application shows that $M_{{\mathit{\theta }}_{\mathrm{e}}}$ has several advantages as a coordinate compared to using latitude: (1) variations in CO₂ with pressure are smaller at fixed $M_{{\mathit{\theta }}_{\mathrm{e}}}$ than at fixed latitude, and (2) the scatter about the mean CO₂ seasonal cycle is smaller when sorting data into pressure–$M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins than into pressure–latitude bins. We have also shown that, at middle and high latitudes, the CO₂ seasonal cycles that are resolved in the airborne data (binned by $M_{{\mathit{\theta }}_{\mathrm{e}}}$ but not pressure) are very similar to the cycles observed at surface stations at the appropriate latitude, with a phase lag of ∼ 2 to 3 weeks. At lower latitudes, CO₂ cycles in the airborne data (binned similarly by $M_{{\mathit{\theta }}_{\mathrm{e}}}$) are less consistent with surface data, as expected due to slow transport and diabatic processes within the Hadley circulation. For characterizing the patterns of variability in airborne CO₂ data, we expect the advantages of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ over latitude will be greatest for sparse datasets, allowing data to be binned more coarsely with pressure or elevation while still resolving features of large-scale variability, such as seasonal cycles or gradients with latitude.

As a second application, we use $M_{{\mathit{\theta }}_{\mathrm{e}}}$ to compute the Northern Hemisphere tropospheric average CO₂ from the HIPPO and ATom airborne campaigns by integrating CO₂ over $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces. With a small correction for systematic biases induced by limited hemispheric coverage of the HIPPO and ATom flight tracks, we report a seasonal amplitude of 7.8 ± 0.14 ppm and a downward zero-crossing on Julian day 173 ± 6.1. This hemispheric average cycle may prove valuable as a target for validation of models of surface CO₂ exchange.

Our analysis also clarifies that computing hemispheric averages with the $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method depends on adequate spatial coverage. The coverage provided by the HIPPO and ATom campaigns appears more than adequate for computing the average seasonal cycle of CO₂ in the Northern Hemisphere, and the errors for this application remain small if the coverage is limited to either above or below 600 mbar or reduced to retain only 1 % of the measurements. Most critical is maintaining coverage in latitude or $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces. The $M_{{\mathit{\theta }}_{\mathrm{e}}}$ integration method of computing hemispheric averages assumes that the tracer is uniformly distributed and instantly mixed on θ_e ($M_{{\mathit{\theta }}_{\mathrm{e}}}$) surfaces. We have shown that systematic gradients in CO₂ are resolved with pressure at fixed $M_{{\mathit{\theta }}_{\mathrm{e}}}$, which reflects the finite rates of dispersion on θ_e surfaces. Further improvements to the integration method seem possible by integrating separately over different pressure levels, taking account of the different mass fraction in different pressure bins (e.g., Fig. 5). The need is especially relevant for high $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins which are less completely mixed and which tend to intersect the Equator or have separate surface branches. For these $M_{{\mathit{\theta }}_{\mathrm{e}}}$ bins, it would be more appropriate to integrate over M_θ in the upper and lower atmosphere separately. This complication is of minor importance for computing the mass-weighted average CO₂ cycle, because the cycle of CO₂ is small in these air masses.

The definition of $M_{{\mathit{\theta }}_{\mathrm{e}}}$ requires horizontal and vertical boundaries for the integration of dry air masses. We use the dynamic tropopause (based on potential vorticity units) and the Equator as boundaries, which is appropriate for integrating tropospheric inventories in a hemisphere. Other boundaries may be more appropriate for other applications. For example, $M_{{\mathit{\theta }}_{\mathrm{e}}}$ could be computed from the lowest θ_e surface in the Southern Hemisphere with a latitude cutoff at 30^∘ S, to apply to airborne observations only over the Southern Ocean. On the other hand, the boundary choice only influences $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces that actually intercept the boundaries, making the choice less important at high latitudes in the lower troposphere (lowest $M_{{\mathit{\theta }}_{\mathrm{e}}}$ surfaces). Some tropospheric applications may also benefit by integrating over dry potential temperature (θ) rather than θ_e.

Based on our promising results for CO₂, we expect that $M_{{\mathit{\theta }}_{\mathrm{e}}}$ may be usefully applied as a coordinate for mapping and computing atmospheric inventories of many tracers, such as O₂ ∕ N₂, N₂O, CH₄, and the isotopes of CO₂, whose residence time is long compared to the timescale for mixing along isentropes. $M_{{\mathit{\theta }}_{\mathrm{e}}}$ may also prove useful in the design phase of airborne campaigns to ensure strategic coverage. Our results show that, to study the seasonal cycle of a tracer on a hemispheric scale, it is critical to have well-distributed sampling in $M_{{\mathit{\theta }}_{\mathrm{e}}}$.

Appendix A: Temporal variation in $M_{{\mathit{\theta }}_{\mathrm{e}}}$

Following Walin's derivation for cross-isothermal volume flow in the ocean (Walin, 1982), we show how ${\dot{M}}_{{\mathit{\theta }}_{\mathrm{e}}}=\frac{\partial }{\partial t}{M}_{{\mathit{\theta }}_{\mathrm{e}}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)$ can be related to energy and mass fluxes. We start by deriving the relationship for M_θ (based on potential temperature θ) but later generalize to apply to $M_{{\mathit{\theta }}_{\mathrm{e}}}$.

Table A1 Definition of variables.

All definitions are summarized in Table A1, and Fig. A1 is the schematic diagram of mass and energy flux.

Figure A1 Illustration of terms defined in Table A1. Shaded area denotes the region R(θ,t) with θ^′ lower than θ, which is the area of mass integration to yield M_θ. The curve denotes a given θ or M_θ surface.

All mass and heat fluxes into region R(θ,t) are counted as positive. The heat fluxes through the tropopause, Equator, and surface of region R(θ,t) can be divided into an advective (F(θ,t)) and a turbulent (D(θ,t)) component. Integrating over the tropopause and equatorial boundary, we have

$$\begin{array}{}\text{(A1)}& {\displaystyle }{Q}_{\mathrm{T}}(\mathit{\theta },t)=C_{\mathrm{pd}}\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial F_{\mathrm{T}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}{\mathit{\theta }}^{\prime }\mathrm{d}{\mathit{\theta }}^{\prime }+\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial D_{\mathrm{T}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime },\text{(A2)}& {\displaystyle }{Q}_{\mathrm{E}}(\mathit{\theta },t)=C_{\mathrm{pd}}\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial F_{\mathrm{E}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}{\mathit{\theta }}^{\prime }\mathrm{d}{\mathit{\theta }}^{\prime }+\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial D_{\mathrm{E}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime },\text{(A3)}& {\displaystyle }{Q}_{\mathrm{I}}(\mathit{\theta },t)=C_{\mathrm{pd}}\cdot F_{\mathrm{I}}(\mathit{\theta },t)\cdot \mathit{\theta }+D_{\mathrm{I}}(\mathit{\theta },t),\end{array}$$

where C_pd is the heat capacity of dry air in units of J kg⁻¹ K⁻¹.

Based on the continuity of mass and energy for region R(θ,t), we obtain

$$\begin{array}{}\text{(A4)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}{M}_{\mathit{\theta }}(\mathit{\theta },t)& =F_{\mathrm{T}}(\mathit{\theta },t)+F_{\mathrm{E}}(\mathit{\theta },t)+F_{\mathrm{I}}(\mathit{\theta },t)\\ & =\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial F_{\mathrm{T}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime }+\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial F_{\mathrm{E}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime }\\ & +F_{\mathrm{I}}(\mathit{\theta },t),\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}{C}_{\mathrm{pd}}{\displaystyle \frac{\partial }{\partial t}}& \underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial M_{\mathit{\theta }}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}{\mathit{\theta }}^{\prime }\mathrm{d}{\mathit{\theta }}^{\prime }=Q_{\mathrm{T}}(\mathit{\theta },t)+Q_{\mathrm{E}}(\mathit{\theta },t)\\ & +Q_{\mathrm{I}}(\mathit{\theta },t)+\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial Q_{\mathrm{s}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime }\\ & +\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial Q_{\mathrm{int}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime }.\end{array}\end{array}$$

Substituting Eqs. (A1) to (A3) into Eq. (A5) and differentiating with respect to θ yields

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}{C}_{\mathrm{pd}}\mathit{\theta }{\displaystyle \frac{\partial }{\partial t}}& {\displaystyle \frac{\partial M_{\mathit{\theta }}(\mathit{\theta },t)}{\partial \mathit{\theta }}}=C_{\mathrm{pd}}\mathit{\theta }\left({\displaystyle \frac{\partial F_{\mathrm{T}}\left(\mathit{\theta },t\right)}{\partial \mathit{\theta }}}+{\displaystyle \frac{\partial F_{\mathrm{E}}\left(\mathit{\theta },t\right)}{\partial \mathit{\theta }}}\right.\\ & \left.+{\displaystyle \frac{\partial F_{\mathrm{I}}\left(\mathit{\theta },t\right)}{\partial \mathit{\theta }}}\right)+C_{\mathrm{pd}}{F}_{\mathrm{I}}(\mathit{\theta },t)+{\displaystyle \frac{\partial Q_{\mathrm{diff}}(\mathit{\theta },t)}{\partial \mathit{\theta }}}\\ & +{\displaystyle \frac{\partial Q_{\mathrm{s}}(\mathit{\theta },t)}{\partial \mathit{\theta }}}+{\displaystyle \frac{\partial Q_{\mathrm{int}}(\mathit{\theta },t)}{\partial \mathit{\theta }}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{Q}_{\mathrm{diff}}(\mathit{\theta },t)& =\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial D_{\mathrm{T}}\left({\mathit{\theta }}^{\prime },t\right)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime }+\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial D_{\mathrm{E}}\left({\mathit{\theta }}^{\prime },t\right)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime }\\ & +D_{\mathrm{I}}(\mathit{\theta },t).\end{array}\end{array}$$

Differentiating Eq. (A4) with respect to θ and multiplying C_pd⋅θ yields

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}{C}_{\mathrm{pd}}\mathit{\theta }{\displaystyle \frac{\partial }{\partial t}}{\displaystyle \frac{\partial M_{\mathit{\theta }}(\mathit{\theta },t)}{\partial \mathit{\theta }}}& =C_{\mathrm{pd}}\mathit{\theta }\left({\displaystyle \frac{\partial F_{\mathrm{T}}\left(\mathit{\theta },t\right)}{\partial \mathit{\theta }}}+{\displaystyle \frac{\partial F_{\mathrm{E}}\left(\mathit{\theta },t\right)}{\partial \mathit{\theta }}}\right.\\ & \left.+{\displaystyle \frac{\partial F_{\mathrm{I}}\left(\mathit{\theta },t\right)}{\partial \mathit{\theta }}}\right).\end{array}\end{array}$$

Subtracting Eq. (A8) from Eq. (A6), we obtain

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{C}_{\mathrm{pd}}{F}_{\mathrm{I}}(\mathit{\theta },t)& =-{\displaystyle \frac{\partial Q_{\mathrm{diff}}(\mathit{\theta },t)}{\partial \mathit{\theta }}}-{\displaystyle \frac{\partial Q_{\mathrm{s}}(\mathit{\theta },t)}{\partial \mathit{\theta }}}\\ & -{\displaystyle \frac{\partial Q_{\mathrm{int}}(\mathit{\theta },t)}{\partial \mathit{\theta }}}.\end{array}\end{array}$$

Equation (A9) divided by C_pd plus Eq. (A4) yields

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}{M}_{\mathit{\theta }}(\mathit{\theta },t)& =-{\displaystyle \frac{\mathrm{1}}{C_{\mathrm{pd}}}}\left({\displaystyle \frac{\partial Q_{\mathrm{diff}}(\mathit{\theta },t)}{\partial \mathit{\theta }}}+{\displaystyle \frac{\partial Q_{\mathrm{s}}(\mathit{\theta },t)}{\partial \mathit{\theta }}}\right.\\ & \left.+{\displaystyle \frac{\partial Q_{\mathrm{int}}(\mathit{\theta },t)}{\partial \mathit{\theta }}}\right)+\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial F_{\mathrm{T}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime }\\ & +\underset{-\mathrm{\infty }}{\overset{\mathit{\theta }}{\int }}{\displaystyle \frac{\partial F_{\mathrm{E}}({\mathit{\theta }}^{\prime },t)}{\partial {\mathit{\theta }}^{\prime }}}\mathrm{d}{\mathit{\theta }}^{\prime }.\end{array}\end{array}$$

Equation (A10) illustrates the temporal variation in M_θ, where Q_int includes radiative heating (i.e., sum of shortwave and longwave heating), dissipation of the kinetic energy of turbulence, and latent heat release due to evaporation and condensation.

To modify Eq. (A10) to apply to $M_{{\mathit{\theta }}_{\mathrm{e}}}$ rather than to M_θ, it is necessary to replace all θ with θ_e and additionally account for the following:

1. Condensation and evaporation is conserved on the θ_e surfaces, but the gaining and losing of water vapor through surface evaporation and water vapor transport contributes to θ_e. This contribution can be computed as the product of latent heat of evaporation and the extra water vapor content. Thus, the surface contribution (Q_S) needs to include both sensible heating of the atmosphere (Q_sen) and the water vapor flux from the surface into the atmosphere (Q_evap). Similarly, the diffusion term within the atmosphere (Q_diff) needs to include both heat and water vapor ($Q_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$).

2. Internal heating (Q_int) needs to exclude latent heat releasing due to evaporation and condensation of liquid water, which cancel in θ_e, but it still needs to include heating from ice formation, which does not cancel in θ_e. We subtract this ice component from the rest of the internal heating, yielding two terms $Q_{\mathrm{int}}^{\prime }$ and Q_ice, with $Q_{\mathrm{int}}=Q_{\mathrm{int}}^{\prime }+Q_{\mathrm{ice}}$.

Therefore, we can write the temporal variation in $M_{{\mathit{\theta }}_{\mathrm{e}}}$ as

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}{M}_{{\mathit{\theta }}_{\mathrm{e}}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)& =\underset{-\mathrm{\infty }}{\overset{{\mathit{\theta }}_{\mathrm{e}}}{\int }}{\displaystyle \frac{\partial F_{\mathrm{T}}\left({\mathit{\theta }}_{\mathrm{e}}^{\prime },t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}^{\prime }}}\mathrm{d}{\mathit{\theta }}_{\mathrm{e}}^{\prime }+\underset{-\mathrm{\infty }}{\overset{{\mathit{\theta }}_{\mathrm{e}}}{\int }}{\displaystyle \frac{\partial F_{\mathrm{E}}\left({\mathit{\theta }}_{\mathrm{e}}^{\prime },t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}^{\prime }}}\mathrm{d}{\mathit{\theta }}_{\mathrm{e}}^{\prime }\\ & -{\displaystyle \frac{\mathrm{1}}{C_{\mathrm{pd}}}}\left({\displaystyle \frac{\partial Q_{\mathrm{diff}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}}}+{\displaystyle \frac{\partial Q_{\mathrm{sen}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}}}\right.\\ & +{\displaystyle \frac{\partial Q_{\mathrm{evap}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}}}+{\displaystyle \frac{\partial Q_{\mathrm{int}}^{\prime }\left({\mathit{\theta }}_{\mathrm{e}},t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}}}\\ & \left.+{\displaystyle \frac{\partial Q_{\mathrm{ice}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}}}+{\displaystyle \frac{\partial Q_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left({\mathit{\theta }}_{\mathrm{e}},t\right)}{\partial {\mathit{\theta }}_{\mathrm{e}}}}\right).\end{array}\end{array}$$