Concurrent heat waves and their linkage to large-scale meridional heat transports through planetary-scale waves

We investigate concurrent heatwaves across the Northern Hemisphere through the linkage between extremes in Meridional Heat Transport (MHT) and in hemispheric land surface temperature (LST). MHT provides fundamental insight on how large-scale thermodynamics links to atmospheric large-scale dynamics, especially through the action of the eddy planetary-scale circulation in the mid-latitudes. The phase and amplitude of these waves can in fact favor the simultaneous occurrence of heatwaves in remote regions, but how this relates to the amount of heat carried by them has not yet been discussed. In this work, we find that the conditional occurrence of extremely weak MHT and extremely warm hemispheric LSTs is significantly more frequent than other conditional occurrences, both in Summer (JJA) and in Winter (DJF). We argue that the combination of extremely weak, in some cases equatorward, MHTs and warm LSTs in JJA are associated with the reversal of the MHT contribution by zonal wavenumber 3, which is in turn associated with the frequency of atmospheric blocking in western Eurasia and the intensity of blockings and jet stream over the Northwestern Pacific. In DJF, the weak, albeit never equatorward, MHT – warm LST events are characterized by a suppression of the climatologically dominant wavenumber 2, which weakens the overall MHT. The flow is anomalously zonal across much of North America and Eurasia, with reduced frequency of atmospheric blockings and downstream displacement of the jet stream, advecting moist and mild air eastward into the continents. Overall, such dynamical pattern correspond to abnormally warm and widespread temperatures in North America, Eastern Europe, and China. The conditional occurrence of extremely weak MHTs and warm hemispheric LSTs is found to be related to between 30 % and over 40 % of extremely warm hemispheric LST days in both seasons.

1 Introduction

A growing body of literature provides evidence of the hemispheric signature of regional heatwaves (cfr. Van Loon and Thompson, 2023). In other words, heatwaves often do not occur in isolation, but rather multiple heatwaves co-occur in geographically remote regions – also referred to as spatially compounding. The interest in these events is motivated by their socio-economic implications: in a highly interconnected world, their impacts can be larger than those of isolated extremes, notably for actors with correlated exposure (Mills, 2005; Tigchelaar et al., 2018; Kornhuber et al., 2020; Vogel et al., 2019; Gaupp et al., 2020). Moreover, co-occurring heatwaves are increasing at a faster pace than all heatwaves (Messori et al., 2024).

The large-scale atmospheric circulation plays a key role in synchronising occurrence of heatwaves across large distances. The literature has specifically highlighted the role of Rossby wave patterns with specific zonal wavenumbers (Petoukhov et al., 2013; Kornhuber et al., 2017, 2020; Strigunova et al., 2022; White et al., 2023). Kornhuber et al. (2020) identifies two specific wave configurations (wave-5 and wave-7) that are responsible for locked, high-amplitude meanders in the jet stream, leading to persistent hot and dry conditions in regions far apart. Wave resonance has been suggested as a triggering mechanism for the formation and amplification of such quasi-stationary wave patterns, especially in boreal Summer (Petoukhov et al., 2013; Kornhuber et al., 2017). The resonant amplification of planetary-scale waves is forced by underlying conditions, such as land-sea contrast, orography and sea-surface temperature anomalies. Quasi-stationary large-scale wave patterns have also been discussed in the context of recurrent Rossby wave packets (Röthlisberger et al., 2019), and have been identified through local wave activity (e.g. Chen et al., 2015; Huang and Nakamura, 2016), or lowpass-filtered mean meridional wind speed at 250 hPa (Röthlisberger et al., 2016a, b; Bartusek et al., 2022). Diagnostics of the waveguidability of the upper tropospheric jet stream (Wirth et al., 2018) have been used to investigate the co-occurrence of surface extremes, such as the heatwaves (Wirth and Polster, 2021; White et al., 2022). Related work has used the notion of quasi-resonant amplification (Petoukhov et al., 2013). The choice of waveguidability definition plays an important role. Ray tracing and the associated “turning latitudes” present multiple limitations when applied to transient waves, and analyses with barotropic flow models showed that the strength of the meridional gradient of potential vorticity better reflects how the jet stream serves as a waveguide for the zonal propagation of such waves (Wirth, 2020; Wirth and Polster, 2021; Segalini et al., 2024).

Several works have emphasized the role of blocking for the occurrence of heatwaves (cfr. Schaller et al., 2018), distinguishing among the different reasons for heatwave development, such as adiabatic, diabatic heating, advection (Röthlisberger and Papritz, 2023). However, the different response of summer blockings and heatwaves to climate change, the latter increasing in frequency (e.g. Christidis et al., 2015; Vogel et al., 2019), while the former remaining mainly unchanged (e.g. Woollings et al., 2018), suggests that the linkage between blockings and heatwaves is not trivial, and that thermodynamic factors should be also taken into account (e.g. Spensberger et al., 2020), especially concerning concurrent heatwaves.

The literature on large-scale heatwaves and their concurrence has therefore primarily focused on large-scale atmospheric dynamical drivers (Lembo et al., 2024). The study of thermodynamics has focussed on local aspects such as subsidence and diabatic heating, overlooking potential large-scale thermodynamic processes. Atmospheric meridional heat transports (MHT) in the mid-latitudes exhibit an intermittent behavior (Swanson and Pierrehumbert, 1997; Messori and Czaja, 2013, 2014, 2015; Messori et al., 2017; Lembo et al., 2019), and sporadic episodes of very strong or weak transport at hemispheric scale are associated with the occurrence of distinct regional and hemispheric weather regimes (Lembo et al., 2022). There is thus a direct link between extremes in the MHT and the atmospheric circulation at multiple spatial scales, which in turn modulates surface temperature extremes. Consistently with this picture, Lembo et al. (2019) and Messori et al. (2017) demonstrated that MHT extremes result from the interplay of waves at both planetary and synoptic scales, with the relative importance of the different scales varying between boreal Summer and Winter.

Here, we hypothesise that concurrent heatwaves correspond to an anomalous meridional temperature distribution on a hemispheric scale, and thus are reflected in the variability of large-scale MHT. This perspective complements that of the large-scale circulation dominating the literature, and for the first time provides a robust dynamic and thermodynamic context for concurrent heatwave episodes. To test the above hypothesis, we analyse concurrent heatwaves in terms of hemispheric land-surface temperature (LST) anomalies, and use the MHT to link their occurrence to the large-scale atmospheric circulation. We focus on both boreal Summer (June–August; JJA) and Winter (December–February; DJF).

The paper is structured as follows: in Sect. 2 data and methods are described. Methodologies for blocking and jet stream detection are outlined in Sect. 2.2.1. Subsequently, the MHT wavenumber decomposition (Sect. 2.2.2), the selection of extremes (Sect. 2.2.3) and significance tests (Sect. 2.2.4) are outlined. Results are presented in Sect. 3. In Sect. 3.1, the different frequency of unconditional and conditional occurrences are described, while in Sect. 3.2 conditional occurrences of extremely weak MHT and extremely warm LST anomalies in boreal Summer are presented. Similar occurrences for boreal Winter are described in Sect. 3.3. An interpretation of what our results imply in terms of linkages between concurrent heatwaves and the general atmospheric circulation is given in Sect. 4, while a summary and key conclusions are drawn in Sect. 5.

