Impact of precipitation mass sinks on midlatitude storms in idealized simulations across a wide range of climates

The combination of precipitation formation and fallout affects atmospheric flows through the release of latent heat and through the removal of mass from the atmosphere, but because the mass of water vapor is only a small fraction of the total mass of Earth's atmosphere, precipitation mass sinks are often neglected in theory and models. However, a small number of modeling studies suggest that water mass sources and sinks can intensify heavily precipitating weather systems. These studies point to a need to more systematically verify the impact of neglecting precipitation mass sinks, particularly for warmer and moister climates in which precipitation rates can be much higher. In this paper, we add precipitation mass sources and sinks to an idealized general circulation model and examine their effects on steady-state midlatitude storm track statistics. The model has several idealizations, including that all condensates immediately fall out of the atmosphere, and is run across a wide range of climates, including very warm climates. We find that modifying the model to include mass sources and sinks has no detectable effect on midlatitude variability or extremes, even in climates much warmer and moister than the modern. However, we find that a 10-fold exaggeration of mass sources and sinks is sufficient to produce more intense midlatitude weather extremes and increase surface pressure variance. This result is consistent with theoretical potential vorticity analysis, which suggests that the dynamical effects of mass sources and sinks are much smaller than the dynamical effects of accompanying latent heating unless mass sinks are artificially amplified by at least a factor of 10. Finally, we use simulations of “tropical cyclone worlds” to attempt to reconcile our results with earlier work showing stronger deepening in a simulation of a tropical cyclone case study when precipitation mass sinks were included. We demonstrate that abruptly “turning on” mass sources and sinks can lead to stronger transient deepening in some individual storms (consistent with results of past work) but weaker transient deepening in other storms, without modifying the steady-state statistics of storms in equilibrium with the large-scale environment (consistent with our other results). Our results provide a firmer foundation for using general circulation models that neglect moist mass sources and sinks in climate simulations, even in climates much warmer than today, while leaving open the possibility that their inclusion might lead to short-term improvements in forecast skill.

1 Introduction

In midlatitudes, large-scale atmospheric flows are dominated by extratropical cyclones and anticyclones. These disturbances play an important role in climate and the general circulation of the atmosphere through transports of heat, moisture, and angular momentum (Held, 1975; Pierrehumbert and Swanson, 1995; Chang et al., 2002). Additionally, they organize extratropical precipitation (Pfahl and Wernli, 2012) and are the primary source of day-to-day weather variability in midlatitudes (Peixoto and Oort, 1992) and, consequently, have significant impacts on society. Understanding the strength of midlatitude cyclones requires understanding the mechanisms that drive their growth and decay, and much of what we know is based on theory for dry atmospheres (e.g., Charney, 1947; Eady, 1949; Lorenz, 1955; Salmon, 1980; O'Gorman and Schneider, 2008a; O'Gorman, 2010). However, Earth's atmosphere is not dry, and the presence of moisture also affects the dynamics of midlatitude weather systems.

