## Abstract

The present study analyzes the spatial and temporal variability of zonally integrated meridional atmospheric heat transport due to transient eddies in a hierarchy of datasets. These include a highly idealized two-layer model seeded with point geostrophic vortices, an intermediate complexity GCM, and the European Centre for Medium-Range Weather Forecasts interim reanalysis (ERA-Interim) data. The domain of interest is the extratropics. Both the two-layer model and the GCM display a pronounced temporal variability in the zonally integrated meridional transport, with the largest values (or pulses) of zonally integrated transport being associated with extended regions of anomalously strong local heat transport. In the two-layer model these large-scale coherent transport regions, termed “heat transport bands,” are linked to densely packed baroclinic vortex pairs and can be diagnosed as low-wavenumber streamfunction anomalies. In the GCM they are associated with both the warm and cold sectors of midlatitude weather systems. Both these features are also found in ERA-Interim: the heat transport bands match weather systems and occur primarily in the storm-track regions, which in turn correspond to planetary-scale climatological streamfunction anomalies. The authors hypothesize that the temporal variability of the zonally integrated heat transport is partly linked to oscillatory variations in the storm-track activity but also contains a background red noise component. The existence of a pronounced variability in the zonally integrated meridional heat transport can have important consequences for the interplay between midlatitude dynamics and the energy balance of the high latitudes.

## 1. Introduction

Meridional heat transport by the atmosphere and oceans has been the subject of a vast body of literature and has been long acknowledged as crucial to both the time-mean structure (e.g., Budyko 1969; Sellers 1969; Stone 1978) and time-varying component (e.g., Bjerknes 1964) of Earth’s climate. More recent efforts have used complex GCMs in order to diagnose the response of heat transport to anthropogenic forcing (e.g., Zelinka and Hartmann 2012) and its role in polar amplification (e.g., Pithan and Mauritsen 2014).

The present study builds upon the smaller number of papers that have highlighted the fundamental statistical and dynamical features of the heat transport. Swanson and Pierrehumbert (1997) first noted that atmospheric heat transport by transient motions can display extreme events, orders of magnitude larger than the typical climatological values. Messori and Czaja (2013, 2014) extended this analysis concluding that, at any given location in the extratropics, net seasonal heat transport is effectively set by a few extreme days every season. Similar conclusions hold when considering zonally integrated meridional heat transport, whose probability density function (PDF) shows large infrequent transport pulses (Messori and Czaja 2015, hereafter MC15). These pulses are not characterized by a uniform increase in transport around the full latitude circle but rather by regions of locally enhanced transport. Crucially, a systematic understanding of the characteristics and drivers of these pulses is still lacking.

A number of idealized models, aiming to reproduce the large temporal variability of the transport, have been proposed. These range from Swanson and Pierrehumbert’s (1997) Lagrangian advection to Messori and Czaja’s (2013) simple sinusoidal wave superposition to Ambaum and Novak’s (2014) storm-track oscillator. The simplicity of these models suggests that the sporadic nature of heat transport is a very fundamental feature of our atmosphere. However, these efforts have focused on the transport at fixed locations, without addressing its spatial characteristics.

There is therefore a clear knowledge gap in our understanding of zonally integrated heat transport pulses in the atmosphere, both in terms of conceptual models and of identifying the underlying dynamical features in the real atmosphere. In the present study, we aim to address this gap by first characterizing the pulses in a highly idealized model and then using more realistic datasets to link them to large-scale atmospheric features. This both complements and extends the previous analyses and provides a coherent overview of the variability of zonally integrated meridional heat transport in the atmosphere. The hierarchy of models/datasets and the methodology we use are described in sections 2 and 3. The starting point is a quasigeostrophic potential vorticity (PV) model, based on the two-layer, point-vortices heton concept (Hogg and Stommel 1985a,b). This environment, although highly idealized, allows the simultaneous investigation of both the temporal sporadicity and spatial variability of the transport (section 4). The second part of the analysis focuses on an intermediate complexity GCM run in an aquaplanet configuration (section 5). This more realistic framework enables a better understanding of the dynamical features underlying the zonally integrated heat transport pulses. Finally, we investigate the pulses in the ERA-Interim dataset (section 6). In section 7 we highlight the strengths and limitations of the idealized setups and relate the insights they afford to the features identified in the reanalysis. Concluding remarks are provided in section 8.

## 2. Models and data

### a. The heton model

The first and most idealized model we discuss (section 4) is a setup based on Hogg and Stommel’s heton concept (Hogg and Stommel 1985a). The model is a two-layer system, seeded with point geostrophic vortices. These vortices deform the interface between the two layers and are therefore associated with both a thermal anomaly and a streamfunction. The magnitude of these two quantities depends on the strength of the vortices. Two vortices in different layers of the model with PV of opposite sign can form a baroclinic pair and effect a heat transport; such pair is therefore termed a heton. Depending on the sign of the PV anomalies, hetons can be either hot or cold (see Fig. A1). The model and its dynamics are described in detail in a number of studies (e.g., Hogg and Stommel 1985a; Legg and Marshall 1993). Further details are provided in appendix A.

If more than one heton is seeded in the same domain, it will interact with all the neighboring hetons. As a result of these interactions, the vortices may split from their original pair and recouple with other vortices in the same or in a different layer. In the far field, hetons with same temperature repel; this is not true of the near field, and clusters of vortices can therefore form. These effectively behave like a composite vortex. Since the strength of the hetonic interactions, and hence the resulting velocity fields, are proportional to the intensity of the vortices themselves, composite vortices can move very rapidly across the model domain. This behavior was first described by Hogg and Stommel (1985b), who called it hetonic explosion.

Since the heton model is based on a finite number of discrete vortices rather than a continuous PV distribution, the meridional heat transport can be computed from vortex motions. The quantized nature of the model provides a stimulating analogy to the real atmosphere where, as discussed in the introduction, the transport is highly discontinuous and fundamentally sporadic in nature. The integrated hetonic heat transport across a closed boundary is found as a function of the number of vortices crossing said boundary, the strength of the vortices, and the rate of seeding of new vortices within the domain (see appendix B). A limitation of this technique is that it does not allow for a local heat transport to be computed.

The experiment described here considers a circular domain, meant to represent an idealized polar region (see also appendix C). Cold hetons are seeded within the domain at every time step at random locations to reproduce the thermal effect of spatially inhomogeneous outgoing longwave radiation [for local and mesoscale processes affecting this, see, for example, Ramanathan et al. (1989) and Raval et al. (1994)]. The vortices are then free to interact and cross the domain boundary. The intensity of the vortices decays over time, analogously to the relaxation of a thermal anomaly. The domain size is chosen to roughly coincide with a spherical cap bounded by the 60° latitude circle. The whole setup approximates to a large, distributed baroclinic vortex, much like the original heton experiments performed by Hogg and Stommel (1985b). The hope is to capture the large-scale statistical features of the heat transport witnessed in the real atmosphere in terms of hetonic heat transport into the cap. Heat transport is computed as the zonally integrated transport across the seeding domain boundary and across circular boundaries with various other radii. Distances in the model are expressed in terms of Rossby radii of deformation *R*_{D}, where one radius of deformation is approximately equal to the typical length scale for vortex interactions (see Table C1). In this framework, the seeding domain has a radius of 4 deformation lengths. The seeding rate, strength, and decay rate of the vortices are selected such that the mean transport into the seeding domain is close to real-world values. In the model, such transport is calculated to be just above 3.5 PW; as point of comparison, the annual-mean net meridional atmospheric heat transport across 60°N is of order 3 PW (Fasullo and Trenberth 2008). Further details on the experimental setup and a brief discussion of the sensitivity tests performed are provided in appendix C. The analysis focusses on a 1-month-long equilibrium simulation; this is sufficiently long to provide smooth distributions of the transport while being short enough to investigate in detail the zonally integrated transport extremes in section 4b.