2 Data and Methods

2.1 Data

We use ERA5 reanalysis data (Hersbach et al., 2020) over 1979–2022, at a 6-hourly time resolution for the computation of MHT and at daily time resolution for LST anomalies. MHT is computed for the 30–60° latitudinal band (hereafter “mid-latitude channel”), using three-dimensional fields of air temperature (T), specific humidity (q), geopotential height (z), and zonal and meridional velocities (u and v, respectively). The fields originally produced on a Gaussian grid have been retrieved from the Copernicus Climate Change Service with an interpolated longitude per latitude rectangular grid having 0.5×0.5° horizontal resolution and 137 model vertical levels. We also analyse composite anomalies of the 2-dimensional fields: sea-level pressure (SLP), LST, sea-surface temperature (SST), total precipitation, geopotential height at 500 hPa (z₅₀₀) and zonal velocity at 250 hPa (u₂₅₀).

2.2 Methods

2.2.1 Blocking and jet stream diagnostics

For the computation of blocking diagnostics, we refer to the classic Tibaldi-Molteni index (Tibaldi and Molteni, 1990) with slight modifications introduced in Scherrer et al. (2006), as described in Davini et al. (2012). The meridional gradient reversals are computed at every gridpoint as:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{GHGS}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})={\displaystyle \frac{z_{\mathrm{500}}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})-z_{\mathrm{500}}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_S)}{{\mathit{\varphi }}_{\mathrm{0}}-{\mathit{\varphi }}_S}}\text{(2)}& {\displaystyle }\mathrm{GHGN}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})={\displaystyle \frac{z_{\mathrm{500}}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_N)-z_{\mathrm{500}}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})}{{\mathit{\varphi }}_N-{\mathit{\varphi }}_{\mathrm{0}}}},\end{array}$$

where ϕ₀ is set to be ranging between 30 and 75° N, whereas ${\mathit{\varphi }}_{\mathrm{s}}={\mathit{\varphi }}_{\mathrm{0}}-\mathrm{15}$ and ${\mathit{\varphi }}_N={\mathit{\varphi }}_{\mathrm{0}}+\mathrm{15}$. At every latitude λ₀ a blocking event is counted if at any latitude ϕ₀ between 30 and 75° N there is at least a latitude point for which the following conditions are satisfied:

$$\begin{array}{}\text{(3)}& \mathrm{GHGS}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})>\mathrm{0};\mathrm{GHGN}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})<-\mathrm{10}\mathrm{m}(\mathit{°}\mathrm{lat}{)}^{-\mathrm{1}}\end{array}$$

The frequency of blocking as a function of longitude is therefore the ratio between the number of blocking events and the total days within a season.

The blocking intensity follows from Wiedenmann et al. (2002) and is closely related to the aforementioned definition of the blocking frequency. Following the Davini et al. (2012) formulation, once a blocking event is detected at gridpoint (λ₀,ϕ₀) we first define RC(λ₀,ϕ₀) as:

$$\begin{array}{}\text{(4)}& \mathrm{RC}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})={\displaystyle \frac{\frac{z_{\mathrm{500}}({\mathit{\lambda }}_W,{\mathit{\varphi }}_{\mathrm{0}})+z_{\mathrm{500}}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})}{\mathrm{2}}+\frac{z_{\mathrm{500}}({\mathit{\lambda }}_E,{\mathit{\varphi }}_{\mathrm{0}})+z_{\mathrm{500}}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})}{\mathrm{2}}}{\mathrm{2}}},\end{array}$$

where ${\mathit{\lambda }}_W={\mathit{\lambda }}_{\mathrm{0}}-\mathrm{60}$ and ${\mathit{\lambda }}_E={\mathit{\lambda }}_{\mathrm{0}}+\mathrm{60}$, so that RC denotes the average of the geopotential height 60° of longitude upstream and downstream of the blocking. The blocking intensity BI(λ₀,ϕ₀) is then the adimensional value:

$$\begin{array}{}\text{(5)}& \mathrm{BI}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})=\mathrm{100}\left[{\displaystyle \frac{z_{\mathrm{500}}({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})}{RC({\mathit{\lambda }}_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})}}-\mathrm{1}\right]\end{array}$$

It shall be noted that the precondition to have non-zero BI is that a blocking event is actually detected, according to the Tibaldi-Molteni methodology.

For the jet stream intensity, we follow the procedure outlined in Parker et al. (2019). First, we low-pass filter the u₂₅₀ fields over time at every gridpoint with a 10 d Lanczos filter (Duchon, 1979). Then, a zonal average is computed within NH. The jet speed intensity and location at each longitude λ is defined as the maximum velocity detected among all latitude gridpoints and the associated latitude ϕ.

The rationale behind adopting blocking frequency/intensity and jet speed as diagnostics of the hemispheric scale atmospheric circulation follows from the arguments by Dorrington and Strommen (2020), evidencing that the jet stream variability projects on the weather regimes. A proper consideration of this aspect allows to consider additional types of blockings that complete addressing the modes of variability of the large-scale atmospheric circulation accomplished through the traditional weather regimes classification (e.g. Fabiano et al., 2021; Dorrington et al., 2022). Accordingly, we complement the regime analysis of the MHT extremes outlined in Lembo et al. (2022) by introducing blockings and jet maximum speed, relating hemispheric heatwaves to MHT extremes through planetary-scale waves.

2.2.2 Computation of MHT and wavenumber decomposition

The derivation of the MHTs and the zonal wavenumber decomposition follows Graversen and Burtu (2016); Lembo et al. (2019, 2022). We hereby briefly summarize their methodology.

First, at every vertical level, the total energy E is computed as the sum of the kinetic energy, $K=\frac{\mathrm{1}}{\mathrm{2}}\left(u^{\mathrm{2}}+v^{\mathrm{2}}\right)$, and the moist static energy, H:

$$\begin{array}{}\text{(6)}& H=L_{\mathrm{v}}q+c_{p}T+gz\end{array}$$

where q is specific humidity; T is air temperature; z is geopotential height; and the constants L_v, c_p and g are the latent heat of vaporization, the specific heat capacity at constant pressure and the gravity acceleration, respectively.

The meridional transport of E is integrated across a given circle of latitude and over the vertical extent of the atmosphere:

$$\begin{array}{}\text{(7)}& \oint \underset{p_{\mathrm{s}}}{\overset{\mathrm{0}}{\int }}vE{\displaystyle \frac{\mathrm{d}p}{g}}\mathrm{d}x=\oint \underset{p_{\mathrm{s}}}{\overset{\mathrm{0}}{\int }}v(H+K){\displaystyle \frac{\mathrm{d}p}{g}}\mathrm{d}x\end{array}$$

where p_s is the surface pressure. The transport is positive when directed northwards.

