Linking midlatitudes eddy heat flux trends and polar amplification

Eddy heat fluxes play the important role of transferring heat from low to high latitudes, thus affecting midlatitude climate. The recent and projected polar warming, and its effects on the meridional temperature gradients, suggests a possible weakening of eddy heat fluxes. We here examine this question in reanalyses and state-of-the-art global climate models. In the Northern Hemisphere we find that the eddy heat flux has robustly weakened over the last four decades. We further show that this weakening emerged from the internal variability around the year 2000, and we attribute it to increasing greenhouse gases. In contrast, in the Southern Hemisphere we find that the eddy heat flux has robustly strengthened, and we link this strengthening to the recent multi-decadal cooling of Southern-Ocean surface temperatures. The inability of state-of-the-art climate models to simulate such cooling prevents them from capturing the observed Southern Hemisphere strengthening of the eddy heat flux. This discrepancy between models and reanalyses provides a clear example of how model biases in polar regions can affect the midlatitude climate.


INTRODUCTION
One of the largest signals of climate change in recent and future decades is polar amplification: the stronger low-level warming of the poles relative to lower latitudes. 1,2 Such an amplification acts to reduce the meridional temperature gradient, and several studies have suggested that the recent Arctic amplification has affected midlatitudes extreme weather. [3][4][5][6] However, while no one disputes the existence of Arctic amplification in recent decades, there is still much debate regarding its effects on midlatitudes, since such effects are difficult to separate from the large internal climate variability. [7][8][9][10][11][12][13][14][15][16] Here, we investigate a new aspect of the possible connection between high-latitude temperature changes and the midlatitude circulation: the eddy heat flux, which, as shown below, exhibits robust and clear trends over the last several decades.
Eddy heat fluxes have large climatic impacts at midlatitudes. Not only do they play an integral role in transferring heat from low to high latitudes, but also in driving the mean meridional circulation, and in the initial stages of baroclinic eddy life cycles. [17][18][19] To better understand the behavior of midlatitude eddies, previous studies have tried to relate eddy fluxes to the gradient of the mean temperature. 20,21 For example, using arguments from linear baroclinic instability theory, several studies [22][23][24][25] have tried to relate changes in Eady growth rate to changes in the fluxes. The relation between the eddies and the mean gradient has also been studied using simple diffusive closures, 19,[26][27][28] where the poleward eddy fluxes are assumed to be proportional to the mean meridional gradient, for example, v 0 T 0 / ÀT y , where v is meridional velocity, T is temperature, the subscript y denotes meridional derivative, and over-bar and prime denote mean and eddy terms, respectively. Such closures are also commonly used in baroclinic adjustment theories, where the eddy fluxes act to stabilize the baroclinically unstable flow, and keep it marginally supercritical to baroclinic instability. 29,30 Eddy heat fluxes are known to be more sensitive to lower level changes in temperature gradient rather than to upper level changes, 31,32 as long as the baroclinicity (temperature gradients) is concentrated in the lower levels of the atmosphere, and is not controlled by changes in static stability. 33,34 Thus, the recent (and projected) anthropogenic-induced Arctic amplification would imply a decline in Northern Hemisphere (NH) meridional eddy heat flux, as a smaller heat transport is required to maintain a weaker temperature gradient. In the Southern Hemisphere (SH), on the other hand, the recent multi-decadal Southern Ocean cooling would imply a strengthening of eddy heat flux. The aim of this work is thus to examine the recent trends in midlatitude eddy heat flux, their connection to recent changes in high-latitude temperatures, and the role of anthropogenic emissions in those trends.
Note, that we do not aim to elucidate the mechanisms of polar amplification in this manuscript, which might be affected, for example, by moisture fluxes, 23,35-38 but simply to link the recent polar temperatures trends to the trends in midlatitudes eddy heat fluxes. Our results also differ from previous studies who analyzed the recent trends in NH storm tracks. 39,40 Unlike eddy heat fluxes, which constitute the early development of synoptic eddies and are mostly linked to the meridional temperature gradient, storm tracks, and eddy kinetic energy (EKE) constitute the mature stages of synoptic eddies and are sensitive to both the meridional and vertical temperature gradients 33,34 (see Discussion for more details).