One way moisture influences midlatitude cyclones is through latent heat release by condensation and precipitation. Condensation can enhance lower-level cyclonic potential vorticity (PV) anomalies and surface circulations and, in some cases, entirely replace the effect of low-level advection of PV across a meridional gradient (Kuo et al., 1991; Snyder and Lindzen, 1991; Davis and Emanuel, 1991; Davis, 1992; Stoelinga, 1996; Ahmadi-Givi et al., 2004; Joos and Wernli, 2012; Tamarin and Kaspi, 2016). Additionally, latent heating can erode upper-level cyclonic PV slightly downstream of the peak condensation and, in doing so, contribute to the westward phase tilt required for mutual Rossby wave amplification (Grams et al., 2011; Schemm et al., 2013). These effects matter: Emanuel et al. (1987) showed that baroclinic growth rates in a two-dimensional semigeostrophic model more than doubled when latent heating was included, and the growth rate can be even higher for diabatic Rossby vortices in which latent heating is central to the dynamics (Kohl and O'Gorman, 2022). O'Gorman (2011) showed how a reduced effective static stability due to latent heating affects a number of aspects of the general circulation.

Although the importance of latent heating for midlatitude dynamics is well established, many atmospheric models neglect a different moisture effect that plays a role in dynamics and mass budgets: mass transports by hydrometeors (Trenberth, 1991; Trenberth et al., 1995; Qiu et al., 1993; Lackmann and Yablonsky, 2004; Wacker et al., 2006). Precipitation mass sinks reduce hydrostatic surface pressure (Spengler et al., 2011), and this effect will be counteracted to some extent (depending on length scale) by convergent flow, leading to vortex stretching. Like latent heating, falling hydrometeors act as sources and sinks of PV by moving mass across isentropic surfaces (Schubert et al., 2001; Lackmann and Yablonsky, 2004). However, while latent heat comprises about 2.5 % of the total energy of the atmosphere (Peixoto and Oort, 1992), water vapor contributes only about 0.25 % of the total mass of the atmosphere (Trenberth and Smith, 2005), and changes in atmospheric mass from water vapor sources and sinks are ignored in many meteorological applications.

Nonetheless, some modeling studies have suggested that water vapor mass sources and sinks play a non-negligible role in the dynamics of heavily precipitating weather systems. In simulations of an extratropical cyclone with a mesoscale model, Qiu et al. (1993) found that including precipitation mass sinks increased peak precipitation rates by about 10 %. A similar effect was found for simulations of Hurricane Lili (2002): when mass sinks associated with precipitation were included, peak 10 m wind speeds increased by between 5 and 15 kn, and the storm deepened by an additional 2–5 hPa (Lackmann and Yablonsky, 2004, LY04 hereafter). In the latter study, LY04 compared the production of PV by latent heating and mass sinks and found that, although local PV tendencies produced by latent heating were much larger than those produced by mass sinks, their vertical structures differed. Latent heating in Lili produced PV tendencies with vertically variable signs; mass sinks, on the other hand, produced single-signed PV increases. As a consequence, spatially integrated PV tendencies from the mass sink were non-negligible compared to those from latent heating, and LY04 argued that cancellation between opposite-signed PV tendencies produced by latent heating might allow the mass sink to play a comparatively non-negligible role in generating balanced circulations. Having both positive and negative PV tendencies from latent heating may decrease the effect of latent heating on storm strength either through transport, which can bring diabatically generated PV of one sign into a region of generation of the opposite sign (e.g., Schubert and Alworth, 1987; Büeler and Pfahl, 2017), or through inversion, which implies that winds are related to spatially smoothed PV anomalies (Hoskins et al., 1985).

Given that midlatitude baroclinic cyclones are significantly impacted by moisture, and in light of evidence that the precipitation mass sink may be important for the dynamics of some heavily precipitating systems, this paper investigates whether sources and sinks of water mass play a quantitatively significant role in weather variability and extremes in midlatitude storm tracks. Our work builds on Qiu et al. (1993) and LY04 but makes a number of novel contributions. First, we examine the importance of the mass sink both in climates similar to modern Earth and in much warmer and moister climates where mass sinks due to condensation and precipitation become significantly larger. Second, we focus primarily on the effect of mass sources and sinks on statistics of many storms in statistical equilibrium, not on the transient evolution of individual weather systems. Third, we develop theoretical results that provide insight into how PV sources and sinks, vertical transport, and inversion combine to influence the relative importance of latent heating versus mass sources and sinks. Finally, we consider the effects of including mass sources and sinks in tropical cyclone (TC) world simulations, providing a single-model comparison of the relative importance of mass sources and sinks for tropical and extratropical storms.

We proceed in three parts. First (Sects. 2 and 3), we examine statistical measures of midlatitude variability and extremes in simulations across a broad range of climates using a widely used idealized atmospheric general circulation model (Frierson et al., 2006, 2007; O'Gorman and Schneider, 2008b) that we modify to include sources and sinks of water mass. We find little difference between simulations with and without mass sources and sinks, even in climates much warmer than modern Earth, but do observe changes indicative of an increase in the strength of the strongest midlatitude cyclones when we exaggerate mass sources and sinks by a factor of 10. In Sect. 4, we connect simulation results with arguments based on isentropic potential vorticity, with a focus on explaining why effects from mass sources and sinks become apparent when mass sinks are exaggerated by a factor of 10 (and not, e.g., by a factor of 1 or 100). Finally, we attempt to reconcile our results with LY04 by examining the effect of mass sources and sinks on the intensity of TCs in TC world simulations in the same idealized general circulation model (GCM) (Sect. 5). We find that adding mass sources and sinks can lead to transient deepening of some individual simulated TCs – similar to the results of LY04 – but other TCs weaken, and systematic effects on statistics of large populations of TCs are relatively weak – similar to our results for midlatitude storm tracks.

2 Simulation design

Our simulations use an idealized moist general circulation model based on the Geophysical Fluid Dynamics Laboratory (GFDL) spectral dynamical core. The model is documented in Frierson et al. (2006) and Frierson et al. (2007), with details largely following O'Gorman and Schneider (2008b); we briefly note some important aspects here.

The model uses two-stream gray radiative transfer with the insolation and longwave and shortwave absorber distributions described in O'Gorman and Schneider (2008b). We access a wide range of climates by varying a parameter α that multiplies the longwave optical thickness distribution from α=1, which gives a global-average surface temperature close to that of Earth today, to α=6, which gives a global-average surface temperature of around 315 K. Surface fluxes of momentum and sensible and latent heat are calculated using standard drag laws. The lower boundary is a mixed layer of 1 m depth whose temperature changes in response to radiation and surface turbulent fluxes. The dynamical core uses a spectral solver with ∇⁸ hyperdiffusion of vorticity, divergence, and temperature.

Condensation and precipitation are produced by a simplified Betts–Miller convection scheme (Frierson, 2007) and a large-scale condensation scheme. The convection scheme relaxes temperatures toward a moist adiabat and moisture toward 70 % relative humidity on a 2 h timescale and ensures enthalpy conservation by modifying the timescale for specific humidity adjustment. The large-scale condensation scheme ensures that grid-scale relative humidity cannot exceed 100 %. Unlike O'Gorman and Schneider (2008b), we include reevaporation of falling condensate, which is assumed to bring subsaturated grid cells to saturation before falling through them. The convection and large-scale condensation schemes both assume a pseudoadiabatic limit in which water is removed from the atmosphere as soon at it condenses.

Our main modification to the model is the inclusion of mass sources and sinks associated with evaporation and condensate fallout, respectively. The changes are described in detail in Appendix A, but we summarize them here. Briefly, we add terms to equations for pressure tendencies at the surface and at interfaces between levels to represent the effects of adding and removing mass from the overlying column (Sect. A1). Our implementation scales pressure tendencies from mass sources and sinks by a parameter γ. Setting γ=0 corresponds to neglecting pressure tendencies from mass sources and sinks, setting γ=1 corresponds to including the correct pressure tendencies from mass sources and sinks, and setting γ>1 exaggerates pressure tendencies from mass sources and sinks. Note that γ modifies pressure tendencies from mass sources and sinks but does not alter humidity tendencies themselves. This means that setting γ=0 implies that sinks of moist mass are exactly offset by sources of dry mass and that setting γ>1 implies that sinks of moist mass are accompanied by additional sinks of dry mass. Numerical errors in our implementation of the mass sinks lead to a slow downward drift in total mass, and so our simulations use a dry-mass fixer that ensures conservation of a fixed dry-mass inventory of 10⁵ Pa (Sect. A3). For simplicity, we refer to the dry-mass inventory as “dry mass” and the additional atmospheric mass due to water vapor as “moist mass” for the remainder of the paper. We do not include energy or moisture fixers in our simulations, as neither are required to achieve statistical steady states. Finally, we include minor changes to the calculation of the reference moist adiabat used by the convection scheme to improve consistency with the model dynamical core (Sect. A4), but these are unrelated to the addition of mass sources and sinks.

The form of the mass sink term used in our simulations relies on the pseudoadiabatic assumption, which is consistent with the model physics but nonetheless approximate. In the pseudoadiabatic limit, water mass disappears immediately upon condensation, and the resulting mass loss has an immediate impact on the hydrostatic pressure field at all levels below the level where condensation occurs. An additional consequence of our use of a pseudoadiabatic model is that our simulations contain no representation of the effects of condensate loading on the hydrostatic pressure field. In reality, mass sinks are not associated with condensation itself but with subsequent hydrometeor fallout (e.g., Ooyama, 2001), and condensate loading from cloud water and sedimenting hydrometeors contributes to the total hydrostatic pressure. As shown in Appendix B, the pseudoadiabatic mass sink term can be derived from a mass sink term that includes condensate loading and explicit hydrometeor fallout by assuming that timescales for precipitation formation and fallout are infinitely fast. Because our study focuses primarily on the dynamics of midlatitude eddies and tropical cyclones (large-scale weather systems that evolve over timescales of days), we think that treating precipitation formation and fallout as infinitely fast is a reasonable simplifying assumption. Our use of the pseudoadiabatic limit is likely to affect some aspects of our simulations (for example, the tropics may be nearly neutral to pseudoadiabatic ascent, rather than to reversible adiabatic ascent as found in Xu and Emanuel, 1989), but we do not think that the use of a full cloud microphysics scheme or prognostic equations for condensed water species are essential for our study, which focuses on the existence or absence of the mass sink itself. Nevertheless, we would welcome follow-up studies that examine whether a more sophisticated treatment of clouds and precipitation alters any of our conclusions.

We use a horizontal resolution of T170 (nominal grid spacing of about 80 km at the Equator) and 30 unevenly spaced sigma levels for all simulations, and we access a wide range of climates by running sets of simulations with $\mathit{\alpha }=(\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5},\mathrm{6})$. All of our simulations are spun up from an isothermal rest state for 700 d, and statistics are calculated over the following 300 d based on instantaneous snapshots taken every 6 h. We refer to α=1 simulations, which produce global-mean surface air temperatures close to modern values, as “control” simulations, and simulations with α>1 as “warm” simulations. At each α, we run simulations with three difference treatments of mass sources and sinks. One treatment neglects mass sources and sinks (γ=0, “no MS”), one includes mass sources and sinks (γ=1, “MS”), and one exaggerates mass sources and sinks by a factor of 10 (γ=10, “10× MS”). Total atmospheric mass remains equal to the dry-mass inventory in all no-MS simulations, increases by as much as 3 % in warm MS simulations, and increases by between about 3 % and 25 % (depending on temperature) in 10×-MS simulations (Fig. 1). The mass increases in the MS and 10×-MS simulations because of the source of mass from evaporation, which comes into balance with the precipitation mass sink at statistical steady state, and this evaporation mass source is larger again in the 10×-MS simulations.

Figure 1 Ratio of moist atmospheric mass to dry atmospheric mass, plotted as a function of global-mean air temperature in the lowest model level. The horizontal gray line is a visual aid for comparing moist mass in the warmest MS simulation and control 10×-MS simulation.

3 Impact of the mass sink on midlatitude storm tracks

Our simulations, in both control and warm climates and for all treatments of the mass sink, produce midlatitude storm tracks with pronounced eddy activity and relatively quiescent tropics (Fig. 2). The character of the eddy field, as indicated by near-surface wind speeds, changes significantly between the control climate (α=1, Fig. 2a–c) and the warmest climate (α=6, Fig. 2d–f). In the control climate, strong near-surface winds are widespread across midlatitude baroclinic zones. In the warmest climate, eddies appear less vigorous: near-surface winds are strong only in localized eddies and relatively weak elsewhere. Generally weaker near-surface winds in warmer climates are a symptom of reduced midlatitude eddy kinetic energy, which decreases in the warmest climates because weaker meridional temperature gradients and higher midlatitude dry static stability reduce mean available potential energy (O'Gorman and Schneider, 2008a). Regions of strong near-surface winds may appear more localized in warmer climates because of changes to the dominant mode of midlatitude instability (O'Gorman et al., 2018; Kohl and O'Gorman, 2022). In control simulations, precipitation is strongest in frontal zones, where precipitation rates of 1–3 mm h⁻¹ are typical (Fig. 3a–c). In the warmest simulations, precipitation rates are highest near the center of cyclones and often reach values above 10 mm h⁻¹ (Fig. 3d–f).

We focus on statistics calculated over midlatitude baroclinic zones, defined following O'Gorman and Schneider (2008a) as regions within 15^∘ latitude of the maximum of the vertically integrated eddy potential temperature flux $\stackrel{\mathrm{‾}}{v^{\prime }{\mathit{\theta }}^{\prime }}\cos \mathit{\varphi }$, where v is meridional velocity, θ is potential temperature, ϕ is latitude, $(\cdot {)}^{\prime }$ is an eddy quantity defined relative to the surface-pressure-weighted zonal and time mean $\stackrel{\mathrm{‾}}{(\cdot )}$ along σ surfaces, and the vertical integral is taken from the surface to the lowest level where the zonal- and time-mean temperature lapse rate drops below 2 K km⁻¹. The baroclinic zones track regions of strong eddy activity and shift poleward in warmer climates (Fig. 2). We compute 95 % confidence intervals for all statistics by resampling raw model output in time using block bootstrapping with a block size of 10 d. We report confidence intervals based on 200 block bootstrap realizations, each with 300 d of data. The confidence intervals are robust to changes in the number of bootstrap realizations; using 100 or 50 bootstrap realizations modifies them only slightly.

Figure 2 Wind speed in the lowest model level from the last saved snapshot of simulations in the control climate with α=1 (a–c) and the warmest climate with α=6 (d–f) with no mass sinks (a, d), realistic mass sinks (b, e), and mass sinks exaggerated by a factor of 10 (c, f). Midlatitude baroclinic zones defined based on the meridional eddy potential temperature flux are delimited by dashed lines, and rectangles delimit regions for which precipitation rates are shown in Fig. 3.

Figure 3 Surface precipitation rates from the last saved snapshot of simulations in the control climate with α=1 (a–c) and the warmest climate with α=6 (d–f) with no mass sinks (a, d), realistic mass sinks (b, e), and mass sinks exaggerated by a factor of 10 (c, f). Regions shown are outlined with rectangles in Fig. 2.

Figure 4 shows vertically averaged eddy kinetic energy (EKE), near-surface EKE, and surface pressure variance averaged over the midlatitude baroclinic zones. Vertically averaged EKE is averaged over all model levels, and near-surface EKE is defined as EKE in the lowest model level. Vertically integrated and near-surface EKE both decrease with warming (O'Gorman and Schneider, 2008a), and values from simulations with different treatments of the mass sink are statistically indistinguishable. Like EKE, surface pressure variance decreases with warming and is statistically indistinguishable in no-MS and MS simulations. However, in 10×-MS simulations, surface pressure variance increases with warming and is substantially higher than in no-MS and MS simulations at α=3 and above (global-mean surface temperatures warmer than ∼305 K). No-MS and 10×-MS simulations can differ in surface pressure variance even though they do not differ in EKE because large mass sinks produce variations in surface pressure with an unbalanced component that is not associated with large-scale fluid flow.

Figure 4 Vertically integrated EKE (a), near-surface EKE (b), and surface pressure variance (c) averaged over midlatitude baroclinic zones and plotted as a function of global-mean surface air temperature. Vertical lines indicate 95 % confidence intervals computed using block bootstrap resampling as described in the text. Horizontal lines indicating 95 % confidence intervals for global-mean surface temperature are also included but are too small to be visible.

Figure 5 shows extreme percentiles of near-surface wind speeds, surface pressure anomalies, and (instantaneous) precipitation in midlatitude baroclinic zones. In all cases, percentiles are calculated as a function of latitude (aggregating over longitude and time), then combined in area-weighted latitudinal averages over the baroclinic zones. Precipitation percentiles are calculated for all times, not just times when rain occurs. The 99th percentile near-surface wind speeds, defined as wind speeds in the lowest model level, decrease with warming and are statistically indistinguishable in no-MS and MS simulations but are systematically higher in 10×-MS simulations. Extreme surface lows, calculated as the first percentile of surface pressure anomalies relative to zonal and time averages, are likewise statistically indistinguishable in no-MS and MS simulations but systematically deeper in 10×-MS simulations in warm climates. The behavior of 99.9th percentile rain rates mirrors that of extreme wind speeds and surface lows: no-MS and MS simulations are again statistically indistinguishable, but 10×-MS simulations produce systematically stronger precipitation in warm climates. The 99.9th percentile precipitation rates increase with warming, as expected from simple arguments about the influence of warming on extreme precipitation (O'Gorman and Schneider, 2009; Muller et al., 2011). Lower precipitation percentiles (e.g., 99th percentile precipitation rates) decrease with warming in the warmest climates (not shown), suggesting that increases in 99th percentile precipitation rates may be limited to some extent by energetic constraints on mean precipitation (O'Gorman and Schneider, 2008b; Pendergrass and Hartmann, 2014a, b) or by decreases in mean midlatitude eddy kinetic energy (O'Gorman and Schneider, 2008a).

Figure 5 Extreme near-surface wind speeds (a), surface pressure anomalies (b), and precipitation (c) calculated over midlatitude baroclinic zones and plotted as a function of global-mean surface air temperature. The 99th and 99.9th percentiles in panels (a) and (c) highlight particularly strong near-surface wind and precipitation; the 1st percentiles shown in panel (b) highlight particularly strong surface lows. Vertical lines indicate 95 % confidence intervals computed using block bootstrap resampling as described in the text. Horizontal lines indicating 95 % confidence intervals for global-mean surface temperature are also included but are too small to be visible.

In summary, statistics calculated over midlatitude baroclinic zones are consistently statistically indistinguishable in no-MS and MS simulations but noticeably different between no-MS and 10×-MS simulations, particularly in warm climates. Differences between no-MS and 10×-MS simulations suggest that extreme disturbances and the unbalanced component of surface pressure anomalies become more intense with a 10-fold-exaggerated mass sink: midlatitude surface pressure variance increases, extreme near-surface winds become faster, extreme surface lows become deeper, and extreme precipitation strengthens. In contrast to other statistics, mean midlatitude EKE does not show a statistically significant difference between no-MS and 10×-MS simulations, possibly because of constraints imposed by the mean available potential energy that apply to the mean EKE of baroclinic zones (O'Gorman and Schneider, 2008a) but not to individual weather systems or to unbalanced pressure variations. The overall weight of evidence points to minor differences in storm track variability and extremes between no-MS and MS simulations but noticeable increases in the strength of extreme disturbances and surface pressure variability in 10×-MS simulations.

Our results, which show no detectable differences between no-MS and MS simulations even in climates much warmer than the modern, are seemingly inconsistent with the results of LY04, which document a clear increase in TC intensity in a simulation of Hurricane Lili with mass sinks included. This inconsistency, plus the fact that we do find differences between no-MS and 10×-MS simulations, prompts us to better understand the underlying theory and ask the following question: is there a reason to expect the dynamical effects of unexaggerated mass sinks to be dominated by other processes but to become at least marginally detectable when exaggerated by a factor of 10? In the next section, we use some simple theoretical analysis of isentropic potential vorticity equations to argue that there is.

4 Isentropic potential vorticity theory

In this section, we analyze isentropic potential vorticity (IPV) equations to probe why flows in midlatitude baroclinic zones are insensitive to the inclusion of mass sources and sinks but affected by 10-fold-exaggerated mass sources and sinks. Following LY04, our approach is to compare the effects of mass sinks with the effects of latent heating, the latter being a convenient benchmark that generally is considered to be important to include in weather and climate models. We first show that the magnitude of IPV sources produced by latent heating is larger than that produced by mass sinks by a factor of about 200 (Sect. 4.1). This factor is influenced by thermodynamic constants and tropospheric stratification but is not directly climate dependent because latent heating and mass sinks increase in tandem with warming. We then analyze the evolution of IPV anomalies produced by sustained latent heating with a first-baroclinic-mode vertical structure, which produces an IPV source in the lower troposphere and an IPV sink in the upper troposphere. We show that the combined effects of latent heating and vertical advection produce, at long timescales, a positive IPV anomaly throughout most of the depth of the atmosphere and that mass sinks increase the positive IPV anomaly only by about 1 part in 200 (Sect. 4.2). At short timescales, cancellation between opposite-signed IPV anomalies produced by latent heating can limit the strength of induced flow because the inversion is non-local, but we show that flows induced by latent heating remain stronger than flows induced by mass sinks by a factor of 100 or more, except for disturbances with horizontal length scales much larger than a deformation radius (Sect. 4.3).

For consistency with our simulations, our IPV analysis assumes a pseudoadiabatic limit where condensation falls out of the atmosphere immediately upon formation, which influences the form in which mass sinks appear as IPV source terms. As discussed previously, we view our use of a pseudoadiabatic limit as a reasonable choice given the focus of our work. However, we include in Appendix D a derivation of an IPV equation that does not assume pseudoadiabatic fallout (in which case the IPV source from mass sinks involves the divergence of the hydrometeor sedimentation mass flux), and we show that the pseudoadiabatic IPV equation used in this section can be recovered by assuming that precipitation formation and fallout are infinitely fast.

4.1 IPV tendencies from latent heating and mass sinks

We begin with frictionless momentum and continuity equations in potential temperature coordinates on an f-plane (Holton and Hakim, 2012):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathbf{v}}{\partial t}}+{\mathrm{\nabla }}_{\mathit{\theta }}\left({\displaystyle \frac{\mathbf{v}\cdot \mathbf{v}}{\mathrm{2}}}+M\right)+({\mathit{\zeta }}_{\mathit{\theta }}+f)\mathbf{k}\times \mathbf{v}=-\dot{\mathit{\theta }}{\displaystyle \frac{\partial \mathbf{v}}{\partial \mathit{\theta }}},\text{(2)}& {\displaystyle \frac{\partial \mathit{\sigma }}{\partial t}}+{\mathrm{\nabla }}_{\mathit{\theta }}\cdot \left(\mathit{\sigma }\mathbf{v}\right)+{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\mathit{\sigma }\dot{\mathit{\theta }}\right)=\mathit{\sigma }\dot{q_{\mathrm{v}}}.\end{array}$$

θ is potential temperature, v is the horizontal velocity, ∇_θ is the horizontal gradient on potential temperature surfaces, ${\mathit{\zeta }}_{\mathit{\theta }}=\mathbf{k}\cdot {\mathrm{\nabla }}_{\mathit{\theta }}\times \mathbf{v}$ is the relative vorticity, and f is the Coriolis parameter. $M=c_{p}T+gz$ is the Montgomery streamfunction, with c_p being the specific heat capacity of dry air, T the absolute temperature, g the acceleration from gravity, and z geometric height. $\mathit{\sigma }=-g^{-\mathrm{1}}\partial p/\partial \mathit{\theta }$ is the isentropic density, with p being the pressure. $\dot{\mathit{\theta }}$ is a diabatic potential temperature tendency, used here to represent the effect of latent heat release, and $\dot{q_{\mathrm{v}}}$ is the specific humidity tendency from condensation and precipitation fallout. We consider a pseudoadiabatic limit where all condensation is immediately removed from the atmosphere, consistent with model simulations described earlier, so $\dot{q_{\mathrm{v}}}$ corresponds directly to a mass sink. A derivation of the continuity equation including mass sources and sinks is provided in Appendix C; the version used here approximates $\mathrm{1}-q_{\mathrm{v}}\approx \mathrm{1}$.

Following Holton and Hakim (2012), we take $\mathbf{k}\cdot {\mathrm{\nabla }}_{\mathit{\theta }}\times$ of the momentum equation to obtain

$$\begin{array}{}\text{(3)}& \begin{array}{rl}({\displaystyle \frac{\partial }{\partial t}}& +\mathbf{v}\cdot {\mathrm{\nabla }}_{\mathit{\theta }})\left({\mathit{\zeta }}_{\mathit{\theta }}+f\right)+\left({\mathit{\zeta }}_{\mathit{\theta }}+f\right){\mathrm{\nabla }}_{\mathit{\theta }}\cdot \mathbf{v}\\ & =-\mathbf{k}\cdot {\mathrm{\nabla }}_{\mathit{\theta }}\times \left(\dot{\mathit{\theta }}{\displaystyle \frac{\partial \mathbf{v}}{\partial \mathit{\theta }}}\right),\end{array}\end{array}$$

then expand the right-hand-side derivative to rewrite the result as

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\left({\mathit{\zeta }}_{\mathit{\theta }}+f\right)+\left({\mathit{\zeta }}_{\mathit{\theta }}+f\right){\mathrm{\nabla }}_{\mathit{\theta }}\cdot \mathbf{v}=\mathbf{k}\cdot \left({\displaystyle \frac{\partial \mathbf{v}}{\partial \mathit{\theta }}}\times {\mathrm{\nabla }}_{\mathit{\theta }}\dot{\mathit{\theta }}\right),\end{array}$$

where

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}={\displaystyle \frac{\partial }{\partial t}}+\mathbf{v}\cdot {\mathrm{\nabla }}_{\mathit{\theta }}+\dot{\mathit{\theta }}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\end{array}$$

is the material derivative expressed in potential temperature coordinates. Using the continuity equation to eliminate the stretching term gives

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{D}P}{\mathrm{D}t}}=P\left({\displaystyle \frac{\partial \dot{\mathit{\theta }}}{\partial \mathit{\theta }}}-\dot{q_{\mathrm{v}}}\right)+{\mathit{\sigma }}^{-\mathrm{1}}\mathbf{k}\cdot \left({\displaystyle \frac{\partial \mathbf{v}}{\partial \mathit{\theta }}}\times {\mathrm{\nabla }}_{\mathit{\theta }}\dot{\mathit{\theta }}\right),\end{array}$$