Following Graversen and Burtu (2016), at every timestep and for every latitude, the zonal wavenumber decomposition is performed as:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\widehat{\mathcal{F}}}_{\mathrm{0}}(t,\mathit{\varphi })=D\sum _{i=\mathrm{1}}^M{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{a}_{i,\mathrm{0}}^m{a}_{i,\mathrm{0}}^{E}k=\mathrm{0}\text{(9)}& {\displaystyle }{\widehat{\mathcal{F}}}_k(t,\mathit{\varphi })=D\sum _{i=\mathrm{1}}^M{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(a_{i,k}^m{a}_{i,k}^E+b_{i,k}^m{b}_{i,k}^E)k=\mathrm{1},\mathrm{\dots },N\end{array}$$

where the sum is over all M hybrid model levels, i is the hybrid-level index, and a and b are the Fourier coefficients defined as:

$$\begin{array}{}\text{(10)}& {\displaystyle }{a}_{i,k}^{\mathrm{\Psi }}(t,\mathit{\varphi })={\displaystyle \frac{\mathrm{2}}{D}}\int {\mathrm{\Psi }}_i(t,\mathit{\varphi },\mathit{\lambda })\cos \left({\displaystyle \frac{k\mathrm{2}\mathit{\pi }\mathit{\lambda }}{d}}\right)\mathrm{d}\mathit{\lambda }\text{(11)}& {\displaystyle }{b}_{i,k}^{\mathrm{\Psi }}(t,\mathit{\varphi })={\displaystyle \frac{\mathrm{2}}{D}}\int {\mathrm{\Psi }}_i(t,\mathit{\varphi },\mathit{\lambda })\sin \left({\displaystyle \frac{k\mathrm{2}\mathit{\pi }\mathit{\lambda }}{d}}\right)\mathrm{d}\mathit{\lambda }\end{array}$$

Here Ψ is either the mass flux between two hybrid half levels $m=v\frac{\mathrm{d}p}{g}$ or E, $\mathrm{d}p=p_{i+\mathrm{1}/\mathrm{2}}-p_{i-\mathrm{1}/\mathrm{2}}$ is the pressure difference between hybrid half levels, D=2πRcos (ϕ) with R being the Earth’s radius, k being the zonal wavenumber, and t, ϕ, λ the time, latitude and longitude, respectively, and d=360°. N, i.e. the highest considered zonal wavenumber, is taken to be 20, following Lembo et al. (2022). A correction for the k=0 wavenumber accounting for the zonal- and vertical-mean meridional mass-flux variability is applied following Liang et al. (2018): the zonal mean meridional heat transport is obtained by removing from the zonal mean transport for k=0 a term accounting for energy transport associated with the “instantaneous” total meridional mass flux. This term is averaged out in the long-term mean and dependent on a zero-energy reference level. A detailed explanation of this procedure is provided in Lembo et al. (2019).

2.2.3 Selection of extremes

The selection of MHT extremes follows Lembo et al. (2022) with a few additional details explained in Faranda et al. (2024) and references therein. It is based on the Peak over Threshold approach of Extreme Value Theory (EVT,  Pickands, 1975; Coles et al., 2001), after deseasonalization (subtracting the long-term daily mean), and detrending (removing a linear trend) of the data and declustering of extremes. We use an objective declustering method based on the “extremal index”. We first estimate the extremal index, the inverse of the mean cluster size, with the intervals method (Ferro and Segers, 2003). Then we treat threshold exceedances that are separated by less than r below threshold values as belonging to the same cluster. r is the n_cth largest inter-exceedance time, where n_c is the number of clusters estimated based on the extremal index. The MHT extremes are defined on a daily timescales for 1° wide latitude bands within the mid-latitude channel. An EVT-based convergence algorithm is then used to determine suitable percentile thresholds for the selection of extremes in each latitudinal band. This algorithm displays the GPD (Generalised Pareto Distribution) shape parameter as a function of increasing threshold values. The optimal threshold value is then the one at which the shape parameter is already stable against further increasing threshold values and has the lowest uncertainty (Coles et al., 2001; Galfi et al., 2017; Lembo et al., 2022; Faranda et al., 2024). For the declustering and estimation of the GPD shape parameter, we use the “ismev” (Heffernan, 2025) and “extRemes” (Gilleland and Katz, 2016) R packages.

Table 1 Percentage of days for which at least one latitude within each 10° wide latitudinal range displays an extreme MHT. The different columns show separately positive and negative MHT extremes in DJF and JJA. The bottom row shows the percentage of days featuring an extreme anywhere in the mid-latitude channel.

Table 1 displays the percentage of days characterized by the occurrence of a total (i.e. including all wavenumbers) MHT extreme in at least one latitudinal band within three 10° wide latitudinal ranges. These percentage values are remarkably homogeneous across latitudes, although there is a somewhat larger variation when comparing extremely strong and weak MHTs. This is possibly related to the different properties of the extreme MHT distributions, whose skewness has a more pronounced latitudinal variations for strong than for weak extremes (cfr. Lembo et al., 2022).

We compute LST extremes by spatially averaging daily fields across the 20–70° N domain. We remove the linear trend from the field at each gridpoint and deseasonalize, considering DJF and JJA separately. Hemispheric warm and cold extremes are then defined as all values beyond the 95th and 5th percentiles, yielding 202 and 198 d for JJA and DJF, respectively.

In each 10°-wide latitudinal band between 30 and 60° N, co-occurring extreme MHT and hemispheric-mean LST are defined as the set of days during which hemispheric LST extremes co-occur with an extreme MHT.

2.2.4 Significance tests

To assess the role of MHT extremes in the occurrence of LST extremes, we compare at each grid point the distribution along time of LST, SST, SLP, z₅₀₀ and u₂₅₀ during co-occurring hemispheric LST and MHT extremes (LST|MHT) with their distribution during LST extremes without concurrent MHT extremes ($\mathrm{LST}|\stackrel{\mathrm{‾}}{\mathrm{MHT}}$). We test whether the difference between the composite mean of LST|MHT occurrences and the composite mean of $\mathrm{LST}|\stackrel{\mathrm{‾}}{\mathrm{MHT}}$ is random or reflects an emerging pattern related to these events by applying a statistical significance test. The same test is also applied to the difference between composite means of the same variables in correspondence of the strongest 5 % hemispheric heatwaves and of their general distribution, i.e. without any conditioning to MHT or selecting days corresponding to extreme hemispheric heatwaves. In all cases, the null hypothesis is that the two means are equal. The testing procedure consists of three steps: first, grid points showing a sizeable deviation from the mean in correspondence of co-occurring LST and MHT extremes are selected; secondly, neighboring grid points are clustered and sufficiently large clusters are retained; the chosen statistical test is applied to the cluster-averaged observable, with a correction to manage false discovery rate due to multiple testing over clusters. The procedure is described in detail in Appendix A. As an output, a mask is provided for each composite map, highlighting regions where the composite mean of anomalies conditioned to MHT is significantly larger or smaller than that of anomalies conditioned to no MHT, or the composite mean of anomalies corresponding to extreme hemispheric heatwaves are significantly different from the mean of the general distribution of anomalies.

3 Results

Table 2 shows counts of LST|MHT events in each 10° latitudinal band. The combination of extremely weak MHTs and warm hemispheric LST extremes (hereafter ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$) is preferred in both JJA and DJF, with 25 events found for DJF, and 27 for JJA. The lower third of the mid-latitude channel features a smaller number of such days than the other latitudinal bands, with all days recorded there being detected also in the 40–50° N latitudinal band.

Table 2 Number of days featuring co-occurring LST|MHT events for all combinations of type and season, within different latitudinal ranges. The percentages of joint occurrences relative to the total number of MHT extreme events are given in brackets. Frequencies exceeding 10 % are highlighted in bold. For the relative frequency of occurrences in the full 30–60° N latitudinal band, numbers are highlighted in bold when 10 % threshold is exceeded in at least one 10-degrees-wide latitudinal band.

3.1 Unconditioned extremely weak MHTs

Given that ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ conditioned events occur more often than other conditional occurrences, as shown in Table 2, we first aim to consider the relevance of these cases with respect to all MHT_↓ occurrences.

Despite emerging among other combinations, the frequency of ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences in the two seasons is relatively marginal with respect to unconditional MHT_↓ and LST₉₅ occurrences. However, the numbers displayed in Table 2 refer to single-day occurrences. While the persistence of MHT extremes is mostly limited to one or a few days (cfr. Lembo et al., 2019), hemispheric LST extremes in many cases span over consecutive days. 16 out of the 27 combined weak MHT-warm LST events in JJA exceed the threshold for extremely warm temperatures for at least 5 consecutive days – a widely used definition for heatwaves (Solomon, 2007; Radinović and Ćurić, 2012; Coughlan de Perez et al., 2023). If the full length of these 16 heatwave events is taken into account, a total of 81 heatwave days would be counted, which are at least partly associated with MHT_↓ extreme occurrences. This represents 40.1 % of all LST₉₅ days. In DJF, this frequency amounts to 50 d, representing 25.3 % of all LST₉₅, with 11 out of 25 ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ events belonging to heatwaves lasting more than 5 d. If heatwaves of any duration were taken into account, 55.0 % (37.8 %) extremely warm days would be associated with heatwaves featuring at least one ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ event day in JJA (DJF). This hints at a different seasonal behavior of hemispheric heatwaves, compared to regional heatwaves (cfr. Tuel and Martius, 2024). We interpret this as the wintertime dynamics being more efficient in mixing the low-level atmosphere at the synoptic timescales whereby temperature extremes are less likely to occur. For the sake of consistency between seasons, we consider daily occurrences of extremes associated with LST₉₅ events of any duration and focus on ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences in boreal Summer and Winter.

In order to highlight the common features of combined ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences, it is also informative to investigate whether MHT_↓ has a signature in large-scale patterns of circulation and temperature, hinting to a relation between the thermodynamics and the dynamics of the system. Therefore, we show SST, LST, SLP and u₂₅₀ composite means of MHT_↓ events. We separately considered MHT_↓ events in the 50–60° N (Fig. 1) and in 40–50° N (Fig. B1) latitudinal bands. We focus mostly on the former most northerly band, as this shows the largest effects on the considered fields. Maps for extreme events in the 30–40° N latitudinal band are not shown since they follow those of the neighboring 40–50° N, coherently with findings from previous works (Lembo et al., 2019, 2022).

Considerable differences emerge between the two seasons, regarding the extent of regions showing statistical significant patterns. Mostly, SST and LST anomalies are non-significant in JJA (Fig. 1a, c), with the exception of the southern coast of Alaska, where significant negative anomalies emerge, and scattered regions across Northeastern America, where positive significant anomalies are found. However, both seasons show a general pattern of warming over land in the southern part of the mid-latitudes and cooling in the northern part (Fig. 1a, e), which is consistent with negative mid-latitude MHT anomalies bringing less heat from south to north; the associated warming in southern mid-latitudes may contribute to extreme warming events there. The SLP composites (Fig. 1f) are mostly non significant; positive anomalies are found across most of the Northern Pacific, possibly hinting at the recurrence of a Pacific-North American positive phase (PNA+) (cfr. Lembo et al., 2022), but these are only significant across small parts of the basin close to the coast of North America.

DJF (Fig. 1b, d, f, h) is characterized by mostly positive LST anomalies spreading across eastern Europe, north China, and central North America. Remarkably, no widespread SST anomalies (at least being significant) appear associated with the MHT_↓. The SLP anomalies (Fig. 1f) highlight a northeastward displacement of the Icelandic low towards Scandinavia and Russia, as well as a weakening of the area of low pressure covering the Northern Pacific. Accordingly , the zonal gradients in the fields of z₅₀₀ anomalies (Fig. 1g) evidence an overall decrease, especially across the Pacific and the western Atlantic. This reflects a suppression of meridional heat exchanges by planetary-scale waves, coherently with Lembo et al. (2019) and somewhat counterintuitively matches with hemispherically spread of warm LST anomalies.

Figure 1 Composite over MHT negative extremes within the 50–60° N band of detrended-deseasonalized LST, SLP, SST and z₅₀₀ anomalies with respect to the 1979–2022 baseline period in JJA (a, c, e, g) and DJF (b, d, f, h). The red/blue contours denote significant anomalies with respect to the full statistics of the field over time (see Sect. 2.2.4). Grey lines in (e), (f), (g), and (h) denote the climatological mean after removal of the zonal mean at each latitude for SLP (e, f) and z₅₀₀ (g, h). Dotted lines denote negative values, solid lines positive values. Isolines for z₅₀₀ are equally spaced every 35 m² s⁻² for positive values, and 45 m² s⁻² for negative values. Isolines for z₅₀₀ are equally spaced every 410 Pa for positive values, and 350 Pa for negative values.

3.2 Extremely Weak MHTs conditioned to extremely warm hemispheric LSTs in boreal Summer

We now shift our attention to the composite maps of ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ events during JJA, averaging across the 34 d satisfying these conditions, as outlined in Table 2. First of all, it is worth noticing that the majority of these occurrences are associated with MHT_↓ events spanning a broad latitudinal range in the mid-latitudes. This justifies our choice to consider all events occurring in at least one of the three 10-degrees-wide latitudinal bands.

The LST anomalies (Fig. 2a) are found to be largest over central Asia, amounting to up to 6 K warmer than the climatology. Other significant regions are Eastern Europe, the Middle East, parts of east Asia and the western coast of North America. Maps of LST anomalies for each ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrence are shown in the Supplement (Fig. S2), highlighting regions that feature extremely warm LST values. Remarkably, regions more often affected by warm LST are central Eurasia, Northwestern and Northeastern America, consistently with Fig. 2a.

SST anomalies (Fig. 2c) are mainly found to be significant in the North Atlantic sector, exhibiting a positive-phased tripole structure (similarly to what is found for unconditioned LST₉₅ events in Fig. 6c and discussed later in this section) that mimics to some extent the positive phase of the renowned leading Empirical Orthogonal Function (EOF) of SST anomalies in the North Atlantic (cfr.  Watanabe et al., 1999; Yu et al., 2024). One could hypothesize that this pattern is reflecting an influence of SST patterns on heatwave occurrences in the European region. In fact, it has been previously assessed that the SST anomalies and tendencies in the North Atlantic have an influence on the Summer atmospheric circulation, and as such influence near-surface temperatures over Europe (Krüger et al., 2023). In addition, it is relatively well understood that the two leading EOF modes of Summer atmospheric circulation in the North Atlantic-European (NEA) sector, i.e. the Summer North Atlantic Oscillation (SNAO) and the Summer Eastern Atlantic (SEA) patterns, are closely related to the SST tripole mode (cfr. Folland et al., 2009; Osborne et al., 2020).

The SLP composites in Fig. 2b suggest that ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ could be associated with a SNAO+ phase, with the expansion of a ridge from the Atlantic ocean into Northern Europe (consistently with unconditioned LST₉₅ events in Fig. 6e). Extratropical SST anomalies in the Pacific are thought to be less related to the atmospheric circulation, unlike tropical Pacific anomalies (Thomson and Vallis, 2018). This is well reflected in Fig. 2c. SLP composites in the Pacific-North American sector closely align to the patterns of anomalies projected onto the leading mode of North American heatwave variability (Yu et al., 2023, cfr. also Fig. 4a), hinting at a recurrent pattern pairing the large-scale atmospheric dynamics to heatwaves occurring in the region. The dipole pattern featuring high pressure systems in the western part of North America, and low pressure system over the eastern part is also a typical occurrence associated with heatwaves, as recently outlined in Chapman et al. (2025). An in-depth analysis of a notable example of these combined ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences (namely, the concurrent heatwave associated with the Northwestern Pacific heatwave in late June 2021) is provided in the Supplementary material (Sect. S2).

Composite mean of U₂₅₀ anomalies (Fig. 2b) hints at an overall increase of the subtropical jet stream maximum speed in the Pacific and Atlantic domains. This is reflected in the probability distribution functions (PDF) of the jet stream intensity over the two basins (Fig. 3a, c) as compared to the climatological distribution, to the distributions within MHT_↓ samples and within LST₉₅ samples. The distribution associated with ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences exhibits a significant tail towards extremely intense jet stream in the Pacific, that is coherent with the climatological distribution and the distribution for MHT_↓ and LST₉₅. However, the median value is set to be significantly larger than the climatology, although smaller than the median value for LST₉₅ events. Similarly in the Atlantic, the distribution is peculiarly asymmetric towards weak jet stream, while having larger median value than the climatology and than the MHT_↓ distribution.

Given the relevance of atmospheric blocking for the development of heatwaves in Summer, we separately consider composites of geopotential fields in the free atmosphere. To this aim, Fig. 4a shows composite maps of detrended z₅₀₀, while Fig. 5a, c shows the blocking frequencies as a function of longitude and the instantaneous measure of blocking intensity (see Wiedenmann et al., 2002; Davini et al., 2012) for the ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences in JJA. z₅₀₀ anomalies evidence a strengthening of the ridges characterizing the azonal climatological fields over western North America, the central Atlantic and western Eurasia. On the other hand the troughs across east Asia and the central Pacific are weakened. The longitudinal section of blocking frequencies highlights that the western Eurasian ridge indeed signals the presence of blockings when hemispheric heatwaves and weak MHT co-occur. To a lesser extent, these combined events also feature increased occurrence of blockings in Eastern Eurasia and in the eastern Atlantic, features that are shared with other unconditioned LST₉₅ events (cfr. Figs. 5a and 6g). The composite mean map of instantaneous BI, although relatively patchy, denotes that the ridge over western North America stands out during episodes of intense blocking, though not associated with a significant increase in frequency.

The comparison of LST₉₅ (Fig. 6a, c, e, g) and conditioned ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ composites highlights some peculiar features of the latter occurrences of hemispheric heatwaves, distinguishing them from other hemispheric heatwaves. The Eurasian warming is to a great extent unchanged in the two cases. This is not unexpected; in fact (cfr. Röthlisberger and Papritz, 2023) diabatic heating dominates heatwave occurrences over continental and dry regions, where the solar warming related to cloud-free conditions and consequently the absence of local cooling due to evaporation play a major role into the heating of continental masses. On the other hand, adiabatic warming and the role of advection are important where storm tracks are climatologically located. This is the case of the Northwestern America, where orography through downsloping winds and the advection of previously heated warm and humid masses from the oceans are crucial, and to a lesser extent of southeastern Europe, where subsidence embedded in the climatologically dominating Summer anticyclones, play a major role. Therefore, dynamics are important in shaping the occurrence and intensity of heatwaves in this region through concurrent events, as it is evident from the differing patterns compared to Fig. 6. These patterns have been found to be somehow typical of concurrent heatwave occurrences, as argued in Chapman et al. (2025). However, the 2021 case (see Supplement, Sect. S2) highlights that such dynamic features combine with local thermodynamic conditions in enhancing the LST warm anomalies, through diabatic heating and subsidence.

We finally complement the analysis of conditioned ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ by looking at the spectrum of MHT anomalies as a function of zonal wavenumber and latitude for all days related to a heatwave with a ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrence (Fig. 7a). The climatological spectrum is superimposed on the spectrum of anomalies to show whether the magnitude of the changes are sufficient to revert the sign of the transport. In the high latitudes, wavenumbers k=2 and k=3 feature strong negative (southward) anomalies (up to −0.25 PW), that indeed revert the MHT from positive to negative, i.e. determining a southward heat transport for these scales. At higher wavenumbers, however, the positive k=5 anomaly north of 50° N enforces the poleward transport by 50 %, adding 0.15 PW to a baseline climatological 0.3 PW. The climatological maxima of the JJA spectrum are confined at the tropical edge of the k=2 meridional section, where negative anomalies are also found, peaking at about 1 PW. We finally notice that the positive anomaly at k=5 can be understood as a stretching of the sloping maxima of MHT anomalies in the k=5–8 range as seen in the climatology, with maxima going from higher to lower latitudes for increasing wavenumber. We interpret these as the signature of the planetary-scale amplification of the Rossby part of the wave spectrum, peaking at k=5–7 in JJA, in response to heatwave occurrence, as found in Kornhuber et al. (2020) and coherently with what found in Lembo et al. (2022).

Figure 2 Composite mean anomaly maps of (a) LST, (b) SLP, (c) SST and (d) u₂₅₀, averaged across all days featuring an extremely weak MHT and an extremely warm hemispheric LST anomaly in JJA. The red/blue contours denote significant anomalies with respect to the statistics of all LST₉₅ daily events. Grey lines in (b), (d) denote the climatological mean after removal of the zonal mean at each latitude for SLP (b, f) and u₂₅₀ (d, h) Dotted lines denote negative values, solid lines positive values. Isolines for SLP are equally spaced every 410 Pa for positive values, every 350 Pa for negative values. For u₂₅₀ isolines are every 10 m s⁻¹ for negative and positive values.

Figure 3 PDFs of maximum jet speed intensity over specified domains in the Northern Hemisphere for (a) the Pacific basin in JJA, (b) the Pacific basin in DJF, (c) the Atlantic ocean in JJA, (d) the Atlantic ocean in DJF. The Pacific region is comprised in the (160–240) longitudinal range, the Atlantic region in the (280–360) longitudinal range. The PDFs are standardized and obtained via kernel smoothing function at 100 equally spaced points in the range of jet speed intensities computed as described in Sect. 2.1. Grey solid lines refer to the population of all events, black solid lines to the ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences, grey dashed and dashed-dotted lines to MHT_↓ and LST₉₅ occurrences, respectively.

Figure 4 Composite mean anomaly maps of z₅₀₀ detrended-deseasonalized anomalies for conditioned to the occurrence of ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ in (a) JJA and (b) DJF. The red/blue contours denote significant anomalies with respect to the statistics of all LST₉₅ daily events. Grey lines denote the climatological mean z₅₀₀ after removal of the zonal mean z₅₀₀ at each latitude. Dotted lines denote negative values, solid lines positive values. Isolines are equally spaced every 35 m² s⁻² for positive values, every 45 m² s⁻² for negative values. Blank areas denote anomalies that are larger than the color range.

Figure 5 Longitudinal sections of blocking frequencies for (a) JJA, (b) DJF. Grey and black solid lines represent longitudinal sections for all events and for ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences, respectively, while Grey dashed and dashed-dotted lines denote MHT_↓ and LST₉₅ occurrences, respectively. Composite mean maps of BI index for days featuring ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences in c) JJA and d) DJF. Non-negative values are only displayed in gridpoints where at least one day of occurrence exhibits a blocking. BI is a non-dimensional index ranging between 0 and 10 (cfr. Wiedenmann et al., 2002; Davini et al., 2012).