### b. The aquaplanet GCM

The intermediate complexity climate model used is the MITgcm (Marshall et al. 1997a,b). The data are analyzed after regridding onto a regular 2° × 2° latitude–longitude grid. The model is run as a hydrostatic aquaplanet atmosphere-only model on 25 equally spaced pressure levels. With the exception of the radiation scheme, the physics module is based on Frierson et al. (2006) and Frierson (2007) and includes the amendments made in O’Gorman and Schneider (2008). In the radiation scheme, optical depths are calculated from the specific humidity field using the three spectral region radiation scheme described in Geen et al. (2016). In keeping with the idealized setup of the simulation, the sea surface temperatures are fixed to the control profile from the Aqua-Plant Experiment (APE) Atlas (Neale and Hoskins 2000). A perpetual-equinox top-of-atmosphere insolation, with no seasonal or diurnal cycle, is used (O’Gorman and Schneider 2008).

Ten years of daily (1200 UTC) model data are analyzed. In view of the absence of a seasonal cycle, all days of the year are taken into consideration. Furthermore, no distinction is made between the two hemispheres, since the simulation’s setup is entirely symmetric.

### c. The reanalysis data

The present paper utilizes the European Centre for Medium-Range Weather Forecasts interim reanalysis (ERA-Interim) product (Dee et al. 2011). Daily (1200 UTC) fields with a resolution of 0.7° latitude × 0.7° longitude are considered. MC15 also analyzed 6-hourly data and found that the spatial and temporal variability of the transport can be adequately characterized with daily values. The period taken into consideration includes 22 December–February (DJF) and 22 June–August (JJA) time series, from June 1989 to February 2011. This period is selected to ease the comparison with previous analyses. In the interest of conciseness, unless otherwise specified the four season–hemisphere combinations are analyzed together. Past studies have shown that they all display qualitatively similar heat transport variability characteristics (e.g., Messori and Czaja 2013).

## 3. Methodology

### a. Meridional heat transport by transient motions

The analysis of the heat transport in the GCM and the reanalysis focuses on the 850-hPa level. This is the level of the peak heat transport by transient eddies, and previous analyses found that it provides a good indication of the overall transport statistics (Messori and Czaja 2013). Figures only refer to the 850-hPa level; for completeness, tables also provide upper-level values at 300 hPa.

The transient-eddy transport is computed as the product of meridional velocity *υ* and moist static energy *H* (herein also termed MSE) temporal anomalies (e.g., Neelin and Held 1987). These are defined as departures from the linearly detrended seasonal mean and are denoted by a prime. Velocity (and hence heat transport) is positive poleward in both hemispheres. The Northern Hemisphere (NH) and Southern Hemisphere (SH) domains considered in the analysis span 30°–89°N and 30°–89°S, respectively (87°N and 87°S in the aquaplanet). For ease of reference, throughout the paper the transient-eddy transport is simply called heat transport.

In the GCM and the reanalysis, the term local event refers to the transport value at a single grid point. The zonally integrated transport is computed by integrating local *υ*′*H*′ values over the width of their respective grid boxes and then summing all values around a given latitude. The transport is then normalized by layer thickness and is expressed in watts per 1000 hPa—namely, the meridional transport in watts that would occur if the flux were realized at all vertical levels. The values obtained can also be directly compared to the results of the heton model, which are expressed in petawatts and can be interpreted as a vertically integrated meridional heat transport.

In all three models, extreme events (termed local extremes for local transport and zonal extremes for zonally integrated transport) are defined as values of *υ*′*H*′ that exceed the 95th percentiles of the respective distributions for the full domain and time period considered. MC15 have shown that the qualitative results are largely independent of the exact percentile chosen as threshold.

For the GCM and the reanalysis, this study also discusses heat transport bands. A band at a given latitude is identified as a continuous zonal region displaying (i) meridional heat transport values above the median value of the local heat transport distribution and (ii) at least one local extreme. This selection process is illustrated schematically in Fig. 1. A band lasts for more than 1 day if, within one latitude grid box, a second band is identified that matches more than half of the longitudes included in the original band. The spatial extent of the bands is quantified in terms of wavenumbers. A band would represent the positive half of the wavelength of an ideal, perfectly symmetric wave. The wavenumbers discussed in the text are therefore scaled by a factor of 2 relative to the extent of the band. This 1:2 scaling does not translate perfectly to a real-world signal but provides nonetheless an indication of the spatial scale of the events described here. In the heton model, local heat transport and, hence, local extremes cannot be quantified numerically. However, the vortices can travel both in pairs/small groups or a cluster and in large, densely packed arrays. If one interprets a single vortex pair moving out of the domain as a large local thermal anomaly, these different behaviors therefore allow for well-separated and banded local events to be quantified in a qualitative fashion.

Finally, we discuss the most likely value (MLV) of PDFs. This is defined as the central value of the bin with the highest frequency of events. The exact value of the MLV clearly depends on the choice of bins, but it is nonetheless a robust indicator if used to compare the order of magnitude of the most frequent value of a variable to the one of its most extreme realizations (MC15).

### b. Atmospheric fronts

To provide a physical–dynamical interpretation of the heat transport extremes discussed in the text, we identify atmospheric fronts using a metric introduced in Sheldon (2015) and later adopted in Geen (2016) and Parfitt et al. (2016). Front intensity is defined as

where *ζ*_{925} is isobaric relative vorticity and *|*∇*T*_{925}*|* is the isobaric temperature gradient. *F*, the nondimensional counterpart of *F**, is then obtained through division by typical scales for temperature gradient [1 K (100 km)^{−1}] and vorticity (10^{−4} s^{−1}, a reference Coriolis parameter value at midlatitudes). A front is then identified at a given location if *F* ≥ 0.1 (aquaplanet) or *F* ≥ 1 (ERA-Interim). These thresholds are chosen to provide spatial and frequency characteristics similar to those of alternative frontal detection methods (Geen 2016; Parfitt et al. 2016).

### c. Atmospheric rivers

In the discussion of our results in section 7, we relate local heat transport extremes to atmospheric rivers (ARs). These are long narrow bands of large integrated water vapor flux and account for a large part of the meridional moisture transport in the midlatitudes (Zhu and Newell 1998). Here, we identify ARs following the methodology of Zhu and Newell (1998), recently adopted also in Baggett et al. (2016). Starting from the horizontal vertically integrated water vapor flux **Q,** we define the flux associated with ARs as

Here, *Q*_{t} is a threshold value defined as

where is simply the zonal average of *Q* at each latitude, while *Q*_{MAX} is the maximum value of *Q* at each latitude. The prefactor of 0.3 has been chosen after extensive testing in previous studies to best discriminate the filamentary ARs from the climatological fluxes (Zhu and Newell 1998; Baggett et al. 2016).

### d. Wavelet spectra