where

$$\begin{array}{}\text{(7)}& P={\displaystyle \frac{{\mathit{\zeta }}_{\mathit{\theta }}+f}{\mathit{\sigma }}}\end{array}$$

is the IPV. The source terms involving $\dot{\mathit{\theta }}$ and $\dot{q_{\mathrm{v}}}$ are similar to source terms for Ertel potential vorticity (EPV) derived by Schubert et al. (2001), and the IPV evolution equation is identical to that given in Hoskins (1991) except that we assume the flow is frictionless and include an IPV source term for mass sources and sinks.

Next, we assume we are at a horizontal maximum in latent heating such that the final term in Eq. (6) vanishes. We further assume that the latent heating profile is given by

$$\begin{array}{}\text{(8)}& \dot{\mathit{\theta }}={\displaystyle \frac{{\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}}{T}}\sin \left(\mathit{\pi }{\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{s}}}{{\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}}}\right),\end{array}$$

where θ_s and θ_t are the potential temperature at the surface and the top of the region of latent heating, and T is a timescale. In the pseudoadiabatic limit, the precipitation mass sink is related to the latent heating rate by

$$\begin{array}{}\text{(9)}& \dot{q_{\mathrm{v}}}=-{\displaystyle \frac{\mathrm{\Pi }}{L_{\mathrm{v}}}}\dot{\mathit{\theta }},\end{array}$$

where L_v is the latent heat of vaporization for water and

$$\begin{array}{}\text{(10)}& \mathrm{\Pi }={\displaystyle \frac{c_{p}T}{\mathit{\theta }}}=c_p{\left({\displaystyle \frac{p}{p_{\mathrm{s}}}}\right)}^{R/c_p}\end{array}$$

is the Exner function, with p_s being the surface pressure and R the gas constant for dry air. Substituting these expressions for $\dot{\mathit{\theta }}$ and $\dot{q_{\mathrm{v}}}$ into Eq. (6) and non-dimensionalizing by setting