Figure 6 Same as Fig. 1 for all days featuring LST₉₅ occurrences.

Figure 7 Spectrum of MHT anomalies (W) averaged in JJA (a) and DJF (b) across all days featuring a ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ event, as a function of wavenumber k and latitude. Solid contours denote the climatological mean spectrum; isolines are labeled in 10¹⁴ W.

3.3 Extremely weak MHTs conditioned to extremely warm hemispheric LSTs in boreal Winter

We now turn our attention to ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ events in the boreal Winter (DJF). Figure 8, analogously to Fig. 2, displays composite mean anomaly maps of LST, SLP, SST and u₂₅₀ for such daily occurrences in DJF. We identify several salient features:

• The LST anomalies span the majority of continental regions, peaking in the high latitudes, i.e. central Siberia and Canada;

• Scattered cold SST anomalies are present over the central Atlantic, while the western tropical Pacific exhibits large positive anomalies;

• SLP anomalies display extended low pressure systems over the continents, coherently with a combined AR/PNA- pattern (cfr. Lembo et al., 2022);

• u₂₅₀ anomalies show a weakening of the jet stream at the eastern edge of North America and possibly a strengthening downstream over North Africa and the Iberian peninsula. Similarly, the jet stream at the eastern edge of Eurasia is found to be weakened, whereas it is stronger off the western coast of North America;