Part of the analysis looks at a Morlet wavelet transform of zonally integrated meridional heat transport. Wavelets provide a time–frequency picture of the spectral power of a time series while the global wavelet spectrum, which is obtained by combining the spectral power at all time steps in the dataset, provides very similar results to the more traditional Fourier spectrum. Following Grinsted et al. (2004), we normalize the heat transport time series to have zero mean and unit standard deviation prior to applying the transform. The significance level of the wavelet power spectrum is computed using a chi-squared test, while the significance of the global spectrum follows Torrence and Compo’s time-averaged test (Torrence and Compo 1998). For further details on the properties of the Morlet wavelet and on wavelet transforms in general, the reader is referred to Goupillaud et al. (1984), Torrence and Compo (1998), and Grinsted et al. (2004).

## 4. A minimal physics model of heat transport

We begin by analyzing the results of the heton experiment. As discussed in section 2a, hetonic heat transport is computed as the zonally integrated transport into the domain of choice and cannot be related to transport values at a single location. We aim to address two basic questions:

Temporal variability: Do the top percentiles of the distribution of model heat transport (here termed extremes) account for significant portions of the overall integral of the distribution (section 4a)?

Spatial variability: If so, do they systematically match specific dynamical features in the model (section 4b)?

### a. The temporal variability of hetonic transport

Figure 2a displays the time series of the zonally integrated hetonic transport. The horizontal dashed line marks the 95th percentile of the curve. The very complex dynamics that govern large numbers of vortices lead to a highly variable heat transport, with several well-pronounced peaks. To determine whether these peaks play a significant role, we compute the contribution of the top five percentiles of the heat transport distribution to the net meridional transport (Table 1). Note that the transport is computed across circular domains with different boundaries, starting from the size of the heton seeding domain (4*R*_{D}). The fraction of the transport accounted for by the extremes is seen to generally increase with increasing radius, although this trend is reversed for radii beyond twice that of the seeding domain.

It is instructive to compare the values shown in Table 1 to the weight of the corresponding events in a Gaussian distribution. To obtain these values, Gaussian distributions with the same means and standard deviations as the model data are constructed. The top five percentiles of the Gaussian curves are then selected, and their weight relative to the integral of the Gaussians is computed as a percentage. The resulting values are shown in the second data column of Table 1. Even though the values in the two columns appear to be similar, a random Monte Carlo sampling procedure shows that they are statistically different at the 99% confidence level for all radii considered. Note that, if instead of computing a percentile threshold for the Gaussian distributions the same numerical threshold found for the hetonic distributions is applied, the weight of the extreme events in the Gaussian case is between 1% and 6% (Table 1).

The important role of the upper percentiles can clearly be seen in the PDF of the hetonic transport (Fig. 2b). In the figure, the transport is computed across a domain boundary of 7*R*_{D}, where the simulation presents a large contribution from the extremes. This radius was chosen, even though it is not the one with the largest extreme event contribution for the simulation discussed here, because it has the highest average contribution across the different sensitivity tests performed (see appendix C). The PDF reflects the transport’s high variability and displays a well-defined MLV (marked by the continuous vertical line) and a long positive tail leading to a large positive skewness (2.2). These results provide a basis for calling the top percentiles of the hetonic transport distributions extremes.

### b. The spatial variability of hetonic transport

We next link the temporal variability highlighted in section 4a to specific dynamical features in the heton model. We find that the transport’s variability is largely determined by the number of vortices crossing the domain boundary (see appendix B), pointing to heat transport extremes as time steps at which large numbers of vortices leave the domain. This could correspond to the movement of either a large number of individual vortices or to one or more large vortex clusters. The latter case would be analogous to the hetonic explosions described in section 2a. To test whether hetonic explosions correspond to heat transport extremes, we look at snapshots of the model domain at the time steps during which the heat transport peaks. These are marked by numbers in Fig. 2a. In the interest of conciseness, only the largest peak, namely number 5, will be discussed in detail here.

The position of the vortices and the corresponding upper- and lower-layer streamfunctions for the present case study are shown in Figs. 3a–c, respectively. There is a dense cluster of vortices, marked by the thick ellipse, which detaches itself from the bulk of the other vortices. The individual vortices within the cluster move synchronously, effectively behaving as a composite vortex. The streamfunction (Figs. 3b,c) highlights this coherent behavior. The flow initially displays a wavenumber 3 disturbance. One of the troughs then begins to grow, extends, and finally detaches itself from the main circulation pattern and exits the domain. This feature is found to drive a transport of roughly 6.0 PW.^{1} The local maximum of the heat transport extreme is ~7.5 PW. The hetonic explosion is therefore the primary driver of the heat transport extreme, accounting for approximately 80% of the time step’s net heat transport value. An analysis of the other transport events marked in Fig. 2a comes to very similar conclusions (not shown).

The large heat transport values are systematically associated with hetonic explosions, which evolve from low-wavenumber disturbances developing along the rim of the seeding domain. The link between hetonic explosions and heat transport extremes also explains why the weight of the extremes increases when the radius of the heat transport domain is increased. The vortex clusters often form outside the seeding domain, meaning that large heat transport events associated with hetonic explosions are best captured at large radii. We conclude that the heton model highlights both the intrinsically sporadic nature of the heat transport’s temporal variability and the important role of coherent vortex clusters in driving zonally integrated transport extremes.

## 5. Heat transport in an aquaplanet

We begin by showing that the intermediate complexity GCM described in section 2b), also reproduces a temporally sporadic zonally integrated transport. Next, we take advantage of the fact that the GCM allows us to compute both local and zonally integrated transport to investigate the following:

the link between local and zonal extreme event frequency (section 5a),

the existence of geographically extensive regions of above-average transport and their role in driving zonally integrated extremes (section 5b), and

the dynamical origins of these features (section 5c).

### a. The zonal-mean view and local extremes

We compute the zonally integrated meridional heat transport as described in section 3a. In constructing the PDFs, the integral of the heat transport around a given latitude circle, on a given day, is treated as a single data point. Figure 4a displays the PDF for the transport over the range 30°–87°N and 30°–87°S. The distribution displays a prominent MLV at near-zero values and a long, thin positive tail, which leads to a skewness of 0.91. The contributions of the top percentiles to both the poleward-only and the net seasonal transports—respectively taken to be the integral of the positive portion and of the whole transport distribution for the full spatial and temporal domains being analyzed—are displayed in Table 2. As was the case for the heton model, these values are significantly higher than what would be expected for a Gaussian variable.

To understand the origin of these zonal pulses we compare the frequency of local extremes on zonal extreme days versus all other days. Virtually all zonal extremes correspond to days/latitudes displaying at least one local extreme, while 49% of zonally nonextreme days/latitudes display no local extremes. The PDFs in Fig. 4b display the number of local extremes around a full latitude circle on days/latitudes corresponding to zonal extremes (white bars) and on all other days (gray bars), excluding days/latitudes with no local extremes. The upper bound on the number of local extremes that can occur at a given latitude is given by the number of grid boxes around each latitude circle—namely, 180. PDF for the extreme zonal days peaks at larger values than that for all other days, suggesting that the zonal extremes are typically associated with a heightened frequency of local extremes. The two distributions are statistically different under a two-sample Kolmogorov–Smirnov test, with the null hypothesis of the same parent distribution rejected at the 5% significance level.

### b. The spatial variability of meridional heat transport