$$\begin{array}{ll}{\displaystyle }t& {\displaystyle }\to T\stackrel{\mathrm{̃}}{t}\\ {\displaystyle }\mathit{\theta }-{\mathit{\theta }}_{\mathrm{s}}& {\displaystyle }\to \left({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\right)\stackrel{\mathrm{̃}}{\mathit{\theta }}\\ {\displaystyle }P& {\displaystyle }\to {\displaystyle \frac{f}{{\mathit{\sigma }}_{\mathrm{0}}}}\stackrel{\mathrm{̃}}{P}\\ {\displaystyle }\mathrm{\Pi }& {\displaystyle }\to c_p\stackrel{\mathrm{̃}}{\mathrm{\Pi }},\end{array}$$

where σ₀ is a reference isentropic density, gives

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{D}\stackrel{\mathrm{̃}}{P}}{\mathrm{D}\stackrel{\mathrm{̃}}{t}}}=\stackrel{\mathrm{̃}}{P}\left[\mathit{\pi }\cos \left(\mathit{\pi }\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)+{\displaystyle \frac{c_p\left({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\right)}{L_{\mathrm{v}}}}\stackrel{\mathrm{̃}}{\mathrm{\Pi }}\sin \left(\mathit{\pi }\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)\right].\end{array}$$

The ratio of the first and second bracketed terms in Eq. (11), which gives the ratio of the IPV source from latent heating to the IPV source from mass sinks, scales like

$$\begin{array}{}\text{(12)}& A={\displaystyle \frac{\mathit{\pi }{L}_{\mathrm{v}}}{c_p\left({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\right)}}\approx \mathrm{200},\end{array}$$

using 40 K as an estimate of θ_t−θ_s that is appropriate for a latent heating profile that extends through the entire depth of the midlatitude troposphere. The expression for and approximate value of A – a simple measure of the relative magnitude of PV tendencies from latent heating and mass sinks – are the main results of this section. The value of A suggests that IPV anomalies produced by latent heating are likely to be about 200 times larger than IPV anomalies produced by mass sinks. Additionally, the fact that A depends only on thermodynamic constants and the potential temperature range within the region of latent heating suggests that its value is unlikely to vary strongly with climate. Both features are consistent with there being no detectable differences in cyclone statistics in our simulations without a mass sink and with the mass sink included. Our estimated value of A≈200 is also consistent with lower-tropospheric (Ertel) potential vorticity sources reported by LY04 for Hurricane Lili: in their simulation, lower-tropospheric PV sources from latent heating peak at around 25 PV units (PVU) h⁻¹, while lower-tropospheric sources from mass sinks peak at around 3–4 PVU d⁻¹. PV sources from mass sinks reported by LY04 are somewhat larger near the melting line, likely because of microphysical effects not considered in our theory, but remain about 70–80 times smaller than PV sources from latent heating.

Finally, we note that the value of A provides a plausible explanation for simulation results showing quantitatively detectable (albeit still minor) differences in midlatitude cyclone statistics between simulations without a mass sink and simulations with a 10-fold-exaggerated mass sink. With a 10-fold exaggeration of mass sinks relative to latent heating, our estimate of the ratio of IPV tendencies from latent heating and mass sinks is reduced to $A/\mathrm{10}\approx \mathrm{20}$. This value, although still significantly larger than 1, is small enough to suggest that the dynamical effects of mass sinks may not remain completely insignificant compared to the effects of latent heating: a 5 % ($\mathrm{1}/\mathrm{20}$) increase to the diabatic PV source is not a first-order effect but is not obviously negligible either.

4.2 Evolution of IPV anomalies

Our scale analysis in the previous section suggests that IPV tendencies from precipitation mass sinks are a factor of about 200 smaller than IPV tendencies produced by associated latent heating. However, IPV tendencies produced by latent heating and mass sinks have very different characteristic vertical structures: first-baroclinic latent heating produces an IPV source in the lower troposphere and an IPV sink in the upper troposphere, whereas the associated mass sink produces an IPV source everywhere. In terms of the effect on the induced winds, having different signs of IPV sources from latent heating may lead to weaker winds because the inversion from IPV to winds is spatially non-local and because vertical transport of IPV may bring diabatically generated IPV into a region where the diabatic generation is of a different sign. This effect is highlighted by LY04, who emphasize the role of spatial cancellation in limiting volume-integrated potential vorticity tendencies from latent heating. In this and the following section, we examine the ramifications of having regions with different-signed IPV sources from latent heating for the relative importance of latent heating versus mass sinks. We start by examining the temporal evolution of potential vorticity anomalies produced by the idealized latent heating and mass sink profiles introduced in the previous section.

Beginning with Eq. (11) and neglecting horizontal advection gives

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{P}}{\partial \stackrel{\mathrm{̃}}{t}}}=& -\sin \left(\mathit{\pi }\stackrel{\mathrm{̃}}{\mathit{\theta }}\right){\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{P}}{\partial \stackrel{\mathrm{̃}}{\mathit{\theta }}}}+\stackrel{\mathrm{̃}}{P}\mathit{\pi }\cos \left(\mathit{\pi }\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)\\ & +\stackrel{\mathrm{̃}}{P}{\displaystyle \frac{c_p\left({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\right)}{L_{\mathrm{v}}}}\stackrel{\mathrm{̃}}{\mathrm{\Pi }}\sin \left(\mathit{\pi }\stackrel{\mathrm{̃}}{\mathit{\theta }}\right).\end{array}\end{array}$$

We retain vertical advection (first term on the right-hand side) because vertical advection of anomalous potential vorticity plays an important role in setting potential vorticity profiles in mature cyclones (Schubert and Alworth, 1987; Büeler and Pfahl, 2017). The non-dimensional Exner function $\stackrel{\mathrm{̃}}{\mathrm{\Pi }}=(p/p_{\mathrm{s}}{)}^{R/c_p}$ should in principle vary in time as a circulation develops and p(θ) changes, but because $\stackrel{\mathrm{̃}}{\mathrm{\Pi }}$ only varies from about 1 to 0.5 as p varies from p_s to $p_{\mathrm{s}}/\mathrm{10}$, we approximate $\stackrel{\mathrm{̃}}{\mathrm{\Pi }}$ as fixed in time. We still need to specify its vertical structure, however, and for simplicity we assume that pressure varies linearly from p_s at θ=θ_s to (1−f_p)p_s at θ=θ_t, giving $p=p_{\mathrm{s}}(\mathrm{1}-f_p\stackrel{\mathrm{̃}}{\mathit{\theta }})$ and

$$\begin{array}{}\text{(14)}& \stackrel{\mathrm{̃}}{\mathrm{\Pi }}=(\mathrm{1}-f_p\stackrel{\mathrm{̃}}{\mathit{\theta }}{)}^{R/c_p}.\end{array}$$

We integrate Eq. (13) starting from $\stackrel{\mathrm{̃}}{P}=\mathrm{1}$ at $\stackrel{\mathrm{̃}}{t}=\mathrm{0}$ (a quiescent initial condition with zero relative vorticity) and ending at $\stackrel{\mathrm{̃}}{t}=\mathrm{2}$. We set ${\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}=\mathrm{40}$ K, corresponding to A≈200, and f_p=0.7, corresponding to latent heating that extends from 1000 to 300 hPa. The solution, which we plot as an anomaly relative to the initial condition and denote ${\stackrel{\mathrm{̃}}{P}}^{\prime }$, initially shows positive lower-tropospheric and negative upper-tropospheric IPV anomalies of equal depths (Fig. 6a). The opposite-signed anomalies are produced by latent heating (second term on the right-hand side of Eq. 13). However, the depth of the negative upper-tropospheric anomaly rapidly decreases, and negative IPV anomalies grow much less quickly than positive IPV anomalies. Vertical advection (first term on the right-hand side of Eq. 13), which transports positive IPV anomalies upward into regions with an IPV sink from latent heating, is responsible for decreasing the depth of the negative IPV anomaly. The factor of $\stackrel{\mathrm{̃}}{P}$ multiplying IPV sources and sinks (second and third right-hand side terms in Eq. 13) is responsible for limiting the growth of negative IPV anomalies because it causes $\stackrel{\mathrm{̃}}{P}$ to decay exponentially toward 0 in response to a sustained IPV sink, whereas positive IPV anomalies can grow without bound.

Figure 6 Evolution of IPV in a precipitating column and the contribution from the precipitation mass sink. Panel (a) shows the IPV anomaly ${\stackrel{\mathrm{̃}}{P}}^{\prime }$ that develops in response to latent heating, mass sinks, and vertical advection. Panel (b) shows the difference between ${\stackrel{\mathrm{̃}}{P}}^{\prime }$ and the IPV anomaly ${\stackrel{\mathrm{̃}}{P}}_{\text{noMS}}^{\prime }$ that develops with only latent heating and vertical advection included. This difference quantifies the influence of the mass sink on ${\stackrel{\mathrm{̃}}{P}}_{\text{both}}^{\prime }$. Panel (c) shows the magnitude of the relative difference between IPV anomalies that develop with and without the mass sink.

We quantify the influence of mass sinks on ${\stackrel{\mathrm{̃}}{P}}^{\prime }$ by integrating Eq. (13) a second time with the mass sink IPV source (third right-hand side term) removed. Subtracting the resulting IPV anomaly, which we denote as ${\stackrel{\mathrm{̃}}{P}}_{\text{noMS}}^{\prime }$, from ${\stackrel{\mathrm{̃}}{P}}^{\prime }$ shows that the mass sink contributes a positive IPV anomaly throughout the entire depth of the atmosphere (Fig. 6b). However, the contribution of the mass sink is small compared to the total IPV anomaly. Plotting $\left|\right({\stackrel{\mathrm{̃}}{P}}^{\prime }-{\stackrel{\mathrm{̃}}{P}}_{\text{noMS}}^{\prime })/{\stackrel{\mathrm{̃}}{P}}^{\prime }|$ (i.e., the magnitude of the relative difference between anomalies produced with and without mass sinks) shows that mass sink contributes only a fraction of about $\mathrm{0.5}\mathit{\%}\approx \mathrm{1}/A$ of the total positive IPV anomaly. In the lower troposphere ($\stackrel{\mathrm{̃}}{\mathit{\theta }}<\mathrm{0.5}$), a positive IPV tendency from upward transport and an IPV source from latent heating both contribute to the positive IPV anomaly and are both significantly larger than the IPV source from mass sinks. In the upper troposphere ($\stackrel{\mathrm{̃}}{\mathit{\theta }}>\mathrm{0.5}$), the combined positive tendency from upward transport plus negative tendency from latent heating eventually becomes positive and significantly larger than the IPV source from mass sinks. There is always a thin layer (initially near $\stackrel{\mathrm{̃}}{\mathit{\theta }}=\mathrm{0.5}$ and higher at later times) where the net IPV anomaly produced by compensating sinks from latent heating and sources from upward transport is small, and the relative difference between anomalies produced with and without mass sinks is fairly large (region between the 5 % contours in Fig. 6c). At long times ($\stackrel{\mathrm{̃}}{t}\gtrsim \mathrm{1}$), however, a positive IPV anomaly extends throughout most of the depth of the atmosphere, and the fractional contribution of mass sinks to the overall PV anomaly is small everywhere except in a thin layer near the top of the domain where the anomaly changes sign. So, despite the initial positive and negative IPV anomalies produced by latent heating being of opposite sign in the lower and upper troposphere, the long-term evolution of the IPV field is toward a vertically extensive positive IPV anomaly produced almost entirely by latent heating plus vertical advection, with only a very small fractional contribution from mass sinks.

4.3 Spatial cancellation during inversion

At short times ($\stackrel{\mathrm{̃}}{t}\ll \mathrm{1}$), however, IPV anomalies produced by latent heating do exhibit positive and negative regions of similar magnitude. Because of the non-local relationship between IPV anomalies and induced circulations, this dipole structure has the potential to limit the strength of circulations induced by latent heating. In this section, we use linearized IPV inversions to examine whether the relative strength of circulations induced by latent heating and mass sinks scales with the IPV tendency ratio A or whether spatial cancellation of IPV anomalies during inversion significantly weakens circulations induced by latent heating relative to those induced by mass sinks.

Substituting $\stackrel{\mathrm{̃}}{P}={\stackrel{\mathrm{̃}}{P}}^{\prime }+\mathrm{1}$ into Eq. (13), assuming that ${\stackrel{\mathrm{̃}}{P}}^{\prime }\ll \mathrm{1}$ (valid at short times), and dropping right-hand-side terms that involve ${\stackrel{\mathrm{̃}}{P}}^{\prime }$ gives a short-time approximation to Eq. (13):

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{P}}^{\prime }}{\partial \stackrel{\mathrm{̃}}{t}}}=\mathit{\pi }\cos \left(\mathit{\pi }\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)+{\displaystyle \frac{c_p\left({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\right)}{L_{\mathrm{v}}}}\stackrel{\mathrm{̃}}{\mathrm{\Pi }}\sin \left(\mathit{\pi }\stackrel{\mathrm{̃}}{\mathit{\theta }}\right).\end{array}$$

Assuming that the strength of latent heating and mass sinks decays like a Gaussian with increasing radius r provides expressions for short-time IPV anomalies produced by latent heating,

$$\begin{array}{}\text{(16)}& {\stackrel{\mathrm{̃}}{P}}_{\text{heat}}^{\prime }=\cos \left(\mathit{\pi }\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)\exp (-{\stackrel{\mathrm{̃}}{r}}^{\mathrm{2}}),\end{array}$$

and mass sinks,

$$\begin{array}{}\text{(17)}& {\stackrel{\mathrm{̃}}{P}}_{\text{mass}}^{\prime }={\displaystyle \frac{\mathrm{1}}{A}}\stackrel{\mathrm{̃}}{\mathrm{\Pi }}\sin \left(\mathit{\pi }\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)\exp (-{\stackrel{\mathrm{̃}}{r}}^{\mathrm{2}}),\end{array}$$

where $\stackrel{\mathrm{̃}}{r}=r/L$ and L is the horizontal length scale (Fig. 7a and b). Because we compute relative circulation strengths using a linearized inversion relation, our results depend only on the relative magnitudes of IPV anomalies, and we use solutions to Eq. (15) at $\stackrel{\mathrm{̃}}{t}=\mathrm{1}/\mathit{\pi }$ for simplicity even though the short-time approximation (${\stackrel{\mathrm{̃}}{P}}^{\prime }\ll \mathrm{1}$) breaks down before then.

Figure 7 Short-time IPV anomalies produced by latent heating (a) and mass sinks (b), and azimuthal velocity fields induced by latent heating (c) and mass sinks (d) from a linearized inversion with B=1.

Small-amplitude axisymmetric IPV anomalies and azimuthal velocities are approximately related by a linear elliptic equation (Eq. 32 of Hoskins et al., 1985):

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \left(rv\right)}{\partial r}}\right)+{\displaystyle \frac{f^{\mathrm{2}}}{g{\mathit{\sigma }}_{\text{ref}}}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left({\left({\displaystyle \frac{\partial \mathrm{\Pi }}{\partial p}}\right)}_{\text{ref}}^{-\mathrm{1}}{\displaystyle \frac{\partial v}{\partial \mathit{\theta }}}\right)={\mathit{\sigma }}_{\text{ref}}{\displaystyle \frac{\partial P}{\partial r}}.\end{array}$$

Here, v is the azimuthal velocity, σ_ref is a reference-state isentropic density profile, and $(\partial \mathrm{\Pi }/\partial p{)}_{\text{ref}}$ is the pressure derivative of the Exner function in the reference state. Our earlier assumption that pressure varies linearly with potential temperature means that σ_ref is a constant ${\mathit{\sigma }}_{\mathrm{0}}=g^{-\mathrm{1}}{f}_p{p}_{\mathrm{s}}/({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}})$. Non-dimensionalizing by substituting

$$\begin{array}{ll}{\displaystyle }r& {\displaystyle }\to L\stackrel{\mathrm{̃}}{r}\\ {\displaystyle }v& {\displaystyle }\to fL\stackrel{\mathrm{̃}}{v}\\ {\displaystyle }p& {\displaystyle }\to p_{\mathrm{s}}\stackrel{\mathrm{̃}}{p},\end{array}$$

in addition to earlier non-dimensionalizations, gives

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\partial }{\partial \stackrel{\mathrm{̃}}{r}}}\left({\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{r}}}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{̃}}{r}\stackrel{\mathrm{̃}}{v}\right)}{\partial \stackrel{\mathrm{̃}}{r}}}\right)+{\displaystyle \frac{\mathrm{1}}{B}}{\displaystyle \frac{\partial }{\partial \stackrel{\mathrm{̃}}{\mathit{\theta }}}}\left({\stackrel{\mathrm{̃}}{\mathrm{\Pi }}}_{p,\text{ref}}^{-\mathrm{1}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{v}}{\partial \stackrel{\mathrm{̃}}{\mathit{\theta }}}}\right)={\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{P}}{\partial \stackrel{\mathrm{̃}}{r}}}\end{array}$$

as a non-dimensional inversion relation, with

$$\begin{array}{}\text{(20)}& {\stackrel{\mathrm{̃}}{\mathrm{\Pi }}}_{p,\text{ref}}=(\mathrm{1}-f_p\stackrel{\mathrm{̃}}{\mathit{\theta }}{)}^{R/c_p-\mathrm{1}}.\end{array}$$

B is a non-dimensional Burger number and can be written either as

$$\begin{array}{}\text{(21)}& B={\displaystyle \frac{L_{\mathrm{D}}^{\mathrm{2}}}{L^{\mathrm{2}}}},\end{array}$$

where

$$\begin{array}{}\text{(22)}& L_{\mathrm{D}}={\displaystyle \frac{{\left(g{\mathit{\sigma }}_{\mathrm{0}}R/p_{\mathrm{s}}\right)}^{\mathrm{1}/\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\right)}{f}}\end{array}$$

is the Rossby deformation radius, or as