Unlike in JJA, we do not find a clear signature of prominent extratropical SST anomalies related to ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ events, especially in the Atlantic. The warm anomaly in the tropical Pacific hints at a possible role of El Niño through the mediation of a westerly QBO phase (cfr. Kumar et al., 2022), but we did not find a significant correlation between our combined occurrences and these indices (not shown). The prominent feature emerging from the SLP composites is the excessive zonality over the high latitudes of Eurasia and North America, interrupted by a weakening of the Pacific and Atlantic climatological troughs. The situation is better elucidated by the composite maps of z₅₀₀ anomalies (Fig. 4b), showing that the ridges over the Rockies and western Europe are strengthened, but at the same time confined towards the midlatitudes, coherently with the downstream displacement of the jet. The PDFs of the jet stream maximum velocity for the Atlantic and Pacific domains (Fig. 3b, d) coherently describe a pronouncedly larger median jet speed peak than the climatological average in the Pacific, which is also consistent with the same PDF for LST₉₅ occurrences; the signal, though, is less clear for the Atlantic domain. The frequency and intensity of blockings is described in Fig. 5b, d. We notice that almost all longitudes blockings are less prone to develop compared to climatology, and that this is in line with what occurs in MHT_↓ events. On the contrary, a significant increase in blocking frequency over the eastern Pacific is found in LST₉₅ events, amounting to a 25 % increase. Regarding BI, the most notable feature is that only few regions exhibit at least one blocking occurrence; these are the eastern Pacific, parts of the Atlantic and the adjacent North African regions. Note that in the Atlantic domain there is a clear increase in jet intensity, which is consistent with z₅₀₀ in Fig. 4b and with positive LST over Europe and further into Asia in Fig. 8a. This zonal flow brings weather systems eastward into the Eurasia causing warm anomalies and hampering development of blockings consistently with Fig. 5b. The described pattern is also consistent with the maps of individual ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ displayed in Fig. S2 of the Supplement. In most cases, large parts of central Eurasia are found to be affected by extremely warm LSTs, while North America features LST extremes either in the western or eastern half of the continent.