Southern Hemisphere
It is instructive to start by considering the hemisphere where polar amplification has, to date, not been observed. Figure 1a shows the SH 1979-2017 annual eddy heat fluxes (v 0 T 0 , calculated from daily data, Methods) trends in 13 models of the Coupled Model Intercomparison Project Phase 5 (CMIP5) and in three different reanalyses (Methods). We here use the absolute value of v 0 T 0 , so that positive (negative) values indicate strengthening (weakening) in both hemispheres. CMIP5 models (blue bars) show large spread in v 0 T 0 trends: half of the models simulate a strengthening over the last four decades, while the other half simulate a weakening. As a result, the very small multi-model mean (purple bar) shows no robust trend. In contrast, the reanalyses (green bars) show a robust strengthening of v 0 T 0 , which is only captured by a few models. This discrepancy between reanalyses and models can be further seen in the time evolution, relative to the 1979-1989 period, of v 0 T 0 (Fig. 1c): the strengthening in reanalyses is not captured by climate models, which in fact show a monotonic weakening until the end of the 21st century. We next address this discrepancy, and elucidate why climate models show such a large spread while reanalyses do not.
Model spread could stem from two sources: internal variability of the climate system (which is not averaged out in the individual model runs) or differences in the models' formulations. To determine whether internal variability explains the spread across the models, we make use of the community earth system model (CESM) large ensemble (LE) (Methods). Figure 1b is similar to Fig.  1a, but shows the 1979-2017 trends in v 0 T 0 for 40 individual LE members (red bars), and the same reanalyses (green). Most LE members (36) simulate a weakening of v 0 T 0 . Only four members show a strengthening, and it is considerably weaker than the strengthening in the reanalyses. This suggests that internal variability is likely not the main reason for the large spread across the CMIP5 models, and their discrepancy with the reanalyses, assuming that the variability in the other models is similar to the one of the CESM (the spread across the LE members is half the spread across the CMIP5 models).
We thus next examine the role of the different model formulations in the modeled trends of v 0 T 0 . As discussed in the introduction, one expects a strong coupling between v 0 T 0 and the meridional near-surface air temperature (SAT) gradient: it is thus tempting to relate the discrepancy in v 0 T 0 trends between reanalyses and models to the models' inability to capture the recent multi-decadal cooling of the Southern-Ocean surface temperature. [41][42][43][44] To determine if this relation exists, we start by showing in Fig. 2a the correlation between trends in v 0 T 0 and trends in the meridional gradient of SAT (Δ y SAT, estimated as the difference between low, 20°S-40°S, and high latitudes, 55°S-75°S) in CMIP5 models (blue), LE members (red), and reanalyses (green). We also use the absolute value of Δ y SAT so that positive (negative) values indicate strengthening (weakening) in both hemispheres. The 39-year trends in v 0 T 0 are highly correlated with trends in Δ y SAT, with r = 0.92 across the CMIP5 models (and r = 0.87 when including the LE members and reanalyses as well). Not only a good correlation exists between these two quantities, but CMIP5 models which show positive (strengthening) trends in v 0 T 0 (open blue dots) also show positive (strengthening) trends in Δ y SAT, whereas CMIP5 models which show negative (weakening) trends in v 0 T 0 (filled blue dots) also show negative (weakening) trends in Δ y SAT (most of this high correlation stems from transient, rather than from stationary, eddies, Supplementary Fig. 4). Similarly, most LE members show negative trends in both v 0 T 0 and in Δ y SAT. And finally, reanalyses show positive trends in v 0 T 0 and correspondingly positive trends in Δ y SAT.
Since we are analyzing 13 CMIP5 models, it is conceivable that the spread across the models would change if one would account for a larger subset of models. However, the large spread across the CMIP5 models in Δ y SAT also appears when accounting for a larger ensemble of 38 CMIP5 models ( Supplementary Fig. 5). Thus, given the good correlation between Δ y SAT and v 0 T 0 , the large spread in v 0 T 0 across the 13 CMIP5 models with available daily data, is expected to persist even when accounting a larger ensemble of models.
Next we ask: does the spread in trends in Δ y SAT stem from high or low latitudes SAT trends? To answer this, we decompose the trends in Δ y SAT into trends in low and high latitudes SAT separately (Fig. 2c). This shows that most of the spread in Δ y SAT indeed stems from high-latitude temperature trends: while all models and reanalyses show a comparable and positive lowlatitude warming trend of~0.01 k yr −1 (y-axis in Fig. 2c), models with high-latitude warming stronger than their low-latitude warming (situated below the 1:1 dot-dashed line) show negative trends in Δ y SAT and in v 0 T 0 (filled blue and red dots), whereas models with high-latitude warming weaker than their low-latitude warming (situated above the 1:1 dot-dashed line) show positive trends in Δ y SAT and in v 0 T 0 (open blue dots). Unlike most models, reanalyses show high-latitude cooling trends over the last decades (green dots), which are consistent with the robust positive trends in Δ y SAT and in v 0 T 0 (Fig. 2a). Although causality arguments cannot be made based on correlation alone, these results clearly suggest that the model spread in v 0 T 0 trends, and their partial inability to capture the v 0 T 0 trends in reanalyses, might indeed stem from biases in simulated surface temperature at high latitudes. To further demonstrate this, we show, in Supplementary Fig. 6, that while high-latitude SAT trends have high correlations with trends in v 0 T 0 (r = −0.88 across the CMIP5 models, and when including all data sets together), low latitudes SAT trends have low correlation (r = −0.34 across the CMIP5 models, and r = −0.51 when including all data sets together).
To further examine the model biases in high-latitude temperature, and to ascertain whether the reanalyses are indeed unaffected by such biases (reanalyses might also have biases in model formulations), we next compute temperature trends from observational data sets that are untainted by any model biases: we show the NOAAGlobalTemp, GISTEMP, and HadCRUT4 trends (gray dashed lines and gray dots in Fig. 2a and c, respectively) (Methods). The high-latitude cooling and the positive trends in Δ y SAT in reanalyses are also present in all three observational data sets. Given the good correlation between Δ y SAT and v 0 T 0 , this agreement between the reanalyses and observations further corroborates our interpretation that the robust observed strengthening of v 0 T 0 is not an artifact of the reanalyses, and that the large spread across the models stems from biases in temperature at high southern latitudes.
Finally, we demonstrate that if one corrects the models' surface temperature biases one obtains the same robust strengthening of v 0 T 0 found in the reanalyses. For this we make use of a 9-member ensemble of CESM atmosphere-only runs (global ocean global atmosphere, LE-GOGA) simulations, which prescribed observed surface temperature (Methods). Figure 2b is similar to Fig. 1b but shows the last available 37-year (1979-2015) trends in v 0 T 0 in both the LE-GOGA and reanalyses. Unlike the ocean-atmosphere coupled LE runs, which show a weakening in v 0 T 0 , when the sea surface temperatures (SSTs) and sea-ice are prescribed from observations all members show a strengthening in v 0 T 0 (purple bars), similar to the strengthening in the reanalyses (green bars). The strengthening of v 0 T 0 across all GOGA members can be further seen in the time evolution of v 0 T 0 , which accompanies the evolution in the reanalyses (Fig. 2d). This confirms our interpretation that high-latitude biases in simulated surface temperature affect the midlatitude climate trends.
Northern Hemisphere In contrast to the robust observed strengthening in the SH, v 0 T 0 in the NH reanalyses shows a robust weakening over the last four decades (1979-2017) (green bars in Fig. 3a). Relative to the 1979-1989 period, by 2017 v 0 T 0 has weakened by~6%, and is projected to weaken by~20% by the end of the 21st century (Fig.  3c). In addition, unlike the large model spread in the SH, in the NH most CMIP5 models agree on the sign of trends, and also simulate a robust weakening between 1979 and 2017 (blue bars in Fig. 3a). As a result the multi-model mean (purple bar) also shows a weakening of 0.01 km s −1 yr −1 . As discussed in the introduction, such a weakening is consistent with the reduction in the meridional temperature gradient and warming of the Arctic.  shows that the weakening of v 0 T 0 is also correlated with the reduction in Δ y SAT (estimated as the absolute value of the difference between high, 65°N-85°N, and low, 20°N-40°N, latitudes) with r = 0.62 across the CMIP5 models (which, similar to the SH, mostly stem from transient, rather than from stationary, eddies, Supplementary Fig. 4), and r = 0.51 when including the LE members and reanalyses as well. As in the SH, the spread in Δ y SAT is mostly due to the warming at high latitudes (Arctic), and not at low latitudes (Fig. 4c). The importance of high latitudes in changing the meridional temperature gradient was also documented using reanalyses 45,46 and models under increased greenhouse gases. 47 The lower correlation between v 0 T 0 and Δ y SAT in the NH than in the SH may be related to the fact that the longitudinal distribution of Arctic warming is different than that of v 0 T 0 . Note that the weakening of v 0 T 0 does not imply that eddies do not play an important role in the recent Arctic amplification, since the latter might be driven by the enhanced poleward eddy moisture flux 23,35-38 ( Supplementary Fig. 7).
Is the NH weakening of v 0 T 0 part of the forced response to anthropogenic emissions, or merely part of the internal variability of the climate system? First, the fact that all CMIP5 models project that such a weakening will continue in coming decades (Fig. 3c) suggests that the recent decline in reanalyses constitutes part of the emerged forced response to anthropogenic emissions. Second, in order to quantitatively answer such a question one has to disentangle the forced response from the internal variability. Thus, we again make use of the CESM LE, where the spread across its members is due to internal variability alone, and the mean of the ensemble is the forced response. Figure 3b shows the NH v 0 T 0 1979-2017 trends in the LE members (red) and reanalyses (green). Similar to the reanalyses and CMIP5 models, all LE members (except one) show a weakening in v 0 T 0 . As a result, the mean of the LE (yellow bar) also shows a weakening, which is approximately half of the weakening in the reanalyses. Assuming that the LE realistically simulates the internal variability of the climate system (as demonstrated below), this would indicate that the recent decline in v 0 T 0 is partially (half) due to anthropogenic emissions: the other half is due to internal variability. As the weakening is projected to continue in coming decades across all LE members (Fig. 3d), one suspects that the recent decline might constitutes the emergence of the forced response to anthropogenic emissions. We next analyze this question.
The "time of emergence" has been used in previous studies to identify when a forced signal appears as distinct from the internal variability (the noise). 48,49 While different studies have used different definitions for the signal and the noise, in all studies the time of emergence is estimated as the time when the signal exceeds a certain threshold (usually one standard deviation) of the internal variability, defined as the noise. To assess whether an anthropogenic signal can be detected in the recent trends of v 0 T 0 , we here use two different approaches for estimating the time of emergence.
In the first, following previous studies, 49 we use the time evolution relative to a reference period (here we choose 1979-1989) of the LE in order to define the signal and the noise. The signal is defined as the time evolution of the LE mean, and the noise as the time evolution of one standard deviation across all members. Using this approach the forced signal emerges out of the internal variability by 2009.
In the second approach we estimate the time of emergence for each realization separately, in both the LE and reanalyses, by comparing the signal to a distribution that lacks the forced response. 48,50 The signal is computed as trends over different lengths in each realization, and the noise as one standard deviation across all trends with corresponding lengths in the CESM preindustrial control run (Methods). Following previous studies, 50 we use this same noise for calculating the time of emergence in both the LE and reanalyses. The trends are first calculated for each member and reanalysis over 10 years (from 1979 to 1988) and then over consecutive lengths of trends (from 1979 to 1989, 1990...) until the signal emerges. Not only the LE shows that the signal of the weakening in v 0 T 0 has already emerged, but also the reanalyses: the green vertical lines in Fig. 4b show that in all reanalyses the weakening in v 0 T 0 has already emerged out of the internal variability during the 1990s. By 2017, v 0 T 0 trends in reanalyses could not be explained using the sole presence of internal variability, attesting they constitute the forced signal ( Supplementary Fig. 8). Calculating both the signal and noise of the reanalyses using their time evolution, 48 rather than using the preindustrial control run, yields the same time of emergence. To verify whether the LE adequately simulates the internal variability of NH v 0 T 0 , and thus whether its mean represents the forced response, we compare the LE trend distribution with the trend distribution from the preindustrial run and from the observational LE 51 (Supplementary Fig. 9). These comparisons show that the ensemble size is sufficiently large to capture the variability, and that it overestimates the observed variability. Thus, the emergence of the forced signal in the reanalyses might even have occurred before the late 90s.
Is the forced signal that has been detected in the recent weakening of NH v 0 T 0 anthropogenic or natural? To answer that we make use of a 20-member ensemble, which is identical to the LE simulations, but forced without the time varying greenhouse gases (LE-fixGHG, Methods). We start by showing the LE-fixGHG mean time evolution, relative to the 1979-1989 period, of v 0 T 0 (black line in Fig. 3d). Fixing the greenhouse gases results in no forced v 0 T 0 weakening in NH in recent and coming decades, attesting that the recent observed weakening in v 0 T 0 can be attributed to the increase in greenhouse gases. Similarly, comparing the observed 1979-2017 trends in NH v 0 T 0 with the trends in the LE and LE-fixGHG shows that without the time varying greenhouse gases one cannot explain the recent weakening in NH v 0 T 0 (Fig. 4d).