$$\begin{array}{}\text{(23)}& B={\displaystyle \frac{{\left({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\right)}^{\mathrm{2}}}{\mathrm{\Delta }{\mathit{\theta }}^{\mathrm{2}}}},\end{array}$$

where

$$\begin{array}{}\text{(24)}& \mathrm{\Delta }\mathit{\theta }={\displaystyle \frac{fL}{{\left(g{\mathit{\sigma }}_{\mathrm{0}}R/p_{\mathrm{s}}\right)}^{\mathrm{1}/\mathrm{2}}}}\end{array}$$

measures the vertical penetration of the inversion operator (Hoskins et al., 1985). Thus, we expect spatial cancellation during inversion to reduce the strength of latent-heating-induced circulations the most when Δθ is large compared to θ_t−θ_s or, equivalently, when disturbances are large compared to a deformation radius.

To examine how spatial cancellation during inversion affects the relative strength of circulations induced by latent heating versus mass sinks, we use Eq. (19) to invert ${\stackrel{\mathrm{̃}}{P}}_{\text{heat}}^{\prime }$ (Eq. 16) and ${\stackrel{\mathrm{̃}}{P}}_{\text{mass}}^{\prime }$ (Eq. 17) for azimuthal velocity fields ${\stackrel{\mathrm{̃}}{v}}_{\text{heat}}$ and ${\stackrel{\mathrm{̃}}{v}}_{\text{mass}}$, respectively. We perform inversions numerically on a domain with 1024 grid points between $\stackrel{\mathrm{̃}}{r}=\mathrm{0}$ and $\stackrel{\mathrm{̃}}{r}=\mathrm{5}$ and 256 grid points between $\stackrel{\mathrm{̃}}{\mathit{\theta }}=\mathrm{0}$ and $\stackrel{\mathrm{̃}}{\mathit{\theta }}=\mathrm{1}$, subject to boundary conditions $\stackrel{\mathrm{̃}}{v}=\mathrm{0}$ at the outer boundary and ${\partial }_{\stackrel{\mathrm{̃}}{\mathit{\theta }}}\stackrel{\mathrm{̃}}{v}=\mathrm{0}$ (no pressure perturbations) at the top and bottom boundaries. We assume that $\stackrel{\mathrm{̃}}{v}/\stackrel{\mathrm{̃}}{r}$ (or, equivalently, the relative vorticity ζ_θ) tends toward a constant as r approaches 0. With B=1, this produces a low-level cyclonic circulation and upper-level anticyclonic circulation for ${\stackrel{\mathrm{̃}}{v}}_{\text{heat}}$ and a vertically extensive but much weaker cyclonic circulation for ${\stackrel{\mathrm{̃}}{v}}_{\text{mass}}$ (Fig. 7c and d).

We quantify the relative strength of circulations induced by latent heating and mass sinks by calculating latent-heating-to-mass-sink ratios of the maximum IPV anomaly, maximum azimuthal wind speed, and maximum near-surface azimuthal wind speed for inversions across a broad range of Burger numbers (Fig. 8). The ratio of maximum IPV anomalies ($max\left({\stackrel{\mathrm{̃}}{P}}_{\text{heat}}^{\prime }\right)/max\left({\stackrel{\mathrm{̃}}{P}}_{\text{mass}}^{\prime }\right)$) is independent of B and close to the IPV tendency ratio A. The ratio of maximum azimuthal wind speeds ($max\left|{\stackrel{\mathrm{̃}}{v}}_{\text{heat}}\right|/max\left|{\stackrel{\mathrm{̃}}{v}}_{\text{mass}}\right|$) approaches zero for small B (horizontal length scales much larger than a deformation radius) because of vertical cancellation between opposite-signed regions of IPV anomalies from latent heating. For B∼1, however, vertical cancellation only moderately reduces the strength of circulations induced by latent heating, and the maximum wind speed induced by latent heating is larger than the maximum wind speed induced by mass sinks by a factor of about $A/\mathrm{2}$. For larger B, the ratio of maximum wind speeds asymptotes to a value equal to the ratio of maximum IPV anomalies. Ratios of maximum near-surface wind speeds behave similarly at small and intermediate Burger numbers but diverge at large Burger numbers because IPV anomalies from mass sinks are small near the surface.

Figure 8 Ratios of maximum IPV anomalies (gray), azimuthal wind speeds (solid black), and near-surface azimuthal wind speeds (dashed black) produced by latent heating and the mass sink in inversions with a varying Burger number $B=L_{\mathrm{D}}^{\mathrm{2}}/L^{\mathrm{2}}$.

These results, while likely to be sensitive to some degree to choices about the spatial structure of IPV anomalies, suggest that spatial cancellation during inversion between opposite-signed IPV anomalies is unlikely to significantly reduce the strength of circulations induced by latent heating compared to those induced by mass sinks. Typical length scales for precipitating regions of midlatitude disturbances are of the order of a deformation radius (e.g., extratropical cyclones) or smaller (e.g., across fronts). In these regimes (B≳1), our inversions suggest that circulations induced by mass sinks are likely to be at least ∼100 times weaker than circulations induced by latent heating. This ratio is somewhat smaller than the IPV tendency ratio A but only by a factor of about 2.

Like our initial scale analysis of IPV tendencies, results from idealized calculations of the evolution of IPV in a precipitating column and related inversions are consistent with simulations showing no significant differences between runs without a mass sink and with the mass sink included. The idealized evolution and inversions also suggest that a 10-fold exaggeration of mass sinks should allow them to induce circulations that are within about an order of magnitude of the strength of circulations induced by latent heating, which is at least plausibly consistent with simulation results showing a detectable effect from 10-fold-exaggerated mass sinks. However, our theoretical analysis – which suggests that mass sinks should play a fairly unimportant role in both developing and mature disturbances – is difficult to square with the results of LY04, who demonstrated a noticeable increase in TC intensity in a simulation of Hurricane Lili with an unexaggerated precipitation mass sink included. In the following section, we attempt to reconcile our results with LY04 by testing the effect of mass sources and sinks in TC world simulations.

5 Impact of the mass sink in TC world simulations

We generate TC worlds in the GFDL idealized spectral model by running simulations at a resolution of T170 with the sea surface temperature (SST) fixed to 300 K and the Coriolis parameter set to a value corresponding to 30^∘ latitude everywhere in the model domain. The choice of SST is not critical for TC genesis – simulations with globally uniform SST can generate tropical cyclones at SSTs far above and below those found in TC genesis regions on modern Earth (e.g., Merlis et al., 2016) – but is likely to influence precipitation rates. As discussed below, our choice of SST generates precipitation rates comparable to LY04's Eta model simulations. We replace the gray radiative transfer scheme with fixed radiative cooling of 1.5 K d⁻¹ below the tropopause, defined as the lowest level with temperature less than 200 K, and temperature relaxation to 200 K on a 40 d timescale above the tropopause, and we run simulations with large-scale condensation but no convection scheme. Both changes are necessary to obtain TC worlds: simulations with gray radiation or the simplified Betts–Miller convection scheme (Frierson, 2007) struggle to maintain axisymmetric TC-like vortices. Additionally, we use ∇⁴ hyperdiffusion instead of ∇⁸ hyperdiffusion, increase the hyperdiffusion coefficients for vorticity and temperature by a factor of 5 (relative to our conventional T170 simulations), and increase the hyperdiffusion coefficient for divergence by a factor of 10. Changes to damping parameters are motivated by suggestions from Zhao et al. (2012) that stronger divergence damping may increase the number of TCs produced by global models. We initialize simulations with a quiescent isothermal atmosphere, run simulations at each planetary radius for 600 d without mass sources and sinks (γ=0, no MS), and repeat the simulations with mass sources and sinks (γ=1, MS). Statistics are calculated using instantaneous snapshots saved every 6 h.

Figure 9 Snapshots of near-surface wind speed in no-MS TC world simulations with the planetary radius equal to Earth's radius (a), half of Earth's radius (b), and one-quarter of Earth's radius (c), as well as the total accumulated precipitation over the 3 h period following each snapshot (d–f). Values in the upper-left corner of panels (a–c) show a measure of the maximum surface wind speed in a typical TC, defined as the maximum wind speed for the median TC (ranking all TCs at a given time by their maximum wind speed) averaged in time over the last 400 d of each simulation.

Simulations with this model configuration produce a population of small axisymmetric TC-like vortices (Fig. 9a). At a resolution of T170, which has an effective grid spacing of about 80 km at the Equator and 60 km at midlatitudes, the vortices are only marginally resolved. Computational expense prevents us from easily running TC world simulations at higher resolution, so we instead test robustness to resolution by running additional simulations with the planetary radius reduced to one-half and one-quarter of Earth's radius. Smaller planetary radii produce simulations with better-resolved vortices but lower vortex counts (Fig. 9b and c). Quarter-Earth-radius simulations have an effective resolution of 20 km at the Equator, comparable to the 15 km resolution of Eta model simulations run by LY04. In-storm precipitation rates are likewise comparable in our TC world simulations and LY04's Eta model simulations: 3 h precipitation accumulations of 5–20 mm are widespread around TC cores, and accumulations of 20–60 mm occur over smaller regions (compare Fig. 9d–f and LY04's Fig. 10). Although we compare results in this section primarily to TC simulations analyzed by LY04, we note that the simulations of midlatitude frontal simulations described by Qiu et al. (1993) also use a horizontal resolution (80 km) within the range of our TC world simulations and produce peak precipitation rates (10 mm h⁻¹) comparable to our TC world simulations. Finally, we emphasize that the 20 to 40 km effective grid spacing of quarter- and half-Earth-radius simulations is comparable to resolutions typically used for global modeling of TCs (e.g., Zhao et al., 2009; Merlis et al., 2016; Chavas and Reed, 2019). Simulations at even higher resolution would be useful for assessing robustness, but computational expense prevents us from further decreasing the planetary radius (which requires a shorter time step), and TC size is moreover already comparable to the planetary surface area in the quarter-Earth-resolution simulations.

We analyze TC world simulations using a TC identification algorithm following Cronin and Chavas (2019). For each 6-hourly snapshot, we mask regions with a surface pressure anomaly, measured relative to the zonal-mean surface pressure, less than −20 hPa. We tile the masked pressure field three times in longitude to permit identification of cyclones that cross longitudinal boundaries of the model output and label contiguous masked regions using an eight-connected-component labeling algorithm. For each contiguous masked region, we define the latitude and longitude coordinates of the cyclone center as the surface-pressure-anomaly-weighted centroid of the 13 grid points with the lowest surface pressure anomaly. Unlike Cronin and Chavas (2019), we retain masked regions with less than 13 grid points, in which case we include all grid points in the centroid calculation. We retain cyclones with centroids between 0 and 360^∘ longitude and −80 to 80^∘ latitude, the latter to avoid double-counting storms near the poles. We then define a central pressure anomaly for each cyclone by interpolating the surface pressure anomaly field to the cyclone centroid and a maximum surface wind speed for each cyclone by finding the maximum surface wind speed within the region with a surface pressure anomaly less than −20 hPa. The maximum wind speed within the median TC (ranking all TCs at a given time by their maximum wind speed) is highest in quarter-Earth-radius simulations, indicating that wind speeds in typical TCs increase at higher resolution (Fig. 9a–c). By contrast, global-maximum surface wind speeds are highest in Earth-radius simulations simply because they contain more storms sampled from an intensity distribution with finite width.

Figure 10 shows time series of the central pressure anomaly, averaged across storms, and cyclone count for each TC world simulation. The first identified cyclones appear after several tens of days, and the cyclone count increases rapidly before dropping and equilibrating to a statistically steady value after 200 d. Cyclone count scales approximately linearly with planetary surface area, with about 50 cyclones at steady state in Earth-radius simulations, about 15 at steady state in half-Earth-radius simulations, and 4–5 at steady state in quarter-Earth-radius simulations. The inclusion or omission of mass sources and sinks has no clear effect on TC count. Time-average counts taken over four 100 d blocks between days 200 and 600 are similar, with substantial overlap in ranges of values obtained for different 100 d blocks. Cyclone central pressure anomalies increase in magnitude over the first 100 d and also reach statistical steady states after 200 d. Central pressure anomalies averaged across storms increase slightly with increasing resolution, from about −45 hPa in Earth-radius simulations to about −50 hPa in quarter-Earth-radius simulations. Like cyclone counts, time-average central pressure anomalies are similar with and without mass sources and sinks, with significant overlap between the ranges of time averages taken over four 100 d blocks between days 200 and 600. Additional quarter-Earth-radius experiments without changes to hyperdiffusion parameters produce TCs with slightly weaker central pressure anomalies and likewise provide no evidence that precipitation mass sinks produce stronger TCs, though they suggest that weaker damping may allow mass sinks to increase TC count at the expense of TC intensity (Appendix E).

Figure 10 Time series of average TC central pressure anomalies (a, c, e) and TC count (b, d, f) in no-MS (black) and MS (blue) TC world simulations with the planetary radius equal to Earth's radius (a, b), half of Earth's radius (c, d), and one-quarter of Earth's radius (e, f). The vertical bars at the end of the time series indicate the range of time averages taken over four 100 d blocks between days 200 and 600.

So far, results from our TC world simulations (which show no detectable increase in TC intensity from including mass sources and sinks) are difficult to reconcile with LY04 (who found a clear intensity increase in simulations with mass sources and sinks). The disagreement cannot clearly be attributed to differences in resolution, as we detect no effect from including mass sources and sinks in quarter-Earth-radius simulations with resolutions and precipitation rates comparable to LY04's Eta model simulations. The disagreement is also unlikely to be the result of differences between rapidly deepening versus mature cyclones: LY04's Eta model simulations produced very little deepening over the simulation period, and our TC world simulations show no clear difference between no-MS and MS simulations during the initial deepening period (Fig. 10).

One potentially important difference between our TC world simulations and the simulations presented by LY04 is that our simulations consider properties of a population of cyclones in equilibrium with their environment – and so with intensities potentially constrained by global budgets – while LY04 consider the transient evolution of a single TC in a short-term forecast. To test whether effects from mass sources and sinks are more readily apparent in transient simulations initialized from a common starting state, we restart two simulations from the end of the Earth-radius TC world simulation without mass sources and sinks. In one simulation, we continue to omit mass sources and sinks; in the other simulation, we include them. We run these two simulations for 3 d and save instantaneous snapshots once every hour.

Once these two simulations are complete, we use our cyclone identification algorithm to calculate central pressure anomalies for each cyclone and match corresponding TCs in the two simulations using the following algorithm. For each instantaneous snapshot, we consider a “source” set of cyclones from the simulation with fewer identified cyclones and a “destination” set of cyclones from the other simulation. We match each cyclone in the source set with the nearest cyclone (measured by the great circle distance between centroids) in the destination set. This results in between 46 and 53 matched pairs for each snapshot. The matching is unique (i.e., there are not cases where multiple cyclones in the source set are matched to a single cyclone in the destination set) for the first 56 h and for 9 of the last 16 h. In the remaining 7 of the last 16 h, the number of unique cyclones matched from the destination set is one less than the number of cyclones in the source set, indicating that two cyclones in the source set were matched to the same cyclone in the destination set. Finally, we calculate the pairwise difference between central pressure anomalies for each matched pair.

Figure 11 Transient evolution of pairwise differences in cyclone central pressure over the first 3 d of TC world restart experiments. The median and q=0.1 and =0.9 quantiles are calculated over the entire population of cyclones, and corresponding cyclones in experiments with and without the mass sink are matched as described in the text.

Figure 11 shows time series of the median, 10th percentile, and 90th percentile pairwise difference between central pressure anomalies for each snapshot. The median central pressure difference varies little over the 3 d of simulation, with a slight decrease (less than 1 hPa) during the last 24 h. However, the spread in 10th and 90th percentile central pressure differences increases substantially over the 3 d. The 10th percentile central pressure difference – which captures a pair of cyclones with a larger central pressure anomaly in the simulation with mass sources and sinks – decreases to about −8 hPa by day 3. The 90th percentile central pressure difference – which captures a pair of cyclones with a larger central pressure anomaly in the simulation without mass sources and sinks – departs from 0 slightly more slowly and reaches about 5 hPa by day 3.

This suggests one way to reconcile our results with LY04: including mass sources and sinks can have little effect on the intensity of an equilibrium population of TCs and no obvious effect on the median intensity change in transient simulations of a large population of TCs, but still significantly alter the evolution of individual cyclones over timescales of 1–3 d. Our simulations do not show systematic TC deepening in simulations with mass sources and sinks included but can produce substantial transient deepening for some individual cyclones, similar to results reported by LY04 for Hurricane Lili. To be clear, we have not ruled out the possibility that model differences may bias our results relative to LY04. TCs in our simulations are relatively coarsely resolved, and obtaining TC-like vortices required running simulations without a convection scheme, both of which may affect our representation of TC dynamics. Nevertheless, our analysis of short-term TC world simulations provides a plausible bridge between our results for statistical equilibrium and the single-storm results of LY04.

6 Conclusions

In this paper, we presented idealized simulations and simple theoretical arguments that, together, suggest that water vapor mass sources and sinks are unlikely to have a significant impact on large-scale midlatitude circulations. In simulations with an idealized atmospheric general circulation model, statistics capturing midlatitude variability and extremes changed little when the model's dynamical core was modified to include mass sources and sinks in the mass conservation equation (Sects. 2 and 3). This result held even in very warm climates with global-mean surface temperatures above 315 K and atmospheric water vapor content much larger than that of modern Earth.

Our simulation results were consistent with analysis of an isentropic potential vorticity equation, which showed that PV sources and streamfunction anomalies associated with latent heating are typically larger than those associated with mass sinks by a non-dimensional factor of $\mathit{\pi }{L}_{\mathrm{v}}/\left(c_p\right({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\left)\right)$, where L_v is the latent heat of vaporization of water vapor, c_p is the heat capacity of dry air, and θ_t−θ_s is the potential temperature range over which condensation occurs. For troposphere-deep ascent in Earth's atmosphere (${\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\approx \mathrm{40}$ K), this non-dimensional number is about 200, suggesting that latent heating should produce PV sources and sinks about 2 orders of magnitude larger than those produced by mass sources and sinks. Warmer climates will produce storms with much heavier precipitation rates, but the ratio of PV sources for latent heating versus mass sinks will remain similar. Re-purposing this result for other planets with other condensible species is straightforward: for troposphere-deep ascent on Titan (${\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\approx \mathrm{20}$ K, e.g., Mitchell and Lora, 2016), where methane is the main condensible species ($L_{\mathrm{v}}\approx \mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}$ J kg⁻¹) and the atmosphere is mostly nitrogen (c_p≈10³ $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), $\mathit{\pi }{L}_{\mathrm{v}}/\left(c_p\right({\mathit{\theta }}_{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{s}}\left)\right)\approx \mathrm{80}$, suggesting that PV sources and sinks from troposphere-deep latent heating should be a factor of about 80 larger than those from corresponding mass sinks.

Our simulation results stand in contrast to previous work that found that the precipitation mass sink could enhance precipitation rates in mesoscale weather systems by about 10 % (Qiu et al., 1993) and increase peak surface wind speeds in hurricanes by 5 to 15 kn (Lackmann and Yablonsky, 2004). Based on TC world simulations, we suggest that one way to reconcile our results with previous work is by recognizing that including mass sources and sinks can substantially alter the evolution of individual weather systems on timescales of days without having a systematic effect on populations of many individual storms (Sect. 5).

Our results reveal no urgent need to revisit existing literature on links between climate and statistical properties of midlatitude storm tracks. Our results directly confirm that a reduction in midlatitude EKE with warming at statistical equilibrium on an aquaplanet, documented in previous work (O'Gorman and Schneider, 2008a), remains virtually unchanged in simulations that include mass sources and sinks. Our work did not examine Lagrangian statistics that measure, e.g., the size, intensity, or frequency of midlatitude cyclones themselves (Pfahl et al., 2015), but our simulation results and analysis of PV theory suggest that those are also unlikely to depend significantly on the inclusion of mass sources and sinks. Our work also did not explore climates in which water vapor becomes extremely non-dilute and the character of large-scale flows becomes very different (Pierrehumbert and Ding, 2016). Including mass sources and sinks for condensible species certainly broadens the range of applications for large-scale atmospheric models, and they are included in many models that target simulations of planetary atmospheres (Pierrehumbert and Ding, 2016; Forget et al., 2017; Neary and Daerden, 2018; Chow et al., 2019). The Martian general circulation in particular includes a planetary-scale meridional component that is driven by seasonal condensation of carbon dioxide at the winter pole and is clearly influenced by changes in surface pressure produced by carbon dioxide mass sinks (Pollack et al., 1981, 1990; Haberle et al., 1993; Chow et al., 2019). However, omitting mass sources and sinks seems to be a reasonable simplification for simulations of midlatitude eddies in Earth's atmosphere, even in climates substantially warmer than the modern.

At the same time, some of our results – specifically, our short-term TC world simulations – do indicate that mass sources and sinks can have appreciable effects on the transient evolution of strongly precipitating weather systems, as was suggested by LY04. This raises the possibility that short-term forecasts may benefit from correctly including mass sources and sinks. In our view, however, it remains unclear whether mass sources and sinks are “missing physics” with the potential to provide systematic improvements to forecast skill or whether their inclusion simply provides a small perturbation that eventually alters the trajectory of a chaotic dynamical system. This issue could be addressed in future work using simulations similar to our TC world restart experiments, or additional simulations of TC case studies, with a focus on comparing how trajectories diverge in pairs of simulations with and without mass sources and sinks and in pairs of simulations distinguished by small random perturbations to initial conditions.

Overall, our study offers evidence, from both simulations and theory, suggesting that mass sources from evaporation and mass sinks from precipitation are unlikely to have a significant effect on large-scale midlatitude variability. Our results do not necessarily preclude the possibility that mass sources and sinks may have a more important role on the transient evolution of strongly precipitating mesoscale weather systems, including midlatitude fronts (Qiu et al., 1993) and tropical cyclones (Lackmann and Yablonsky, 2004), as illustrated for the transient evolution of individual tropical cyclones in Sect. 5. Considering the role of precipitation mass sinks in an idealized frontogenesis problem could help to clarify their role in the dynamics of mesoscale midlatitude weather systems. Future work could also examine the effects of mass sources and sinks on the transient evolution and equilibrium statistics of TCs in idealized cloud-system-resolving simulations, which produce TC worlds when run on an f-plane (e.g., Merlis and Held, 2019). The role of mass sinks in the formation, organization, and vigor of individual convective cells could be explored through simulations with cloud-resolving models based on fully compressible non-hydrostatic equations that carefully represent the effects of precipitation fallout on mass, energy, entropy, and momentum budgets (e.g., Ooyama, 2001; Bott, 2008). Finally, models that include a more sophisticated treatment of condensed water (including condensate loading and prognostic hydrometeor fallout) could explore whether processes neglected in our pseudoadiabatic model enhance the dynamical effects of precipitation mass sinks. All are potentially interesting directions for future research.

Appendix A: Model changes

Modifications to the idealized GCM couple sources and sinks of moisture calculated by physics parameterizations to pressure tendencies calculated by the dynamical core while retaining the ability to run simulations into equilibrium across a wide range of climates. We emphasize that the changes are not intended to produce a model fully valid for simulations with strongly non-dilute water vapor, and the model retains approximations – like the use of a humidity-independent heat capacity – that become increasingly inaccurate as the moisture content of the atmosphere increases. Rather, the goal of the changes documented here is to add a representation of the precipitation mass sink, or pressure tendencies resulting from sources and sinks of water vapor, while making relatively few other changes to the model equations.

Incorporating the precipitation mass sink is done in three steps. First, the pressure tendency calculation in the dynamical core is modified to include a term for a water vapor mass source (Sect. A1). Second, mass sources are calculated from specific humidity tendencies produced by the convection, large-scale condensation, and surface flux schemes and passed to the modified pressure tendency calculation (Sect. A2). Finally, a mass fixer is added to the model dynamical core to correct numerical errors in the conservation of dry mass introduced by the inclusion of sources and sinks of moist mass (Sect. A3).

In addition to changes made to incorporate the precipitation mass sink, we make some minor changes to lapse rate calculations in the convection scheme (Sect. A4) to improve consistency with the treatment of thermodynamics and moisture in other parts of the GCM.

A1 Pressure tendency calculations

Pressure tendency calculations are derived from the mass continuity equation. For a parcel with dry-mass density ρ_d; total density ρ; water vapor mixing ratio r_v; and spatial dimensions δx, δy, and δz, mass conservation requires

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\mathit{\rho }\mathit{\delta }x\mathit{\delta }y\mathit{\delta }z={\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}(\mathrm{1}+r_{\mathrm{v}}){\mathit{\rho }}_{\mathrm{d}}\mathit{\delta }x\mathit{\delta }y\mathit{\delta }z={\mathit{\rho }}_{\mathrm{d}}\mathit{\delta }x\mathit{\delta }y\mathit{\delta }z{\displaystyle \frac{\mathrm{D}{r}_{\mathrm{v}}}{\mathrm{D}t}},\end{array}$$

where the second equality follows from the fact that the parcel conserves dry mass. From hydrostatic balance,

$$\begin{array}{}\text{(A2)}& \mathit{\delta }x\mathit{\delta }y\mathit{\delta }z=-{\displaystyle \frac{\mathit{\delta }x\mathit{\delta }y\mathit{\delta }p}{\mathit{\rho }g}},\end{array}$$

where p is pressure and g is the acceleration from gravity, giving

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\delta }x\mathit{\delta }y\mathit{\delta }p}}{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\mathit{\delta }x\mathit{\delta }y\mathit{\delta }p={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{d}}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{D}{r}_{\mathrm{v}}}{\mathrm{D}t}}.\end{array}$$