The upper tails of both distributions in Fig. 4b display very large numbers of local extremes co-occurring at a given latitude. This suggests that local extremes might band together, leading to extensive regions of above-average transport. To quantify this behavior, we analyze heat transport bands, as defined in section 3a and Fig. 1. We find that 54% of all latitude–day pairs display at least one band. Figure 5 displays the PDFs of the number of bands per latitude per day (Fig. 5a), the zonal extent of the bands in wavenumbers (Fig. 5b),^{2} the duration of each band (Fig. 5c), the net meridional heat transport associated with each band (Fig. 5d) and all bands (Fig. 5e) occurring on the same day/latitude, the fraction of the zonally integrated transport for the day/latitude at which each band occurs associated with the band itself (Fig. 5f) and all bands (Fig. 5g) occurring on the same day/latitude, and the fraction of grid boxes belonging to a band (Fig. 5h). All plots exclude days/latitudes with no bands; fractional contributions to the transport also exclude days/latitudes with equatorward transport values. The mean values of these quantities are reported in Table 3. Most of the days/latitudes display between one and five bands, whose longitudinal extent is typically between wavenumbers 6 and 10 and whose duration is typically between 1 and 4 days. Individually, these bands often account for a large portion of the zonally integrated meridional heat transport at the latitude/day on which they occur. If taken cumulatively, the bands occurring on a given day/latitude often approach and even exceed the net zonally integrated transport. Fractional contributions greater than 1 can be explained by the fact that it is not unusual for the local transport to be equatorward (and hence negative in the sign convention used here), while the bands by definition only include grid boxes displaying positive transport values. For both levels, the cumulative fractional contribution of bands to the heat transport exceeds the fraction of grid boxes within a band by a factor of more than 3 (Table 3). This confirms that our definition of bands captures extensive areas of heightened meridional heat transport. We note that the very high fractional contributions found at 300 hPa are mainly driven by a limited number of latitude–day pairs where the net atmospheric heat transport is close to zero.

We next investigate the role of the heat transport bands described above in driving the zonal extremes (Table 3). In terms of the individual contribution, bands occurring on zonal extreme days typically carry more heat poleward than bands occurring during other days; however, since zonal extremes display unusually large net transport values, the fractional contribution of the single bands to the net transport is actually lower. Nonetheless, zonal extremes display on average almost twice the number of bands as nonextreme days, leading to a higher cumulative fractional contribution approaching 1 at 850 hPa.

### c. Dynamical drivers of meridional heat transport bands

The above analysis suggests that the increase in heat transport witnessed on zonal extreme days relative to the climatology can be largely ascribed to banded grid boxes and is driven by an increase in both the number and the magnitude of the bands. Understanding the dynamical origin of these bands can therefore point to the drivers of the sporadic nature of the zonally integrated meridional heat transport. In particular, it is important to understand whether these bands can be likened to large-scale coherent transport structures.

We begin by verifying the sign of the meridional velocity and MSE anomalies driving the bands. A priori, a positive heat transport value might be equally driven by the product of two negative or two positive anomalies. In physical terms, these would correspond to an equatorward wind anomaly associated with unusually cold, dry air or a poleward wind anomaly associated with unusually warm, moist air, respectively. Figure 6a displays the fraction of banded grid boxes associated with a poleward meridional velocity anomaly as a function of latitude. Since the bands only include poleward transport values, this is identical to plotting the fraction of positive *H*′ values. The fractions range from around 55%–60% in the midlatitudes to 100% in the high latitudes. The roughly 60%–40% split seen over a large portion of the domain is consistent with the bulk of the heat transport being associated with the circumglobal storm tracks of the aquaplanet. If this interpretation is correct, the positive *υ*′ cases would correspond to heat transport in the warm sector of the extratropical storms, while the negative *υ*′ would correspond to transport in the cold sector. To test the plausibility of this hypothesis, Figs. 6c and 6d display snapshots of *H*′ and *υ*′*H*′, respectively, at a randomly selected time step. The footprint of the storm track is immediately evident and, as suggested by the fractions discussed above, the cold sectors of the storms effect a sizeable heat transport, only slightly weaker than that associated with the warm sectors. The blue contours mark heat transport bands and can be seen to correspond to both positive and negative *H*′ values. The fronts, which we identify as described in section 3b, are marked by the dark red contours. As expected, they typically lie along the trailing edge of the positive *H*′ anomalies and, hence, at the edge of heat transport bands. Indeed, the PDF of *F* in regions within heat transport bands is significantly shifted toward higher values compared to the distribution of *F* outside transport bands (Fig. 6b), and the two distributions are different at the 5% significance level according to a two-sample Kolmogorov–Smirnov test.

The wavenumber spectra of *υ*′, *H*′, and the transport, displayed in Fig. 7, confirm this pattern. In the midlatitudes, where the storm track is located, the spectra of *υ*′ and *H*′ are dominated by wavenumber 4–7 contributions (Fig. 6c shows seven complete waves over a full latitude circle), while the transport spectrum displays a broadening of the spectral power, with a clear peak around wavenumber 10. This is likely associated with near in-phase *υ*′ and *H*′, since the product of two perfectly in-phase sinusoidal-type variables would have double the wavenumber of its components. We note that the transport also displays significant spectral power at low wavenumbers, between 1 and 5. A qualitatively similar spectrum is seen at 300 hPa (not shown).

## 6. Heat transport extremes in the reanalysis data

The idealized models discussed in sections 4 and 5 point to a strong variability of the zonally integrated meridional heat transport, associated with precise dynamical features. Zonally integrated meridional heat transport in the reanalysis data also shows large, infrequent transport pulses. We briefly describe the temporal variability of the zonally integrated transport in section 6a (the reader is referred to MC15 for a more complete analysis). Next, we focus on the associated transport bands and the underlying dynamical features, which have yet to be understood (section 6b).

### a. The zonally integrated heat transport variability and heat transport bands

The PDF of the zonally integrated atmospheric heat transport in the reanalysisis is reminiscent of the features seen in the aquaplanet distribution (cf. Figs. D1a and 4a). The contributions of zonal extreme days to the net seasonal transport at 850 hPa are again well above what would be expected of a Gaussian variable and in good qualitative agreement with the results for the aquaplanet GCM (cf. Tables D1 and 2). As in the aquaplanet, the zonal extremes correspond to unusually high numbers of local extremes (not shown), pointing to the latter as key drivers of large zonally integrated transport values. We note that the range of latitudes over which we composite the heat transport covers regions with very different dynamics. Nonetheless, highly skewed distributions for both the local (not shown) and zonally integrated heat transport (see Figs. D1c,d) are recovered when specific latitude bands are selected, indicating that the fundamental features of the transport’s variability investigated here are not dependent on the choice of domain.

Next, we repeat the analysis performed on the aquaplanet data in section 5b above. We find that 72% of all latitude–day pairs display at least one band. Figure 8 displays the PDFs of several band metrics and can be compared to the aquaplanet results shown in Fig. 5. Table 4 provides the mean values for different pressure levels and different subsets of days. Most of the days/latitudes display between one and five bands, whose longitudinal extent is typically between wavenumber 5 and 10. The bands mostly last for less than 3 days, slightly below the values seen in the aquaplanet, and generally account for a large portion of the zonally integrated meridional heat transport at the latitude/day on which they occur. If taken cumulatively, the bands occurring on a given day/latitude often exceed the net zonally integrated transport. The mean fractional values for the 300-hPa level are again very sensitive to a limited number of latitude–day pairs where the net atmospheric heat transport is close to zero. A comparison of Figs. 8g and 8h highlights that the cumulative fractional contribution of bands to the heat transport systematically exceeds the fraction of grid boxes within a band. This confirms that the bands match extensive areas of heightened meridional heat transport.