Figure 8 Same as Fig. 2 for all days featuring combined ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences in DJF.

Comparison of composites for LST₉₅ (Fig. 6b, d, f, h) and conditioned ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ events highlights that the large troughs spanning most of Eurasia and central North America are a common feature of hemispheric heatwaves, although the pattern of SLP is somewhat reversed over the Pacific region. As a result, the Western half of North America is less prone to extreme warming in the combined case. The jet stream is in fact in both cases displaced downstream, but in the combined case it peaks southward with respect to the concurrent heatwave occurrences. Possibly, a role in this different behavior is played by the strong negative SST anomalies that are found in Fig. 6d, but an interpretation of this would require additional investigation. The expansion of the Atlantic ridge towards Europe, that is peculiar of ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences, might also be related to the larger warming of western Europe compared to the other concurrent heatwave events.

Finally, the distribution of MHT across the spectrum of different zonal wavenumbers (Fig. 7b) highlights a dominant anomalously negative k=2 contribution in the middle and high latitudes (up to more than −0.5 PW, i.e. a 40 % increase with respect to its climatological value). Unlike in JJA, this is neither sufficient to revert the sign of the transport in that specific wavenumber, nor to determine a total negative MHT. The main outcome of this anomaly is rather a significant reduction in the meridional gradient of MHT at k=2 wavenumber. On the contrary, the meridional gradient is increased for k=3, with positive anomalies at latitudes higher than 45° N. Overall, the combination of k=2 and k=5 negative anomalies leads to an almost vanishing but not negative total MHT (not shown). Note again that the strong MHT reduction is consistent with a general mid-latitude warming and cooling of the high latitudes, as reflected in the z₅₀₀ increase at low and increase at high latitudes (Fig. 4b), and mid-latitude surface warming patterns (Fig. 8a). The role of k=2 waves in shaping the total transport is coherent with what mentioned above regarding the displacement of troughs and ridges, the former ones overstretched along continents at high latitudes, the latter ones confined towards the midlatitudes. The second order effect of the k=3 also has a signature in the circulation, as seen in particular for the z₅₀₀ patterns of anomaly (Fig. 4b), with the strengthening of the western Pacific, European and east Asian ridges.

4 Discussion

Despite the differences between Summer and Winter boreal seasons in the wavenumber-decomposed MHTs and the associated atmospheric circulation (cfr. Lembo et al., 2019, 2022), the scales of the transport related to spatially compounding extremely warm events, namely concurrent heatwaves, are similar (k=3, k=2, and to some extent k=5). These involve a sometimes negative (i.e. equatorward) transport, counteracting an overall positive transport at other scales. However, since the wavenumber-decomposed transports only provide hemispherically-aggregated information, we are not able to identify the regions where the transport is negative or positive for specific wavenumbers. The composite analysis provided in the previous section allows to discern the local thermodynamic factors characterizing concurrent heatwaves from those associated with planetary-scale waves through the action of large-scale atmospheric dynamics.