Replacing $\mathrm{D}\mathit{\delta }x/\mathrm{D}t$ with δu and $\mathit{\delta }u/\mathit{\delta }x$ with ∂_xu and making similar substitutions for y and p gives

$$\begin{array}{}\text{(A4)}& \mathrm{\nabla }\cdot \mathbf{v}+{\displaystyle \frac{\partial \mathit{\omega }}{\partial p}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{d}}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{D}{r}_{\mathrm{v}}}{\mathrm{D}t}},\end{array}$$

where $\mathrm{\nabla }=(\partial /\partial x,\partial /\partial y)$ taken at constant p, $\mathbf{v}=(u,v)$ is the horizontal velocity, and ω is the vertical velocity in pressure coordinates. Replacing the mixing ratio with specific humidity q_v gives

$$\begin{array}{}\text{(A5)}& \mathrm{\nabla }\cdot \mathbf{v}+{\displaystyle \frac{\partial \mathit{\omega }}{\partial p}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-q_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{D}{q}_{\mathrm{v}}}{\mathrm{D}t}}.\end{array}$$

Within the model, we replace ${\partial }_p\mathit{\omega }=-\mathrm{\nabla }\cdot \mathbf{v}$ in pressure tendency calculations with

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\partial \mathit{\omega }}{\partial p}}=-\mathrm{\nabla }\cdot \mathbf{v}+\mathit{\gamma }{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-q_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{D}{q}_{\mathrm{v}}}{\mathrm{D}t}}.\end{array}$$

Setting an appropriate value of γ allows simulations without mass sources and sinks (γ=0), with mass sources and sinks (γ=1), and with exaggerated mass sources and sinks (γ>1). In pressure coordinates, integrating downward from ω=0 at p=0 gives