As in the aquaplanet, the heat transport bands described above provide a strong contribution to the zonal extremes. Bands are more frequent and transport more heat on zonal extreme days than during other days. Since zonal extremes display unusually large net transport values, the cumulative fractional contribution of the bands on these days is slightly below the all-day average (Table 4). Nonetheless, anomalous heat transport at banded grid boxes accounts, on average, for 96% of the additional net heat transport seen on extreme zonal day/latitudes relative to all other day–latitude pairs.

### b. Large-scale coherent transport regions

The above analysis confirms the inferences drawn from the idealized frameworks discussed in sections 4 and 5—namely, that the sporadic nature of meridional heat transport is largely driven by extensive regions of heightened transport. Here we investigate in greater detail the atmospheric phenomena associated with these regions in the reanalysis data.

Figure 9a displays the fraction of banded grid boxes associated with a poleward meridional velocity anomaly as a function of latitude. Up to almost 70°N, the frequency of positive velocity anomalies exceeds that of negative ones by no more than 20%. In the aquaplanet, the approximate 60%–40% split in the midlatitudes was associated with the partitioning of heat transport between the warm and cold sectors of midlatitude weather systems; the same interpretation holds for the reanalysis. Figures 9c and 9d display snapshots of *H*′ and *υ*′*H*′, respectively. As was seen in the aquaplanet model, fronts typically lie on the edges of the heat transport bands. While the majority of the bands lie in or in the proximity of storm-track regions, it is not uncommon to see large bands over the continental landmasses. The PDFs of the frontal metric *F* both within (white) and outside (gray) the bands are shown in Fig. 9b. Banded grid boxes display systematically larger values of *F* than all other grid boxes, and the two distributions are different at the 5% significance level according to a two-sample Kolmogorov–Smirnov test.

Figures 10a and 10b present *υ*′ (contours) and *H*′ (colors) composites at each latitude for all bands associated with positive and negative anomalies, respectively. All bands are centered on the local meridional heat transport maximum, and the *x* axis displays the longitudinal distance from such maximum. As the figure clearly illustrates, both the heat transport maximum and the rest of the band are typically driven by in-phase *υ*′ and *H*′ anomalies. The sign of the anomalies typically reverses beyond the edges of the band, which is consistent with the picture of the boundaries of the bands being marked by frontal systems. The corresponding standard deviations, shown in Fig. D2, suggest that the composites provide a good representation of the typical structure of a band.

In the midlatitudes, where the storm tracks are located, the spectra of *υ*′ and *H*′ are dominated by wavenumber 4–7 contributions, while the transport spectrum displays a broad range of wavenumbers with high spectral power (Fig. 11). The peak around wavenumbers 5–9 is consistent with the zonal extent of the heat transport bands shown in Fig. 8b.

## 7. Discussion and physical interpretation

### a. A hierarchy of models

The highly idealized two-layer model discussed in section 4 considers pairs of point vortices—the hetons—seeded in a circular domain. Individual hetons are baroclinic pairs, which shear the two layers over a radius *R*_{D}. Vortices can cluster, causing a decrease in the system’s potential energy and creating preferred locations for motion away from the model’s seeding domain. Large groups of vortices moving synchronously away from the seeding domain, termed hetonic explosions, result in zonally integrated heat transport extremes. These can also be seen as preferred locations for local conversion of baroclinicity to eddy heat transport, with the caveat that in the heton model we can only compute the zonally integrated meridional heat transport due to any given vortex pair or vortex cluster. From a large-scale perspective, these vortex clusters emerge as streamfunction anomalies developing along the rim of the heton seeding domain. Since the heton model is an equilibrated system, these instabilities span a continuous range of spatial scales. The largest explosions, such as that illustrated in Fig. 3, result in wavenumber 2 or 3 deformations, but higher wavenumber instabilities are also observed (Pedlosky 1985; Helfrich and Send 1988). Pedlosky (1985) also noted that, under the assumption of an initially uniform heton cloud with periodic perturbations with zero phase speed, the radial displacement of the cloud in the top and bottom layers is equal in magnitude and opposite in sign. In other words, since the rim current is a baroclinic flow the vortex clusters in the top and bottom layers are out of phase. We therefore expect a strong vertical shear in the initial stages of a hetonic explosion, as is indeed seen in Fig. 3. As a further caveat, we note that the strongest hetonic explosions largely overcome the climatological zonal flow, while in the real atmosphere the zonal jet can meander and undergo large deviations from the climatology but is almost always clearly visible.

This picture, although highly idealized, is consistent with the results from the ERA-Interim data since it is possible to interpret the zonally integrated heat transport extremes in the reanalysis as being associated with a range of spatial scales, spanning from planetary to synoptic scales (Fig. 11c). Figures 12a and 12b display the ERA-Interim eddy (colors) and total (black contours) streamfunction, averaged over 625–200 and 1000–625 hPa, respectively, for a randomly selected day displaying a zonal extreme in the NH midlatitudes. This is intended to mimic the heton model’s two-layer setup. The heat transport bands at 850 hPa are marked by the blue contours. Figures 12c and 12d show an individual band within the region marked by the black box in Fig. 12b. In the lower layer, it is located on the leading edge of a deep streamfunction trough. The same location corresponds to a crest in the upper layer due to the streamfunction’s westward phase tilt, much like the shear induced by vortex clusters in the two-layer model. At a planetary scale, the streamfunction displays a zonal wavenumber ~3 structure, and there is a region between the eastern seaboard of North America and continental Europe where the streamfunction peaks and troughs are most prominent and heat transport bands cluster (Figs. 12a,b). The latter feature is not unique to the chosen snapshot: extreme zonal days display both a higher than usual number of banded grid boxes (consistent with section 6) and a more pronounced clustering of the bands. This is shown in Figs. 13a and 13b, respectively. Figure 13a shows the distributions of the maximum number of banded grid points per 120° longitude sector (chosen to resemble the scale of the largest hetonic explosions) for extreme (white) and nonextreme (gray) days, excluding days/latitudes with no local extremes. To obtain this, the number of banded grid points is computed for all 120° longitude sectors at all latitudes [so the first sector will span 0°–120°E, the second *r*°E–(120 + *r*)°E, where *r* is the dataset’s longitudinal resolution, etc.]. The maximum value at each latitude and day is then selected. Reasonable variations in the longitudinal extent of the sectors do not qualitatively alter the results, and extreme zonal days always show systematically higher values than their nonextreme counterparts. Figure 13b shows the same quantities in terms of anomalies relative to the latitudinal mean. That is, the mean number of banded grid points per 120° sector at a given latitude and time step is subtracted from the maximum before computing the distributions. We take this quantity as a measure of clustering. Again, extreme days display significantly higher values. Indeed, for both panels the two distributions are statistically different under a two-sample Kolmogorov–Smirnov test, with the null hypothesis of the same parent distribution rejected at the 5% significance level. Since the extent of the individual bands is almost unchanged between extreme and nonextreme zonal days (see Table 4), the above points to an enhanced clustering of heat transport bands during extreme zonal days and can be linked to the existence of storm tracks, as discussed in section 7b below.