In particular, the analysis highlights that it is not only the decreasing intensity of the MHT associated with the planetary-scale wave, but also the phase of the k=3 wave in Summer that favors the emergence of hemispheric heatwaves. In fact, the westward displacement of ridges over continents and the northward tilting of the upstream jet at the same time favors unusual warming in the high latitude continental regions and prompts advection of heat southwards towards mid-latitudes on the continents’ eastern edge. As mentioned before, this is not surprising, especially in Summer, given the localized warming of continental regions due to diabatic heating, most notably in dry environments (Black et al., 2004; Bieli et al., 2015). Changes in frequencies of blockings over Eurasia, and more intense blockings over western North America support this argument, but also the previous analysis of dominant weather regimes (Lembo et al., 2022) suggests the recurrence of ridges over the Western North America, and co-occurring Scandinavian and Siberian blockings being recognized as defining features of MHT_↓ occurrences, that under certain thermodynamic conditions favor the development of concurrent heatwaves. A case study analysis for a notable occurrence (the Northwestern Pacific heatwave in June–July 2021) is provided in the Supplement. The study highlights the role of the negative k=3 contribution in propagating the MHT anomaly southwards; we hypothesize that this propagation occurs downgradient along the Rossby wave spectrum towards higher wavenumbers, and we leave for future study to analyse in-depth how this portion of the spectrum evolves in time during the development of a concurrent heatwave.

In Winter, it is mostly the weakening of the MHT by the k=2 that prompts the enhanced zonality over Eurasia. However, this is a common feature of unconditioned MHT_↓ (Fig. 1f, h; also compare dominant NAO+ and AO weather regimes in Lembo et al., 2022) and LST₉₅ events. The climatological k=2 wave is typically associated with the thermal land-sea contrast, therefore it is not surprising that its phase its geographically locked. The intensification of the k=3 contribution is due to the strengthening of the eastern Pacific and Euro-Atlantic ridges, on one hand damping the warming over the western North America, on the other hand amplifying it over most of Europe. In other words, the combined occurrence of extremely weak MHTs and hemispheric LST warm extremes corresponds to a pattern of hemispheric-wide enhanced zonal circulation, associated with a vanishing heat exchange between low and high latitudes. Temperature anomalies suggest that zonal gradients between land and oceans are damped by this zonal circulation, and at the same time that meridional temperature gradients are to some extent amplified between the Poles and mid-latitudes.

In this work, we deliberately chose to avoid temporal aggregations, therefore we only considered concurrent LST extremes at the daily timescale. In fact, we did not aim here to distinguish the temporal scales at which hemispheric extremes occur from the background conditions that favor their propagation, turning localized extremes into “circum-global teleconnections”. However, while leaving for future study a focused investigation on different temporal aggregations of concurrent heatwaves and on how these reflect on combined occurrences with MHT extremes, we provide here a lagged analysis of combined weak MHT extremes and warm hemispheric LST extreme occurrences. Figure 9 shows counts of MHT extremes as a function of the lag from any LST event, separately for boreal Winter and Summer, for strong and weak MHT extremes. According to a hypothesis test based on a Monte Carlo resampling, for weak MHT extremes (Fig. 9a, c) most of the lags exceeding the significance level are located around lag 0. We also notice, though, that for DJF a few significant lags point towards MHT weak extremes lagging LST extremes by about 25 d. This qualitative analysis generally confirms that the coincident selection of weak MHT extremes and warm hemispheric LST extremes is justified. However, it is also interesting to look at strong MHT extremes (Fig. 9b, d), for which coincident occurrences are rarer than lagged occurrences, and with MHT extremes either leading warm LST occurrences in Summer (Fig. 9b) or lagging them in Winter, in both cases with a 12 d lag. In Summer, this finding is consistent with recent works, relating diabatic heating embedded in moisture transport through synoptic-scale baroclinic perturbations to the downstream development of blocking systems (Pfahl et al., 2015; Steinfeld and Pfahl, 2019; Hauser et al., 2024). We argued here (Sect. 3.2) that ${\mathrm{LST}}_{\mathrm{95}}|{\mathrm{MHT}}_{\downarrow }$ occurrences feature high-latitudinal blockings; it is possible that preceding strong MHT extremes are instrumental in advecting the heat favoring the development of such blockings. We leave for future study to investigate how weak and strong MHT extremes interact through ultra-long (k=2, 3) modes of atmospheric dynamics, favoring the development and propagation of heatwaves along the mid-latitudes.

Figure 9 Counts of days with hemispheric warm LST extremes as a function of lags (in days) with respect to extremely weak MHT days in (a) JJA and (c) DJF, and strong MHT in (b) JJA and (d) DJF. Positive lags mean that MHT extremes lead LST extremes, whereas the opposite holds for negative lags. A dashed thick line denotes the median of counts for random reshuffling with a 3000-iterations Monte Carlo method that preserves clustering and seasonality. The confidence interval at a 90 % significance level is provided in gray. Red markings denote lags for which the number of dots significantly exceed the confidence interval for random reshuffling. An asymmetry index is provided in the subtitle, retrieved as $\mathrm{AI}=(N^{+}-N^{-})/(N^{+}-N^{-})$ with N⁺ being the sum of all counts for each positive lag, N⁻ the sum of all counts for each negative lag. A Monte Carlo two-sided p-value is also provided to denote the significant of the observed asymmetry index.

To our best knowledge, it is the first time hemispheric-scale warm anomalies in Winter are associated with thermodynamic features of the planetary-scale general circulation. Previous works have often linked warm spells in boreal Winter to an excessively zonal flow (Slonosky and Yiou, 2002; Holmberg et al., 2023) or wetter-than-normal conditions (Tuel and Martius, 2024), but a circumpolar extension of these arguments was still missing. Assuming that the Arctic amplification in the future will permanently reduce the meridional temperature gradient, we could hypothesize that this will result in a gradual shift of the k=2 role for MHT extremes from a DJF-like situation, in which enhanced zonal circulation will bring mild air towards the continents, towards a JJA-like situation, with advection of diabatically heated air from the oceans towards the western boundaries of the continents mediated by k=3 suppression or reversal, especially over the Pacific. In this respect, the examination of the notable January 2007 case study (outlined in the Supplement, Sect. S3) highlights that a focused examination of the tropical intrusions and atmospheric rivers toward the western sides of the continents will be crucial in order to disentangle the thermodynamic and dynamic contributions to the MHT extremes. On one hand, diabatic flow is known to be crucial for nonlinear Rossby waves amplification (cfr. Grams and Archambault, 2016). On the other hand, moist processes are fundamental for the development of upstream blocking systems, especially at high latitudes (Pfahl et al., 2015), not only in boreal Summer (Hauser et al., 2024). Heat and latent energy convergence between mid-latitudes and the polar regions is predicted to change significantly with global warming (e.g. Steinfeld and Pfahl, 2019), and therefore a mean state change in the meridional temperature gradient with implications for planetary-scale waves is expected. Understanding how this reverberates on transport and temperature extremes including hemispheric heatwaves is a challenging issue, that is left for a future study.

Finally, our results support previous findings on the typicality in the development of concurrent heatwaves, as highlighted for DJF European heatwaves by Galfi and Messori (2023) and for the 2021 NW Pacific heatwave by Lucarini et al. (2023). That is, that record-breaking events need not result from unprecedented physical drivers, and they can be “typical” with respect to other extremes of comparable intensity. We extend here these arguments about typical concurrent heatwaves to the dominant patterns of the hemispheric circulation.

5 Summary and Conclusions

We investigated the relationship between co-occurring heatwaves in the Northern Hemisphere mid-latitudes and meridional heat transport (MHT). We hereby summarize the main results:

• Extremely weak (in some cases equatorward) MHT and extremely warm hemispheric LST are a more frequent occurrence than other combined occurrences, both in boreal Summer and Winter. The inter-seasonal aspect indicates fundamental physical properties: a weak meridional heat transport leaves more energy for heating in the mid latitudes;

• The combined occurrence of extremely weak MHT and extremely warm hemispheric LSTs has a 11.0 % and 11.5 % relative frequency compared to all weak MHT extreme occurrences in boreal Summer and Winter, respectively. These events are associated with 42.3 % and 30.4 % of extremely warm hemispheric LST days, respectively;

• The combined occurrence of extremely weak MHT and extremely warm hemispheric LSTs in boreal Summer is associated with the reversal of MHT contribution by wavenumber 3, that reverts the sign of the total MHT. An increasing frequency of western Eurasian blocking and intensity of Northwestern Pacific blocking, as well as an increasing intensity of the jet stream in the Pacific region is found to be associated with such occurrences.

• The combined occurrence of extremely weak MHT and extremely warm hemispheric LSTs in boreal Winter is associated with the suppression of the MHT contribution by wavenumber 2. The phase of this wave in fact promotes an enhanced meridional geopotential height gradient, poleward shifted and intensified jet stream, anomalous shift of the low pressure activity towards Scandinavia in the European domain, and a downstream displacement of the jet stream elsewhere, facilitating an overall anomalously zonalised flow, increasing eastward advection of warm and moist air masses across the continents.

The pattern associated with the wavenumber 3 MHT reversal in boreal Summer has been recognized as a typical pattern related to simultaneous heatwave occurrence (Chapman et al., 2025), and it was also found to be related to the renowned record-breaking 2021 Northwestern Pacific heatwave. We argue that the zonality of the atmospheric circulation in boreal Winter, related to the extremely weak MHT, favors an unusual warming of continental regions through the advection of moist and mild oceanic air. Coherently, the blocking frequency at mostly all longitudes is decreased, and the jet stream peak is displaced downstream.

Overall, this work has shown that co-occurring heatwaves at hemispheric scale and the general circulation of the atmosphere under certain conditions can be related to each other by means of the strength of the heat exchanges between the Poles and the Tropics, carried out in the mid-latitudes through planetary-scale eddies. Such eddies are in this respect crucial, and we highlighted that wavenumbers 2 and 3 and their modulation are often instrumental in shaping the simultaneous occurrence of circulation patterns that favor the development of heatwaves at the regional scale. These findings suggest that a careful evaluation of anomalies in such planetary-scale waves may improve our capability to understand concurrent heatwaves, both in boreal Summer and Winter.

Appendix A: Statistical tests for the significance of composite mean anomalies

In principle, a statistical test could be performed at every grid point of a two dimensional field, applying a correction for the false discovery rate (Wilks, 2016). However, in cases involving a large number of simultaneous tests (approximately 130 000 in this study), such corrections result in a substantial loss of statistical power, i.e., the ability of the test to reject the null hypothesis when it is false. Since our primary goal is to identify regions of significant change without requiring fine spatial resolution, we aggregate grid boxes into clusters and perform a single statistical test for each cluster. This approach improves statistical power at the cost of spatial detail.

The clustering methodology can be summarized as follows:

1. grid points for which at least 66 % of the time steps within the LST|MHT events sample exceed or fall below the time-averaged observable for the corresponding grid point are selected;

2. a density-based methodology (Ester et al., 1996) is adopted to cluster neighboring grid points;

3. clusters that are within a radius of ϵ grid points from each other are grouped and considered as a single cluster;

4. only clusters exceeding a given size g_min (in grid points) are selected;

5. a spatial average is performed over each selected cluster or cluster group to obtain time series;

Parameters ϵ and g_min for each observable are given in Table A1. The choice of the parameters is motivated by the different underlying dynamics that we aim to highlight. For instance, surface temperatures are more prone to be affected by complex topography, orography and coastlines than V₂₅₀, hence the smaller g_min and larger ϵ to preserve the spatial structure. Mean values of the time series derived from the LST|MHT and $\mathrm{LST}|\stackrel{\mathrm{‾}}{\mathrm{MHT}}$ samples are compared using a two-sided Student’s t test at a α=0.05 significance level, applying Welch’s modification to account for unequal variances in the two samples. The null hypothesis is that the composite mean of the field in the LST|MHT sample is equal to that of $\mathrm{LST}|\stackrel{\mathrm{‾}}{\mathrm{MHT}}$, while the alternative hypothesis posits that the conditional occurrence of extreme MHT and hemispheric LST events significantly alters the value of the field. The case of V₂₅₀ requires a different procedure: to assess changes in the jet stream, we consider possible differences in the sample 95th percentile of V₂₅₀, which cannot be compared using a Student's t test. To address this, we rely on a non-parametric bootstrap test instead. For each cluster, N=5000 random samples of size n (n being equal to the number of time steps in the LST|MHT sample) are drawn from the full time series. For each bootstrap sample, the empirical 75th percentile is computed. Then, the 75th percentile of V₂₅₀ in the LST|MHT sample is compared to the distribution of the N bootstrap replicates of the same quantity under the null hypothesis. If the LST|MHT value falls in either of the $\mathit{\alpha }/\mathrm{2}$ or $\mathrm{1}-\mathit{\alpha }/\mathrm{2}$ tails of the distribution, we reject the null hypothesis, implying that MHT have a significant effect on V₂₅₀ within the cluster.

Table A1 Parameters for the density-based clustering algorithm.

Since the MHT extremes were selected using declustering (cfr.  Lembo et al., 2022, and Sect. 2.2.3), no pre-whitening or other corrections are applied to address serial correlation. However, since the significance test is repeated for each cluster, and the number of clusters can be of the order of a few dozens, the Benjamini-Hochberg (BH) correction (Benjamini and Hochberg, 1995) with a 0.1 significance level is adopted to control the false discovery rate. This correction is appropriate under the assumption of independent or positively correlated tests, a condition generally satisfied by atmospheric data (Wilks, 2016).

Although a method that allows for generally dependent tests is available, known as Benjamini–Yekutieli correction (Benjamini and Yekutieli, 2001), the latter has been shown to drastically reduce the power of the test in the case of spatially extended climatological data, for which the BH version performs best (Cortés et al., 2020).

Appendix B: Composite mean anomalies related to extremely weak MHT in the 40–50° N latitudinal band

Figure B1 displays composite mean maps of extremely weak MHTs located in the 40–50° N band, similarly to Fig. 1 for the 50–60° N band.

In JJA, no significant anomalies emerge, although the spatial anomaly patterns broadly match those shown in Fig. 1a, c, e, g. Particularly, we notice negative LST anomalies over the southern coast of Alaska, and positive SLP anomalies across the northeastern portion of the Pacific ocean.

In DJF, there is a degree of overlap between significant regions highlighted in Fig. 1b, d, f, h and those shown in Fig. B1e–h. Positive LST anomalies span wide parts of Eastern Europe, China and central North America. SLP anomalies differ from Fig. 1f concerning the latitude of the highest SLP in oceanic ridges, with only the Canadian trough appearing to significantly differ from climatological variability, while the Scandinavian lobe of the circulation having less relevance. Besides that, a trough over eastern Asia appears to be significant, while it is not present in Fig. 1f.

Figure B1 Same as in Fig. 1, for extremely weak MHTs located in the 40–50° N latitudinal band.