DISCUSSION
The recent trends in NH and SH high-latitude temperatures, notably the remarkable warming of the Arctic, and the resulting changes in the meridional temperature gradient, suggest that, in midlatitudes, v 0 T 0 may has also been changing. We analyze four different reanalyses and show that over the last four decades while v 0 T 0 has been weakening in the NH, it has been strengthening in the SH. The weakening in the NH is associated with a decrease in the meridional temperature gradient, which arises from the warming of the Arctic. In contrast, the strengthening in the SH is associated with an increase in the meridional temperature gradient, which arises from the recent multi-decadal high-latitude cooling of SST over the Southern Ocean.
Since most CMIP5 models fail to capture the cooling of Southern Ocean surface temperature, they are unable to capture the strengthening of v 0 T 0 in the SH. Furthermore, the large spread in SH high-latitude warming across the models results in a large spread in v 0 T 0 trends. Prescribing the surface temperature in models, e.g., by preforming atmosphere-only simulations, yields the observed v 0 T 0 trends in models, attesting that the simulated biases in high-latitude surface temperature affect the midlatitudes as well.
In contrast, the recent weakening of v 0 T 0 in the NH is well captured by climate models. The recent weakening is projected to continue in coming decades, suggesting that the observed trends constitute the emergence of the forced response to anthropogenic (and specifically greenhouse gases) emissions. Using the CESM LE we detect an anthropogenic signature in the recent weakening, as the forced anthropogenic signal already emerged out of the internal variability, which itself explains about half of the recent weakening in NH v 0 T 0 . Given the coupling between the recent changes in v 0 T 0 and Arctic SAT, and since the latter have been attributed to anthropogenic emissions, 52,53 our study offers a clear example of how human influences at high latitudes can affect the midlatitudes climate.
Previous studies, which analyzed the recent changes in midlatitudes eddies, have mostly focused on the NH storm tracks. 39,40 Given that v 0 T 0 and EKE represent different stages in the development of eddies, and have different sensitivity to changes in the temperature field, 33,34 the simulated and observed behavior of v 0 T 0 is not seen in recent EKE trends ( Supplementary  Fig. 10); instead, one see a large spread across the CMIP5 models, with the observations falling within the modeled range, but no clear consensus in the models about the sign of the EKE trends, in either hemisphere. Yet, some similarities are also seen between v 0 T 0 and EKE, although these similarities are driven by different physical mechanisms. For example, while we find that prescribing the surface temperature in models corrects the simulated trends in SH midlatitude v 0 T 0 , prescribing the surface temperature was also found to improve the agreement between simulated and observed climatology EKE in the NH. 54,55 And, while we see a hemispheric asymmetry in v 0 T 0 trends over recent decades, which disappears by the end of the 21st century (Figs. 1c, d, and 3c, d), the projected changes in EKE also show hemispherically asymmetric behavior. 24 In summary, if the recent changes in v 0 T 0 (weakening in NH and strengthening in SH) are to continue in coming decades, they will further impact the midlatitudes climate as eddy heat fluxes are a major player in midlatitudes circulation. Moreover, since v 0 T 0 trends are linked to polar temperature trends, it will be important to correct the high-latitude temperature trend biases in the models, specifically in the SH, in order to produce accurate projections of midlatitudes climate.