The picture afforded by Figs. 12 and 13 is therefore consistent with the idealized process by which low-wavenumber hetonic explosions, associated with preferred locations for vortex motion, drive bursts of zonally integrated meridional heat transport. The additional clustering of bands seen on extreme zonal days in ERA-Interim is analogous to larger clusters of individual vortices in the heton model leading to stronger hetonic explosions and thus larger zonally integrated meridional heat transport values. In fact, the hetonic heat transport has a term that depends on the number of vortices leaving the selected domain at a given time step [Eq. (B4b)]. Moreover, the radial velocity with which a vortex cluster breaks away from the rim current—and hence the rate of change of the number of vortices within a given domain—is a monotonically decreasing function of its wavenumber (Hogg and Stommel 1985b). Therefore large clusters will be associated with both a larger number of vortices and a faster motion of such vortices. Naturally, we need to stress that the heton setup is not meant to provide a formal idealization of the real atmosphere, but rather a conceptual model for the process by which a fluid confined to a range of latitudes on a rotating Earth will give rise to instabilities that decrease the system’s potential energy and effect a meridional heat transport.

The aquaplanet GCM, with its more realistic setup, provides a different perspective on the same phenomenon, and a similar analogy to the one discussed above can be drawn with the heton model. Indeed, the heat transport bands in the GCM also display an enhanced tendency to cluster zonally on extreme days relative to nonextreme days (see Figs. 13c,d). As in the reanalysis, the largest zonally integrated transport events in the GCM correspond to large numbers of meridional heat transport bands occurring simultaneously around a given latitude circle. Individual bands typically correspond to the warm or cold sectors of midlatitude cyclones and can be related to the well-known role of baroclinic disturbances in extracting energy from the mean flow by transporting heat poleward (e.g., Blackmon and Lau 1980). Cold sector dynamics account for up to 45% of the heat transport bands at a given latitude, in contrast with the intuitive picture of warm conveyor belts and more generally the circulation associated with the warm sector driving the bulk of the meridional energy transport. This is consistent with Geen et al. (2016) and supports the findings of MC15, who noted that a significant fraction of the local heat transport extremes were associated with negative meridional velocity and moist static energy anomalies. In this regard, the aquaplanet underestimates the contribution of negative *υ*′ and *H*′ anomalies in the high latitudes to the heat transport bands relative to ERA-Interim. This points to land–sea interaction and topography as important factors in determining high-latitude cold-air outbursts.

The wavenumber spectra of the heat transport in both the GCM and reanalysis (Figs. 7c and 11c) highlight the role of a broad range of scales, ranging from wavenumbers 1 and 2 to 5–10. This range mirrors that of the atmospheric motions contributing to the heat transport variability. We interpret the low-wavenumber power as resulting from the constructive interference between transients and the climatological planetary waves (e.g., Garfinkel et al. 2010; Bagget and Lee 2015). Instances of constructive interference have been found to transport heat and moisture to the high latitudes (Goss et al. 2016) and can lead to a more meridionally oriented circulation that diverts synoptic systems poleward (Baggett et al. 2016). These episodes are likely to match our meridional heat transport zonal pulses. The higher wavenumber part of the spectrum instead matches the scale of the individual transport bands. This suggests that the temporal variability and extremes in meridional heat transport by transient eddies are the signature of atmospheric motions spanning a broad range of spatial scales, in agreement with previous analyses (Messori and Czaja 2014). However, we note that there is an important difference in the latitudinal location of the spectral power maxima between the GCM and the reanalysis. The GCM’s circulation is largely dominated by a strong subtropical jet, which likely explains why there is high spectral power in the transport down to 30°–35°N. In ERA-Interim, the more northward peak is likely linked to the midlatitude jet and the contribution of the SH.

### b. Dynamical relevance and implications

While a number of analogies can be drawn between the three datasets, there is a fundamental feature that is unique to the reanalysis: the real world is clearly not zonally symmetric, and the atmospheric motions in ERA-Interim show strong zonal asymmetries, especially in the Northern Hemisphere. Indeed, many of the features seen in the streamfunction snapshots in Fig. 12 are recovered in the climatological mean (not shown). The geographical distribution of the heat transport bands, shown in Fig. 14a, therefore displays preferred zonal locations. The highest frequencies of band occurrence are seen to match the storm-track regions, and the overall distribution resembles closely that of local meridional heat transport extremes (Messori and Czaja 2013). The hatching in Fig. 14 highlights the regions displaying the largest change in band frequency on zonal extreme days versus all other days. Note that, for clarity, the field has been smoothed with a nine-point spatial filter before determining the hatching regions. The zonal extremes are associated with enhanced banding in geographically localized regions, and this accounts for an important fraction of the heightened climatological frequency seen in storm-track regions. This is particularly evident in the North Atlantic, while the correspondence is weaker in the Southern Hemisphere. The extreme zonal days therefore correspond to geographically confined, closely spaced transport bands that, unlike in the idealized models, occur at preferred longitudinal locations.

The collocation of the most active banding regions and the storm tracks is consistent with the role of the warm and cold sectors of midlatitude cyclones discussed in section 7a above. It further supports previous results linking local extremes with positive meridional velocity and MSE anomalies to warm-conveyor-belt-type circulations (MC15). The latter have in turn been associated with the so-called atmospheric rivers—namely, pathways of large water vapor flux in the troposphere (e.g., Zhu and Newell 1998; Bao et al. 2006). A brief analysis shows that there is indeed a close association between ARs and local heat transport extremes. Figure 14b displays the anomaly in the magnitude of vertically integrated horizontal water vapor flux at times and locations associated with a local extreme versus all other days. The anomalies are positive virtually everywhere and, as expected, the highest values are found over the storm-track regions. It is noteworthy that anomalies of around 150 kg m^{−1} s^{−1} are seen in the NH high latitudes and the Arctic basin. We further note that significant differences emerge relative to the map in Fig. 14a, although this is not necessarily surprising since the values in Fig. 14b are unrelated to the frequency of the heat transport extremes. Consistent with the enhanced moisture transport, times and locations associated with local extremes also display a higher rate of occurrence of ARs (i.e., nonzero **Q**_{r} as defined in section 4c above) than all other days, as marked by the hatching in Fig. 14b.

The above analysis further confirms that the variability of the zonally integrated heat transport is driven by a broad range of scales. The individual heat transport bands and local transport extremes are associated with circulation features linked to extratropical cyclones, while the planetary-scale component is associated with varying levels of storm-track activity and the interaction between the different scales (e.g., Baggett et al. 2016). This complements the findings of a recent study (Novak et al. 2015) concerning the wintertime life cycle of the North Atlantic storm track. The authors find that the storm track’s baroclinicity and heat flux are tightly coupled, with high-heat-flux events eroding upstream baroclinicity. This interaction occurs on a range of time scales, leading to fluctuations in both quantities with periods of up to 7–10 days, beyond the typical synoptic scales. A detailed investigation of the spectral features of the zonally integrated transport is beyond the scope of the present paper; we nonetheless verify whether the time scales identified in Novak et al. (2015) are associated with an enhanced transport variability. Figure 15 displays the wavelet spectrum of the wintertime NH zonally integrated transport averaged over 30°–89°N, between 1 and 32 days (see section 3d). The white contour in Fig. 15a marks the cone of influence, which represents the limit beyond which edge effects become important. Darker shades indicate higher spectral power. The black contours mark the 5% significance level. Figure 15b displays the time-averaged spectrum over the full ERA-Interim period considered. The dotted orange line shows a red noise spectrum with the same lag-1 autocorrelation as the transport time series, while the dashed red line shows the 5% significance level. The zonally integrated transport’s spectral power is well above the red noise spectrum at times of up to ~10 days, which is fully consistent with the synoptic and supersynoptic time scales associated with the storm track’s variability. At the same time, we note that the spectral power at these time scales is discontinuous (Fig. 15a), suggesting that the transport’s variability is likely to have additional drivers and is not a simple cyclical process. Indeed, a preliminary analysis of the autocorrelation function of zonally integrated heat transport time series for the Northern Hemisphere winter shows no clear indication of oscillatory behavior (not shown). During summertime, the storm tracks are less active and the meridional heat transport is on average weaker, but there is no reason to expect radically different physical mechanisms to be at play; indeed, the JJA spectrum (not shown) displays heightened spectral power at periods of up to 9 days.