$$\begin{array}{}\text{(A7)}& \mathit{\omega }\left(p\right)=-\mathrm{\nabla }\cdot \underset{\mathrm{0}}{\overset{p}{\int }}\mathbf{v}\mathrm{d}{p}^{\prime }+\underset{\mathrm{0}}{\overset{p}{\int }}\mathit{\gamma }\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}-q_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{D}{q}_{\mathrm{v}}}{\mathrm{D}t}}\right)\mathrm{d}{p}^{\prime }.\end{array}$$

The second term is the contribution to the pressure tendency ω(p) from the source of water vapor described by $\mathrm{D}{q}_{\mathrm{v}}/\mathrm{D}t$. Within the sigma-coordinate dynamical core, the modified expression for ∂_pω (Eq. A6) is used by existing (un-modified) code to calculate surface pressure tendencies and sigma-coordinate vertical velocities that incorporate contributions from mass sources and sinks.

A2 Mass sources and sinks from moist physics

We include mass sources and sinks from the large-scale condensation scheme, the convection scheme, and surface evaporation.

The large-scale condensation scheme produces specific humidity tendencies from condensation and evaporation, and we add resulting mass sources and sinks to the pressure tendency calculation in the dynamical core following Eq. (A6).

The convection scheme also produces specific humidity tendencies but without distinguishing between tendencies from condensation and evaporation (which add and remove mass from parcels) and tendencies from sub-grid-scale transport (which conserves mass following parcels). We assume that specific humidity tendencies are entirely from parcel mass sources and sinks in columns where the convection scheme is running in “deep-convection” mode (which produces some surface precipitation) and entirely from sub-grid transport in columns where the convection scheme is running in “shallow-convection” mode (which produces no surface precipitation). Thus, we add mass sources and sinks to the pressure tendency calculation following Eq. (A6) for columns and times with parameterized deep convection but not parameterized shallow convection.

The surface flux scheme provides the mass flux of water vapor from surface evaporation directly, and we add a mass source to the lowest model level based on the mass flux from surface evaporation.

A3 Mass fixer

The implementation of water vapor mass sources and sinks produces a model that only conserves dry mass up to first order in the size of specific humidity increments, and the resulting numerical errors result in a gradual downward drift of dry atmospheric mass. This downward drift of dry mass occurs even if moist mass sources and sinks from evaporation and precipitation are exactly in balance (see Eq. A12) and prevents stable simulations if left uncorrected. We compensate for unphysical reductions in dry mass produced by numerical errors by using a mass fixer that ensures that the model conserves a fixed inventory of dry mass while leaving an additional component of atmospheric mass free to adjust in response to moist mass sources and sinks.

We can show that the numerical errors enter at second order and tend to reduce dry atmospheric mass by considering a parcel with mass m and specific humidity q_v before and after an update that alters specific humidity by an amount dq and mass by an amount dm. Before the update, the parcel mass is related to specific humidity and dry mass m_d by

$$\begin{array}{}\text{(A8)}& m={\displaystyle \frac{m_{\mathrm{d}}}{\mathrm{1}-q_{\mathrm{v}}}}.\end{array}$$

Because mass changes are computed assuming that parcels conserve dry mass (see Sect. A1), the parcel mass after the update is given by

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{m}^{\prime }& =m+\mathrm{d}m=m+\mathrm{d}\left({\displaystyle \frac{m_{\mathrm{d}}}{\mathrm{1}-q_{\mathrm{v}}}}\right)\\ & =m+{\displaystyle \frac{m_{\mathrm{d}}}{{\left(\mathrm{1}-q_{\mathrm{v}}\right)}^{\mathrm{2}}}}\mathrm{d}{q}_{\mathrm{v}}={\displaystyle \frac{m}{\mathrm{1}-q_{\mathrm{v}}}}\mathrm{d}{q}_{\mathrm{v}},\end{array}\end{array}$$

and the specific humidity after the update is given by

$$\begin{array}{}\text{(A10)}& q_{\mathrm{v}}^{\prime }=q_{\mathrm{v}}+\mathrm{d}{q}_{\mathrm{v}}.\end{array}$$

The implied dry mass of the parcel after the update (which may differ from the parcel's initial dry mass because of numerical errors) is given by

$$\begin{array}{}\text{(A11)}& m^{\prime }={\displaystyle \frac{m_{\mathrm{d}}^{\prime }}{\mathrm{1}-q_{\mathrm{v}}^{\prime }}}.\end{array}$$

Combining Eqs. (A8)–(A11) and solving for $m_{\mathrm{d}}^{\prime }$ gives

$$\begin{array}{}\text{(A12)}& m_{\mathrm{d}}^{\prime }=m_{\mathrm{d}}-{\displaystyle \frac{m}{\mathrm{1}-q_{\mathrm{v}}}}\mathrm{d}{q}_{\mathrm{v}}^{\mathrm{2}}.\end{array}$$

Dry-mass conservation requires $m_{\mathrm{d}}^{\prime }=m_{\mathrm{d}}$, but instead $m_{\mathrm{d}}^{\prime }$ is reduced relative to m_d by an amount that scales with $\mathrm{d}{q}_{\mathrm{v}}^{\mathrm{2}}$.

To obtain stable simulations, we include a mass fixer to enforce conservation of a fixed inventory of dry atmospheric mass. Before each time step, the fixer calculates the current atmospheric dry-mass inventory

$$\begin{array}{}\text{(A13)}& M_{\mathrm{d}}=M-\mathit{\gamma }{M}_{\mathrm{w}},\end{array}$$

where M is the total atmospheric mass, M_w is the mass of water vapor, and γ is the factor used to rescale the water vapor mass source in the pressure tendency calculation. After each time step, the fixer calculates the new dry atmospheric mass

$$\begin{array}{}\text{(A14)}& M_{\mathrm{d}}^{\prime }=M^{\prime }-\mathit{\gamma }{M}_{\mathrm{w}}^{\prime },\end{array}$$

where M^′ and $M_{\mathrm{w}}^{\prime }$ are total and water vapor mass after the time step but before applying corrections. Finally, the fixer multiplies surface pressure by a global multiplicative constant and divides specific humidity by the same constant to set $M_{\mathrm{d}}^{\prime }$ to the dry-mass inventory M_d calculated before the time step while leaving $M_{\mathrm{w}}^{\prime }$ unchanged. Because our model uses sigma coordinates, adjusting the total atmospheric mass by multiplying surface pressure by a single number in all columns changes the mass thickness of each level by an amount proportional to that level's mass thickness. This procedure does not add dry mass in precisely the places where it is lost due to numerical errors, which we view as acceptable because changes in dry mass (proportional to $\mathrm{d}{q}_{\mathrm{v}}^{\mathrm{2}}$) are locally small compared to changes in moist mass (proportional to dq_v), and the mass fixer is used primarily to allow stable long-term simulations that correctly conserve total dry mass.

A4 Lapse rate calculations

In the default version of the model, expressions for lapse rates used to compute reference adiabats used by the simplified Betts–Miller convection scheme differ from lapse rates derived from the governing equations for the resolved flow. In our simulations, we use a modified convection scheme where lapse rates are computed consistently with equations for resolved flow. The inconsistencies corrected by the modified convection scheme are not introduced by the addition of mass sources and sinks, but we document them here for completeness.

Expressions for adiabatic lapse rates in unsaturated and saturated air can be derived from the temperature equation in the model dynamical core,

$$\begin{array}{}\text{(A15)}& {\displaystyle \frac{\mathrm{D}T}{\mathrm{D}t}}={\displaystyle \frac{\mathit{\kappa }{T}_{\mathrm{v}}}{p}}{\displaystyle \frac{\mathrm{D}p}{\mathrm{D}t}}+{\displaystyle \frac{Q}{c_p}};\end{array}$$

the approximate equation used by the model to compute virtual temperature,

$$\begin{array}{}\text{(A16)}& T_{\mathrm{v}}=T(\mathrm{1}+\mathit{ϵ}{q}_{\mathrm{v}})\end{array}$$

with $\mathit{ϵ}=R_{\mathrm{v}}/R_{\mathrm{d}}-\mathrm{1}$; the relationship between specific humidity and temperature changes in the large-scale condensation scheme,

$$\begin{array}{}\text{(A17)}& {\left({\displaystyle \frac{\mathrm{D}T}{\mathrm{D}t}}\right)}_{\text{condensation}}=-{\displaystyle \frac{L_{\mathrm{v}}}{c_p}}{\displaystyle \frac{\mathrm{D}{q}_{\mathrm{v}}}{\mathrm{D}t}};\end{array}$$

the Clausius–Clapeyron equation used to compute saturation vapor pressure,

$$\begin{array}{}\text{(A18)}& e^{*}\left(T\right)=e^{*}\left(T_{\mathrm{0}}\right)\exp \left(-{\displaystyle \frac{L_{\mathrm{v}}}{R_{\mathrm{v}}}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{0}}}}\right)\right);\end{array}$$

the relationship between saturation vapor pressure and saturation specific humidity,

$$\begin{array}{}\text{(A19)}& q_{\mathrm{v}}^{*}(T,p)={\displaystyle \frac{R_{\mathrm{d}}{e}^{*}\left(T\right)}{R_{\mathrm{v}}p-\left(R_{\mathrm{v}}-R_{\mathrm{d}}\right)e^{*}\left(T\right)}};\end{array}$$

and the relationship between saturation vapor pressure and saturation mixing ratio,

$$\begin{array}{}\text{(A20)}& r_{\mathrm{v}}^{*}(T,p)={\displaystyle \frac{R_{\mathrm{d}}{e}^{*}\left(T\right)}{R_{\mathrm{v}}(p-e^{*}(T\left)\right)}}.\end{array}$$

In these equations, T is temperature, $\mathit{\kappa }=R_{\mathrm{d}}/c_p$ is the ratio of the specific gas constant and heat capacity of dry air, Q is a heat source, L_v is the latent heat of vaporization, q_v is specific humidity, R_v is the specific gas constant for water vapor, e^* is saturation vapor pressure, and $q_{\mathrm{v}}^{*}$ is saturation specific humidity.

For an unsaturated adiabat, $Q/c_p=\mathrm{0}$ and q_v is conserved, giving

$$\begin{array}{}\text{(A21)}& {\displaystyle \frac{\mathrm{D}T}{\mathrm{D}\log p}}=\mathit{\kappa }{T}_{\mathrm{v}}\end{array}$$

and

$$\begin{array}{}\text{(A22)}& T\left(p\right)=T\left(p_{\mathrm{0}}\right){\left({\displaystyle \frac{p}{p_{\mathrm{0}}}}\right)}^{\mathit{\kappa }\left(\mathrm{1}+\mathit{ϵ}{q}_{\mathrm{v}}\right)}.\end{array}$$

These expressions replace

$$\begin{array}{}\text{(A23)}& {\displaystyle \frac{\mathrm{D}T}{\mathrm{D}\log p}}=\mathit{\kappa }T\end{array}$$

and

$$\begin{array}{}\text{(A24)}& T\left(p\right)=T\left(p_{\mathrm{0}}\right){\left({\displaystyle \frac{p}{p_{\mathrm{0}}}}\right)}^{\mathit{\kappa }}\end{array}$$

in the calculation of reference temperature profiles below the lifting condensation level (LCL) in the modified version of the convection scheme.

For a saturated adiabat, $Q/c_p={\left(\frac{\mathrm{D}T}{\mathrm{D}t}\right)}_{\text{condensation}}$ and $q_{\mathrm{v}}=q_{\mathrm{v}}^{*}$, and Eqs. (A15)–(A19) give

$$\begin{array}{}\text{(A25)}& {\displaystyle \frac{\mathrm{D}T}{\mathrm{D}\log p}}={\displaystyle \frac{R_{\mathrm{d}}{T}_{\mathrm{v}}-L_{\mathrm{v}}\frac{\partial q_{\mathrm{v}}^{*}}{\partial \log p}}{c_p+L_{\mathrm{v}}\frac{\partial q^{*}}{\partial T}}},\text{(A26)}& {\displaystyle \frac{\partial q_{\mathrm{v}}^{*}}{\partial \log p}}=-\mathit{\alpha }{r}_{\mathrm{v}}^{*},\end{array}$$

and

$$\begin{array}{}\text{(A27)}& {\displaystyle \frac{\partial q_{\mathrm{v}}^{*}}{\partial T}}=\mathit{\alpha }{\displaystyle \frac{L_{\mathrm{v}}{r}_{\mathrm{v}}^{*}}{R_{\mathrm{v}}{T}^{\mathrm{2}}}},\end{array}$$

with

$$\begin{array}{}\text{(A28)}& \mathit{\alpha }={\displaystyle \frac{p(p-e^{*})}{{\left(p-\left(\mathrm{1}-\frac{R_{\mathrm{d}}}{R_{\mathrm{v}}}\right)e^{*}\right)}^{\mathrm{2}}}}.\end{array}$$

Combining gives