METHODS
To examine the recent behavior of the annual stationary and transient eddy heat fluxes in reanalyses and models we make use of daily temperature and meridional wind data, and compute v 0 T 0 , where bar and prime denote zonal and monthly averages, and deviation therefrom, respectively. As v 0 T 0 is maximum at midlatitudes and in the lower part of the troposphere we average it from the surface to 700 mb and between 40°and 70°in the NH, and between 40°and 60°in the SH. We limit the averaging in the SH to 60°in order to avoid artificial near-surface values in pressure interpolated fields over the Antarctic continent ( Supplementary  Fig. 1). However, averaging over a wider region in the SH (i.e., between 40°a nd 70°) does not affect the results (Supplementary Fig. 2). In addition, averaging over the entire atmospheric column yields similar results since most of the changes in v 0 T 0 are at the lower atmosphere ( Supplementary  Fig. 3).

CMIP5 models
We also analyze 13 models that participate in the Coupled Model Intercomparison Project Phase 5 63  Although a few models other than those listed have made daily data needed for calculating v 0 T 0 available, those models show large lowlevel biases in v 0 T 0 , and thus we have excluded them from our analysis.

LE of model simulations
In order to disentangle the forced response to anthropogenic emissions from internal variability, and determine whether recent trends have emerged out of the internal variability, we analyze four experiments with the CESM. The first, is an ocean-atmosphere coupled LE that consists 40 simulations (members) with the same historical and RCP8.5 scenarios as for the CMIP5. 64 The sole difference across the members of the LE is a slight perturbation in the initial condition: each member is initialized with a random difference in atmospheric temperature (O10 À14 K), resulting in different transient responses to the same forcing (internal variability). The second is a nine-member ensemble of atmosphere-only runs (LE-GOGA), which are also forced by the historical and RCP8.5 scenarios between 1880 and 2015. In these simulations the SST and sea-ice are prescribed based on the NOAA ERSSTv4 65 and the Hadley Center HadISST 66 data sets, respectively. The third, is an ocean-atmosphere coupled 1800-year preindustrial control simulation (LE-PI); since the radiative forcing is fixed at year 1850, only internal variability is present in that simulation. The fourth experiment is an ocean-atmosphere coupled LE that consists 20 members with the same historical and RCP8.5 scenarios as for the CMIP5 between 1920 and 2080, but without time-evolving greenhouse gases (LE-fixGHG).

Observations
Monthly mean near-SAT from all above reanalyses and CMIP5 models are validated against three different observed surface temperature data sets: the NOAAGlobalTemp, the GISTEMP v3, 67 and the HadCRUT4. 68 These data sets use a combination of satellite and in situ measurements to produce global surface temperature over land and ocean.