## 8. Conclusions and further thoughts

The present article analyzes the spatial and temporal variability of zonally integrated meridional heat transport in the atmosphere in a hierarchy of three models/datasets: a two-layer quasigeostrophic model, an intermediate complexity aquaplanet GCM, and ERA-Interim. In all three cases the zonally integrated transport is fundamentally sporadic in nature, and its climatological value is sensitive to a small number of extreme events, or pulses—namely, specific latitudes and days at which a very large transport is effected. The aim of this study is to understand the characteristics and underlying drivers of these transport pulses by focusing on large-scale coherent regions of heightened meridional heat transport, which we term bands. MC15 had already linked zonally integrated extremes in reanalysis data to large numbers of closely spaced local extremes. However, analyzing individual local extremes rather than large-scale coherent regions of heightened meridional heat transport made it difficult to provide a physical interpretation for this link and separate the different spatial scales involved. Here, the hierarchy of models we adopt allows for a clearer physical interpretation.

In the heton model, large zonally integrated heat transport values are associated with the formation of large clusters of baroclinic vortices. In the GCM and the reanalysis data the extremes are driven by more frequent, more intense bands, which also display a clustering behavior. The individual bands are associated with the warm and cold sectors of midlatitude weather systems. In the reanalysis, on zonal extreme days bands are more frequent over the major storm-track regions. We further find that part of the variability in the zonally integrated heat transport occurs at time scales consistent with those of the so-called storm-track life cycle (up to 10 days; Novak et al. 2015), superimposed on a background red noise spectrum with a short (~1 day) decorrelation time.

The above considerations lead to several questions for future research. As noted in section 7b, the transport’s high spectral power at supersynoptic time scales is intermittent in time (Fig. 15a). It would be interesting to diagnose whether this is related to the alternation of in-phase versus out-of-phase fluctuations of the different centers of heat flux banding found in both hemispheres (see Fig. 14a) or whether it is the signature of a stochastic component in the variability. Additionally, storm tracks are known to display a significant variability on monthly and longer time scales (e.g., Wettstein and Wallace 2010), which raises the question of whether this might be associated with low-frequency variations in the zonally integrated heat transport. Finally, it would be intriguing to link the clustering of the heat transport bands to traditional metrics for cyclone clustering in the storm tracks.

A parallel avenue for future research would be to relate the spatial and temporal variability of the heat transport we discuss in this study to the energy balance of the high latitudes. A number of studies have highlighted the influence of meridional heat transport and large-scale moisture anomalies on both sea ice and high-latitude temperatures across a broad range of time scales (e.g., Jungclaus and Koenigk 2010; Graversen et al. 2011; Woods et al. 2013; Qian et al. 2015; Graversen and Burtu 2016; Baggett et al. 2016; Woods and Caballero 2016). Because of their large-scale nature and their close link with enhanced moisture transport, it is reasonable to expect that both the local and zonal extremes and the intense heat transport bands could have a significant impact on the polar regions. The enhanced moisture flux in the high northern latitudes seen in Fig. 14b is particularly indicative of such a link. A preliminary analysis indeed shows that zonal extremes occurring north of 60°N lead to above-average temperatures in the Arctic. The mean anomaly in area-averaged 2-m temperature within the Arctic Circle 2 days after the extremes is 1.37 K, which is beyond the 95th percentile of a distribution obtained from random sampling. While beyond the scope of this study, the topic certainly merits a more detailed investigation.

## Acknowledgments

During this research, G. Messori has been funded by the U.K.’s Natural Environment Research Council (RAPID–RAPIT project), Sweden’s Vetenskapsrådet (MILEX project; Grant 2012-40395-98427-17), and the Department of Meteorology of Stockholm University. ERA-Interim data were obtained from the BADC FTP server (at ftp.badc.rl.ac.uk). The aquaplanet and heton model data are available upon request to the authors. The wavelet analysis was performed using code from C. Torrence and G. Compo, available at http://atoc.colorado.edu/research/wavelets/.

### APPENDIX A

#### The Heton Model

The heton model, discussed in section 4 of the main paper, is a two-layer model seeded with point geostrophic vortices. These give rise to deformations of the interface between the two layers, whose nature depends on the sign of the vortices. Any vortex leading to a downward displacement of the interface (namely, a negative PV anomaly in the upper layer or a positive PV one in the lower layer) will be hot. Conversely, any vortex leading to an upward displacement of the interface will be cold. One can further consider combinations of two vortices, each in a different layer. Where both share the same sign, they will form a barotropic pair (Fig. A1a). Where their signs differ, they form a baroclinic pair. The two vortices in a baroclinic pair will both displace the interface in the same direction; one can therefore have hot or cold baroclinic pairs (Fig. A1b). Because they have a well-defined temperature, such pairs will effect a net transport of heat. These heat-transporting baroclinic pairs are termed hetons.

The mathematical framework of the heton model is based on the quasigeostrophic formulation of potential vorticity on an *f* plane. For a full mathematical derivation, we refer the reader to Hogg and Stommel (1985a) and Legg and Marshall (1993). In a cloud of vortices, each vortex will be affected by the presence of all the others. If, in the same model domain, one seeds a number of different pairs, these will therefore interact. At each model time step each vortex will have a velocity (and hence a location) dependent on its interaction with all other vortices. We note that this implies that the model is discretized in time but not in space. Simple setups with very few vortices show that hetons with the same temperature (both hot or both cold) tend to repel if farther away than one Rossby radius of deformation and cluster if closer. On the other hand, hetons with opposite temperatures will attract at far range and split when closer (Hogg and Stommel 1985a). For a uniform array of hetons (or a domain with a reasonably high vortex density, which approximates to a piecewise-constant baroclinic patch), the velocity field at the edges of the array is nonzero, and a rim of baroclinic velocity forms. Individual vortex pairs that attempt to leave the array are split by this rim current and cannot escape. Since, in the near field, hetons with same temperature do not repel, clusters of vortices can form. These effectively behave like a composite vortex and can be strong enough to break through the rim current and propel away from the rest of the hetons. This behavior was first described by Hogg and Stommel (1985b), who called it a hetonic explosion. (Examples of such explosions are shown in Figs. 3 and C1.)

### APPENDIX B

#### Hetonic Heat Transport