$$\begin{array}{}\text{(A29)}& {\displaystyle \frac{\mathrm{D}T}{\mathrm{D}\log p}}={\displaystyle \frac{R_{\mathrm{d}}{T}_{\mathrm{v}}+\mathit{\alpha }{L}_{\mathrm{v}}{r}_{\mathrm{v}}^{*}}{c_p+\mathit{\alpha }\frac{L_{\mathrm{v}}^{\mathrm{2}}{r}_{\mathrm{v}}^{*}}{R_{\mathrm{v}}{T}^{\mathrm{2}}}}}.\end{array}$$

This expression replaces

$$\begin{array}{}\text{(A30)}& {\displaystyle \frac{\mathrm{D}T}{\mathrm{D}\log p}}={\displaystyle \frac{R_{\mathrm{d}}T+L_{\mathrm{v}}{r}_{\mathrm{v}}^{*}}{c_p+\frac{L_{\mathrm{v}}^{\mathrm{2}}{r}_{\mathrm{v}}^{*}}{R_{\mathrm{v}}{T}^{\mathrm{2}}}}}\end{array}$$

in the calculation of reference temperature profiles above the LCL in the modified version of the convection scheme.

We note that the calculation of the LCL itself – done using a lookup table – still assumes that $\mathrm{D}T/\mathrm{D}\log p=\mathit{\kappa }T$ below the LCL. Changing the LCL calculation to use $\mathrm{D}T/\mathrm{D}\log p=\mathit{\kappa }{T}_{\mathrm{v}}$ would have required the current one-dimensional lookup table that estimates LCL pressure to be replaced by a two-dimensional lookup table.

A5 Model verification

Figure A1 shows time series of total atmospheric mass and the relationship between instantaneous precipitation minus evaporation and the surface pressure tendency from mass sources and sinks from global (not TC world) simulations run at a resolution of T42 for the control climate (α=1) without mass sources and sinks (γ=0, no MS), with mass sources and sinks (γ=1, MS), and with mass sources and sinks exaggerated by a factor of 10 (γ=10, 10× MS). All simulations are started from a dry atmosphere with the same initial mass. The total atmospheric mass (as measured by the global-mean surface pressure) remains constant in the simulation without mass sources and sinks, increases slightly in the simulation with mass sources and sinks, and increases by a factor of about 10 more with 10-fold-exaggerated mass sources and sinks (Fig. A1a). After 200–300 d, the atmospheric mass reaches a steady state with no appreciable drift. The surface pressure tendency corresponding to mass sources and sinks is expected from hydrostatic balance to scale with precipitation minus evaporation (P−E) as

$$\begin{array}{}\text{(A31)}& {\left({\displaystyle \frac{\partial p_{\mathrm{s}}}{\partial t}}\right)}_{\text{mass sources and sinks}}=-\mathit{\gamma }g\left(P-E\right)\end{array}$$

with a proportionality factor set by γ and the acceleration from gravity g, and this is confirmed in Fig. A1b.

Figure A1 Time series of global-mean surface pressure (a) and the relationship between precipitation minus evaporation (P−E) and the component of the surface pressure tendency from the mass sink (b) in T42 simulations at α=1. The points in panel (b) are for different grid boxes and instantaneous times at all latitudes, and the diagonal lines in panel (b) show the expected relationship based on hydrostatic balance (see text). Surface pressure tendencies for the 10×-MS simulation are rescaled by a factor of 0.2 in panel (b) to make surface pressure tendencies for MS and no-MS simulations more easily visible.

Appendix B: Relationship between pseudoadiabatic and non-pseudoadiabatic continuity equations

Equation (A1) assumes, consistent with the pseudoadiabatic limit used in our model, that the atmosphere contains only dry air and water vapor. In a more realistic limit, where condensed water is present and parcels gain and lose mass only through the vertical sedimentation of condensed water, mass conservation requires

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\mathit{\rho }\mathit{\delta }x\mathit{\delta }y\mathit{\delta }z=F^{+}\mathit{\delta }x\mathit{\delta }y-F^{-}\mathit{\delta }x\mathit{\delta }y,\end{array}$$

where F⁺ and F⁻ are the sedimentation mass fluxes of condensed water (dimensions of mass per unit area per unit time, defined as positive downward) through the top and bottom of the parcel. Using hydrostatic balance to replace δz with δp gives

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\delta }x\mathit{\delta }y\mathit{\delta }p}}{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\mathit{\delta }x\mathit{\delta }y\mathit{\delta }p=-{\displaystyle \frac{g}{\mathit{\delta }p}}(F^{+}-F^{-}).\end{array}$$

Replacing $\mathrm{D}\mathit{\delta }x/\mathrm{D}t$ with δu, making similar substitutions for y and p, and replacing $(F^{+}-F^{-})/\mathit{\delta }p$ with $\partial F/\partial p$ gives

$$\begin{array}{}\text{(B3)}& \mathrm{\nabla }\cdot \mathbf{v}+{\displaystyle \frac{\partial \mathit{\omega }}{\partial p}}=-g{\displaystyle \frac{\partial F}{\partial p}}.\end{array}$$

The pseudoadiabatic version of the continuity equation (Eq. A4) can be recovered from Eq. (B3) by assuming that timescales for precipitation formation and fallout are infinitely fast. In this limit, the sedimentation flux at pressure level p must be equal to the vertically integrated net condensation at all higher levels because any condensed water produced at higher levels will fall through lower levels within an infinitesimally short time. Because the net condensation rate per unit volume is given by $-{\mathit{\rho }}_{\mathrm{d}}\mathrm{D}{r}_{\mathrm{v}}/\mathrm{D}t$, the sedimentation flux must be

$$\begin{array}{}\text{(B4)}& F\left(p\right)=-\underset{\mathrm{0}}{\overset{p}{\int }}{\mathit{\rho }}_{\mathrm{d}}{\displaystyle \frac{\mathrm{D}{r}_{\mathrm{v}}}{\mathrm{D}t}}{\displaystyle \frac{\mathrm{d}{p}^{\prime }}{\mathit{\rho }g}},\end{array}$$

and differentiating with respect to p gives

$$\begin{array}{}\text{(B5)}& -g{\displaystyle \frac{\partial F}{\partial p}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{d}}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{D}{r}_{\mathrm{v}}}{\mathrm{D}t}}.\end{array}$$

This recovers Eq. (A4) from the right-hand side of Eq. (B3) under the assumption of infinitely fast precipitation formation and fallout.

Appendix C: Continuity equation with mass sources and sinks in potential temperature coordinates

Beginning from Eq. (A3) and using the isentropic density σ to substitute $\mathit{\delta }p=-g\mathit{\sigma }\mathit{\delta }\mathit{\theta }$ gives

$$\begin{array}{}\text{(C1)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\mathit{\delta }x\mathit{\delta }y\mathit{\delta }\mathit{\theta }}}{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\mathit{\sigma }\mathit{\delta }x\mathit{\delta }y\mathit{\delta }\mathit{\theta }={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{d}}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{D}{r}_{\mathrm{v}}}{\mathrm{D}t}}.\end{array}$$

Recall that the material derivative in potential temperature coordinates is

$$\begin{array}{}\text{(C2)}& {\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}={\displaystyle \frac{\partial }{\partial t}}+\mathbf{v}\cdot {\mathrm{\nabla }}_{\mathit{\theta }}+\dot{\mathit{\theta }}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}},\end{array}$$

where v is the horizontal velocity and ∇_θ is the gradient along potential temperature surfaces. Expanding the left-hand-side derivative in Eq. (C1), making substitutions similar to those used to derive Eq. (A4), using the definition of the material derivative in potential temperature coordinates, and replacing r_v with q_v gives

$$\begin{array}{}\text{(C3)}& {\displaystyle \frac{\partial \mathit{\sigma }}{\partial t}}+{\mathrm{\nabla }}_{\mathit{\theta }}\cdot \left(\mathit{\sigma }\mathbf{v}\right)+{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\mathit{\sigma }\dot{\mathit{\theta }}\right)={\displaystyle \frac{\mathit{\sigma }}{\mathrm{1}-q_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{D}{q}_{\mathrm{v}}}{\mathrm{D}t}}.\end{array}$$

Appendix D: Relationship between pseudoadiabatic and non-pseudoadiabatic IPV equations

Beginning from Eq. (B2) and using the isentropic density σ to substitute $\mathit{\delta }p=-g\mathit{\sigma }\mathit{\delta }\mathit{\theta }$ gives

$$\begin{array}{}\text{(D1)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\mathit{\delta }x\mathit{\delta }y\mathit{\delta }\mathit{\theta }}}{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\mathit{\sigma }\mathit{\delta }x\mathit{\delta }y\mathit{\delta }\mathit{\theta }={\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\mathit{\delta }\mathit{\theta }}}(F^{+}-F^{-}).\end{array}$$

Following steps similar to those outlined in Appendix C then gives

$$\begin{array}{}\text{(D2)}& {\displaystyle \frac{\partial \mathit{\sigma }}{\partial t}}+{\mathrm{\nabla }}_{\mathit{\theta }}\cdot \left(\mathit{\sigma }\mathbf{v}\right)+{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\mathit{\sigma }\dot{\mathit{\theta }}\right)={\displaystyle \frac{\partial F}{\partial \mathit{\theta }}}\end{array}$$

as a non-pseudoadiabatic version of the continuity equation in potential temperature coordinates. This is identical to Eq. (2) used to derive the pseudoadiabatic IPV equation but with $\mathit{\sigma }\dot{q_{\mathrm{v}}}$ replaced by $\partial F/\partial \mathit{\theta }$. Following the steps described in Sect. 4.1 with the modified continuity equation gives

$$\begin{array}{}\text{(D3)}& {\displaystyle \frac{\mathrm{D}P}{\mathrm{D}t}}=P\left({\displaystyle \frac{\partial \dot{\mathit{\theta }}}{\partial \mathit{\theta }}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }}}{\displaystyle \frac{\partial F}{\partial \mathit{\theta }}}\right)+{\mathit{\sigma }}^{-\mathrm{1}}\mathbf{k}\cdot \left({\displaystyle \frac{\partial \mathbf{v}}{\partial \mathit{\theta }}}\times {\mathrm{\nabla }}_{\mathit{\theta }}\dot{\mathit{\theta }}\right)\end{array}$$

as a non-pseudoadiabatic version of the IPV equations. Similar to Eq. (21) of Schubert et al. (2001), this version of the IPV equation includes a source term from mass sinks that involves the divergence of the hydrometeor sedimentation mass flux.

The pseudoadiabatic IPV equation (Eq. 6) derived in Sect. 4 can be recovered from Eq. (D3) by assuming that precipitation formation and fallout are infinitely fast. Under this assumption, the divergence of the hydrometeor sedimentation mass flux must be given by Eq. (B5). With the approximation $\mathrm{1}-q_{\mathrm{v}}\approx \mathrm{1}$, used in the derivation in Sect. 4, Eq. (B5) becomes

$$\begin{array}{}\text{(D4)}& -g{\displaystyle \frac{\partial F}{\partial p}}=\dot{q_{\mathrm{v}}}.\end{array}$$

Using $\mathit{\sigma }=-g^{-\mathrm{1}}\partial p/\partial \mathit{\theta }$ gives

$$\begin{array}{}\text{(D5)}& -\mathit{\sigma }g{\displaystyle \frac{\partial F}{\partial p}}={\displaystyle \frac{\partial F}{\partial \mathit{\theta }}}=\mathit{\sigma }\dot{q_{\mathrm{v}}}.\end{array}$$

Substituting $\partial F/\partial \mathit{\theta }=\mathit{\sigma }\dot{q_{\mathrm{v}}}$ into the non-pseudoadiabatic IPV equation (Eq. D3) recovers the pseudoadiabatic IPV equation (Eq. 6) used in Sect. 4.

Appendix E: Impact of precipitation mass sinks in TC world simulations without hyperdiffusion changes

The use of ∇⁴ hyperdiffusion and increased hyperdiffusion coefficients is not essential for producing tropical cyclones in TC world simulations. Quarter-Earth-radius TC world simulations with “default” hyperdiffusion parameters (∇⁸ hyperdiffusion and no increases to hyperdiffusion coefficients relative to conventional T170 simulations) produce TCs with slightly weaker central pressure anomalies compared to quarter-Earth-radius simulations with modified hyperdiffusion (compare Fig. E1a with Fig. 10e) and provide no evidence for systematic TC deepening when mass sinks are added (Fig. E1a). Unlike simulations with modified hyperdiffusion, simulations with default hyperdiffusion parameters suggest that including mass sinks may produce simulations with slightly higher TC counts and slightly weaker average central pressure anomalies (compare averages shown by vertical bars in Fig. E1), possibly indicating that mass sinks help maintain relatively weak TCs. This effect, if real, would be interesting and worthy of further investigation but is distinct from the focus of our work, which is on whether mass sinks produce systematically stronger TCs.

Figure E1 Time series of average TC central pressure anomalies (a) and TC count (b) in quarter-Earth-radius TC world simulations without hyperdiffusion increases. The vertical bars at the end of the time series indicate the range of time averages taken over four 100 d blocks between days 200 and 600.