Since the heton framework is based on discrete vortices, the heat transport can be directly related to vortex motions. The technique applied here is an extension of Legg and Marshall (1993) and was originally applied in an oceanographic context. It considers hetons associated with a thermal anomaly *e*_{0} seeded at a rate *N*_{r} per unit time within a given domain (here meant to represent an idealized polar cap). The vortices are free to exit and reenter the domain. The thermal anomaly can be explicitly diagnosed as

Here *c*_{p} is the specific heat capacity of dry air, *ρ* is a reference air density at midlevels, *L*_{H} is the combined thickness of the two model levels, and *T*′ is a reference temperature difference between the two layers (see Table C1). The integral is performed over all space. Within the model framework, this thermal anomaly is effectively a measure of vortex strength. The simplest setup is to set the vortex strength to a constant value; however, the model allows for vortices with time-varying strength. In our experiments we designate the initial vortex strength *e*_{0} and the strength of the *i*th vortex at a given time step *e*_{i}. The choice of *e* in our experiments is discussed in appendix C; in order to make our computation widely applicable, we leave it here as a variable.

For *N*^{in} being the number of hetons within the domain, and *N*^{out} being the number without, the following will be true at all times:

Here, is the energy change due to the new hetons being seeded (for *N*_{r} expressed in per unit time). Provided the new seeds are cold hetons, in atmospheric terms this would correspond to a radiative cooling term. The term is the energy change due to changes in vortex strength over time. If vortex strength decays, this term could be interpreted as a damping factor for thermal anomalies in the atmosphere. Since vortex pairs may split, at any given time there may be unequal numbers of top- and bottom-layer vortices in the domain. Therefore, *N*^{in} and *N*^{out} are computed as half the number of individual vortices within and without the domain, respectively. Note that, in the above discussion, the seeding domain and the domain relative to which the heat transport is computed are taken to match. It is, however, entirely possible to adopt the same formulation while computing the transport across a domain boundary that does not coincide with the one of the seeding domain.

For the case where the vortices decay over time at a rate *α*_{d}, the corresponding energy change satisfies

Combining the above with Eq. (B2), the following can be derived:

Here, *H*_{T} is the heat flux across the domain boundary. Terms A and B on the left-hand side of Eq. (B4b) represent, respectively, the local change in heat content (i.e., a storage term) and the diabatic processes within the domain (i.e., a cooling term). The heat transport therefore depends on the number of vortices within and without the domain. Vortex clusters crossing the domain boundary will clearly have a large effect on it. Finally, we note that the heat transport effected by a heton can also be diagnosed directly from the streamfunction as a function of vortex separation relative to the deformation radius (Hogg and Stommel 1985a).

### APPENDIX C

#### Heton Experiment Setups

This section provides a more detailed description of the heton experiment discussed in section 4 of the main paper. The circular domain’s size is selected so as to mimic a polar cap bounded by the 60°N latitude circle, and its radius is set to 4 deformation lengths. There is an initial spinup phase, where the number of vortices grows with time. Since vortices decay over time, equilibrium is reached once the initial seeds decay to zero. The new seeds then balance the decay but do not increase vortex numbers further. This setup mimics a constant radiative cooling and the relaxation of thermal anomalies. To get meaningful outputs, the experiment is run at equilibrium for just under 1 month. This is heavily limited by the computational facilities available, since the simulation requires computing the interaction of each vortex with all the others at every time step.

The above experiment presents a large number of model parameters that need to be set. The model time step of just over 4 h is chosen to be one-tenth of the reference time scale for hetonic interactions, itself dependent on the dimensional length scale *λ* and the reference velocity *υ*_{ref}. The vortex lifetime, from time of seeding to time of disappearance, is set to approximately 1 week to match a realistic thermal relaxation time scale in the atmosphere. Table C1 presents a summary of the model parameters. Two additional parameters that need to be set are the vortex strength *s* and the heton seeding rate *N*_{r}. For a fixed radiative cooling rate *H*, once *s* is set so is *N*_{r}, and vice versa. Here, *s* is simply the nondimensional counterpart of *e* [see Eqs. (B1) and (B2)], the vortex strength in joules, while *N*_{r} is in units of hetons per model time step.

In the idealized setup presented here, it is difficult to determine the most appropriate value for the strength of a single vortex. Combinations of *e* and *N*_{r}, however, can be readily related to a radiative cooling rate. The values set in our experiments correspond to a cooling of just under 250 W m^{−2}, which is slightly higher—but on the same order of magnitude—as the radiative cooling rate of the atmosphere in the mid- to high latitudes (e.g., Trenberth and Stepaniak 2004). This value refers to the net radiative cooling of an air column and, together with the optional heton decay rate, effectively sets the mean heat transport in the model’s steady state.

The approach adopted here was to test a range of combinations of these two parameters. For conciseness, only results from the experiment with the fewest hetons (*N*_{r} = 2, *s* = 50) are presented in the main paper. This is chosen because, having fewer hetons, the heton clusters in the plots are more easily visually identifiable. Other experiments have been performed using *N*_{r} = 4, *s* = 25 and *N*_{r} = 10, *s* = 10. The number of combinations tested was severely limited by the computational cost of running simulation with high seeding rates. It was chosen to keep the same linear relationship between *s* and *N*_{r} for all the experiments. Examples of hetonic explosions for all three experiments are shown in Fig. C1. Figures C1a–c display the results for increasing heton number density. The dots represent the positions of individual vortices at three consecutive model time steps. The exploding vortex clusters are marked by the thick ellipses. The thinner contours mark domains with radii of 4 (seeding domain) and 7 deformation lengths.

As illustrated in the figure, the qualitative behavior of the hetons is very similar across all setups. This is further confirmed by the similar fractional contributions found for the extreme heat transport events across the three experiments. Such conclusions agree with the sensitivity tests of similar models performed in the literature, which have shown that the qualitative model behavior is reasonably insensitive to the strength of the individual hetons (e.g., Legg and Marshall 1993). It should, however, be noted that this paradigm only holds under the assumption that the heton number density is sufficient to drive a rim current around the seeding domain. Indeed, experiments with very low seeding rates over large domains do not display any clustering behavior, and the individual hetons behave as fully independent dynamical constructs.

### APPENDIX D

#### Additional ERA-Interim Figures

Table D1 displays the percentage contributions of *υ*′*H*′ zonal extremes to the poleward-only and net meridional heat transport in the ERA-Interim dataset. Figure D1 shows the PDFs of zonally integrated meridional heat transport for different geographical domains in the same dataset. Figure D2 displays the standard deviation of the composite anomalies in meridional velocity and MSE for banded grid boxes.

## REFERENCES

*Advances in Geophysics*, Vol. 10, Academic Press, 1–82, doi:.

## Footnotes

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

Note that this does not contradict the previous statement regarding the inability of our calculation to provide single-point heat transport values. In fact, this is not the transport at a single point but, rather, the transport across the full domain boundary due to a specific subset of vortices. Also note that the figure of 6.0 PW is obtained by qualitatively assigning hetons to the cluster based on the ellipses shown in Fig. 3a.

^{2}

Since the spatial extent of the clusters is initially computed in terms of model grid boxes, a regular bin spacing in terms of wavenumber would result in some bins corresponding to a larger range of grid boxes than others, thus artificially distorting the shape of the distribution. The bins of the PDF in Fig. 5b have therefore been selected to ensure that an equal number of cluster widths in terms of model grid boxes falls within each bin.