Revisiting the reanalysis-model discrepancy in Southern Hemisphere winter storm track trends

Revisiting the SH winter storminess trend discrepancy

We revisit the SH winter storminess trends in reanalyses and CMIP6 and AMIP6 multi-model ensembles (see Methods and Table S1), quantified using vertically integrated bandpass-filtered eddy kinetic energy (hereafter EKE, see Methods) as in previous work⁹. To address (I) and (II), we calculate EKE using a larger number of reanalysis data (expanding from 3 to 8) and model ensemble members (expanding from 16 to 26 for CMIP6 and from 13 to 32 for AMIP6, see Methods). We also ensure that EKE is calculated using the same time frequency and spatial grids: daily-mean, 8 pressure levels, and 1.5° by 1.5° horizontal grids (see Methods).

Expanding the number of reanalysis datasets shows the zonal-mean SH storminess trend exhibits large observational uncertainty (Fig. 1a and see Fig. S1 for time series). The reanalysis trend spread is larger than the multi-model ensemble spread, and not all trends are statistically significant. The largest and the smallest trends come from older generation reanalyses (NCEP2 and MERRA1), whereas the newer generation reanalyses are closer together (ERA5 and JRA3Q). The large spread in reanalysis trends is not significantly affected by start and end dates (Fig. S2), and a large spread is also found using a different metric for storminess (e.g., sea-level pressure variance¹¹, Fig S3). Note that storminess climatology does not exhibit similarly large uncertainty, except for NCEP2 which has about 28% weaker climatology than other reanalyses (Fig. S1).

Six out of eight reanalysis trends, including the newest reanalysis trends, fall outside the 10–90% range of the trend distribution of coupled CMIP6 ensemble which predicts SST (Fig. 1a, Table S2). For five reanalysis trends, no more than one CMIP6 ensemble member simulates a trend as high. Thus, while there is large observational uncertainty, CMIP6 ensemble seems to be underestimating the storminess trends. Although weak, the CMIP6 trend is dominated by greenhouse gas forcing according to the Detection and Attribution Model Intercomparison Project simulations (DAMIP³⁸, see Fig. S4).

Four reanalysis trends fall within the 10–90% range of the trend distribution of prescribed-SST AMIP6 ensemble (Fig. 1a, Table S2). Consistently, comparing the distributions of storminess trends between the two ensembles, the trends in the AMIP6 are significantly greater than those in CMIP6 (p value = 0.01 see Methods), indicating that prescribing observed SST trends significantly increases the storminess trends in the models. Similar results were found in previous work³⁷, where annual-mean atmospheric energy transport trends from transient eddies were greater in the AMIP6 than the CMIP6 ensemble.

The spatial pattern of the storminess trends provides insights into understanding the different trends in CMIP6 and AMIP6 ensembles (Fig. 1c–e). The multi-reanalysis-mean storminess trend is significant across most SH regions including high latitudes of the Indian Ocean, South Pacific, and South Atlantic (Fig. 1c and Fig. S5 for individual reanalysis). However, the CMIP6 ensemble-mean storminess trend in the South Pacific (170°E–60°W) is negligible (Fig. 1d). The AMIP6 ensemble-mean trend better captures the reanalysis trend, especially in the South Pacific, where CMIP6 ensemble shows a near-zero ensemble-mean trend (compare Fig. 1d, e). For more quantitative comparison, storminess trend distributions in the South Pacific are examined for the CMIP6 and AMIP6 ensembles (Fig. 1b). Five reanalysis trends fall outside the 10–90% range of CMIP6 ensemble trend distribution in the South Pacific (Fig. 1b). In contrast, seven reanalysis trends fall within the 10–90% range of AMIP6 ensemble trend distribution (Fig. 1b). Comparing the storminess trend distributions between the two ensembles, the trends in AMIP6 are significantly greater than those in the CMIP6 in the South Pacific (p value S6).

To better understand the role of internal variability related to (II) for SH winter storminess trends, we perform similar analyses using the 50-member Community Earth System Model version 2 Large Ensemble (CESM2-LE) simulations (Fig. 2)^39,40. We also use the 10-member prescribed (observed) SST CESM2 simulations, namely Global Ocean Global Atmosphere (GOGA) simulations (Fig. 2). Both CESM2-LE and GOGA simulations are forced with the same radiative forcing as CMIP6 and AMIP6 ensembles from 1979 to 2014. Here we focus on trends from 1979 to 2013 (see Fig. S6 for CMIP6 and AMIP6). The CESM2-LE simulations also do not capture the observed SST trends (compare Fig. 3a, b).

For CESM2, EKE is defined from monthly anomalies due to data availability (see Methods). The zonal-mean storminess trend in reanalyses quantified using this definition of EKE also exhibits large observational uncertainty (Fig. 2a). It is also notable that the reanalysis trends are greater (especially the statistically significant ones) compared to Fig. 1a by roughly 4–5 times. This is related to using monthly anomalies instead of 2.5–6 day bandpass-filtered anomalies, which includes a broader spectrum of kinetic energy. Nevertheless, changing the EKE definition does not impact the result that storminess trends are greater in the AMIP6 than CMIP6 ensemble (Fig. S7).

Five reanalysis trends, including the newest reanalyses (ERA5 and JRA3Q), fall outside the 10–90% range of the 50-member CESM2-LE trend distribution (Fig. 2a), similar to CMIP6 ensemble in the zonal mean. The CESM2-LE simulations also show a smaller ensemble-mean storminess trend in the South Pacific compared to reanalyses (compare Fig. 2c, d). Similarly, five reanalysis trends fall outside the 10–90% of the CESM2-LE trend distribution in the South Pacific (Fig. 2b). Thus, CESM2-LE simulations underestimate reanalysis trends similar to CMIP6 ensemble, and internal variability cannot fully account for the reanalysis-model mismatch. We also find similar results for CESM2-LE and GOGA comparison (Fig. 2) as the CMIP6 and AMIP6 comparison (Fig. 1, Figs. S6, S7), where the storminess trends in prescribed SST simulations are significantly greater than those in coupled simulations, particularly in the South Pacific.

The analysis above which accounts for (I) and (II) reveals large uncertainty in reanalysis storminess trends in SH. However, most of the coupled simulations (CMIP6 and CESM2-LE) simulate trends that are smaller than most of the reanalyses, including the newest ones, thus they could be underestimating the reanalysis storminess trends. The simulations with prescribed observed SST (AMIP6 and GOGA), on the other hand, show statistically significantly greater storminess trends than the coupled simulations in the zonal mean and the South Pacific. This indicates that (III) coupled simulation’s deficiency in simulating observed SST trends impacts the prediction of SH storminess trends. In what follows, we show that the important factor impacting the storminess trend underestimation in the coupled simulations is their misrepresentation of the SST trends, particularly in the Southern Ocean and tropical Pacific.

Impact of SST trend discrepancies on storminess trends

SST trend discrepancies can impact the storminess through different mechanisms. Southern Ocean SST trends reflect changes in surface fluxes and equatorward ocean energy transport⁴¹. Recent work suggested CMIP6 models underestimate annual-mean SH storminess trends because they do not capture the surface energy flux trends connected to Southern Ocean cooling¹². Furthermore, tropical Pacific SST trends likely impact the SH through Rossby wave teleconnections. More specifically, interannually La Nina leads to a stronger South Pacific storminess^42,43,44, thus a La-Nina-like SST trend would be expected to strengthen the South Pacific storminess trend.

In order to test the hypothesis that SST trend discrepancies contribute to the underestimation of SH winter storminess trends in coupled simulations, we utilize the Southern Ocean (SOPACE) and tropical Pacific (PacPACE) CESM2 pacemaker simulations (see Methods). The pacemaker simulations nudge the Southern Ocean and tropical east Pacific SST anomalies to observations (Fig. 3d, e). Thus, comparing pacemakers to the free-running CESM2-LE simulations allows us to quantify how SH storminess trends in the coupled simulations would change if they correctly simulated the observed Southern Ocean cooling and La-Nina-like SST trend. Since the SOPACE and PacPACE simulations are forced with the same radiative forcing as CESM2-LE simulations, their ensemble-mean differences with the CESM2-LE simulations (Δ_SO and Δ_Pac) quantify the impact of capturing SST trends in the Southern Ocean and the tropical Pacific, respectively (see Methods).

Impact of Southern Ocean SST trend discrepancy on storminess trends

When CESM2 simulations are forced with historical radiative forcings and SST anomalies are nudged to observations in the Southern Ocean (Fig. 3d), zonal-mean storminess trends increase compared to the CESM2-LE simulations (Fig. 4a). Most of the members in SOPACE simulations exhibit larger trends than the CESM2-LE median trend (compare green box and black dashed line in Fig. 4a). More quantitatively, the SOPACE zonal-mean storminess trends are significantly greater than those from the CESM2-LE simulations (p value = 0.02). The strengthening is noticeable across all longitudes in SH, especially at higher latitudes (compare Figs. 2d, 4c).

From an energetic perspective, storminess trends are affected by surface energy flux, top-of-the-atmosphere radiation, and stationary circulation trends (see Eq. 1 in Shaw et al.¹²). Among these factors, Shaw et al.¹² hypothesized the underestimation of annual-mean SH storminess in the CMIP6 models was due to an underestimated surface energy flux trends in CMIP6 models in the SH, which is connected to Southern Ocean cooling trends. The surface energy flux trends reflect ocean heat transport divergence and storage trends⁴⁵. Thus, if the surface energy flux implies equatorward ocean heat transport, it also strengthens the SH storminess because it steepens the atmospheric equator-to-pole energy gradient (i.e., energy is moved away from the South Pole to the equator by the ocean). This influence of surface energy flux trends on the storminess reflecting ocean heat transport trends can be quantified through the moist static energy budget^7,46 with the poleward atmospheric energy transport implied from surface energy flux (F_SFC, see Methods).

The poleward atmospheric energy transport (F_SFC) induced from surface energy flux shows positive ensemble-mean trends across SH in SOPACE simulations (green, Fig. 5a), which is consistent with strengthening SH storminess, similar to ERA5 (blue, Fig. 5a). In contrast, the ensemble-mean trends in CESM2-LE simulations are negative across SH indicating storminess weakening (black, Fig. 5a). Most of the ensemble members in the SOPACE simulations and ERA5 show positive F_SFC trends, whereas most of the members in the CESM2-LE have negative trends (Fig. 5b). Thus, the F_SFC trends from the SOPACE simulations are significantly greater than those from CESM2-LE simulations (p value 5b). Consistent results are found from the CMIP6 ensemble where most of the models simulate negative F_SFC trends (black dashed, Fig. 5b).

Thus, the SOPACE simulations simulate surface energy flux trends consistent with reanalysis trends, which the CESM2-LE simulations struggled to simulate. Consequently, storminess trends are larger in SOPACE than CESM2-LE due to surface energy flux trends that better capture reanalysis trends. The causality is less likely to be the other way around since temperature at the surface is being nudged in the SOPACE simulations.

Note that while surface energy flux and F_SFC trends indicate storminess weakening in CESM2-LE and CMIP6, other factors such as top-of-the-atmosphere radiation trends strengthen the storminess (Fig. S8), leading to weak but positive ensemble-mean storminess trends (e.g., Figs. 1, 2).

While nudging SST anomalies in the Southern Ocean significantly strengthens storminess trends in the zonal mean, its impact in the South Pacific is weak. Quantitatively, the SOPACE and CESM2-LE trend distributions in the South Pacific (Fig. 4b) are not significantly different (p value = 0.29). This suggests that other processes are important for storminess trends in the South Pacific.

Impact of tropical Pacific SST trend discrepancy on storminess trends

When CESM2 simulations are forced with historical radiative forcings and SST anomalies are nudged to observations in the tropical east Pacific (Fig. 3e), South Pacific storminess trends increase compared to the CESM2-LE simulations (Fig. 4b, d). Most of the ensemble members in the PacPACE simulations show larger trends than the CESM2-LE simulation median trend (compare red box and black dashed line in 4b). Thus, the South Pacific storminess trends from the PacPACE simulations are significantly greater than those from the CESM2-LE simulations (p value = 0.02). Note that the zonal-mean trend distribution from the PacPACE simulations is not statistically different (p value = 0.35) from the trend distribution from the CESM2-LE simulations (Fig. 4a).

We hypothesized that the La Nina-like SST trend in the tropical Pacific induces a Rossby wave teleconnection trend to the South Pacific, characterized by weaker subtropical jet and strengthened storminess in the South Pacific, similar to what is seen interannually^42,43,44. Specifically, during La Nina winters, the colder tropical Pacific temperature leads to weaker subtropical jet directly via thermal wind balance. Weaker subtropical jet then leads to upper tropospheric eddy momentum flux convergence, which induces adiabatic descent to warm the midlatitudes (see Eq. 7 in Seager et al.⁴²). The poleward eddy heat flux in the midlatitudes increases to alleviate the adiabatic warming, strengthening the storminess in the South Pacific⁴². As such, the subtropical jet weakening can lead to South Pacific storminess strengthening.

The 250-hPa zonal wind and eddy geopotential (deviation from the zonal mean) trends in the reanalyses and PacPACE simulations show a clear La Nina-like teleconnection pattern that is absent in the CESM2-LE simulations (Fig. 6a–c). In particular, all reanalyses and all members in PacPACE simulations have weakening subtropical jet trends across the South Pacific. The subtropical jet trends in PacPACE simulations are significantly different (p value 6d). This confirms that PacPACE simulations exhibit stronger South Pacific storminess by capturing La Nina-like teleconnection trends in reanalysis.

Combined impact of Southern Ocean and tropical Pacific SST trend discrepancies on storminess trends

To investigate the combined impact of both pacemakers (Δ_Pac + Δ_SO) on the coupled simulations, we create a synthetic large ensemble named SUM with 50 members, which is defined as:

Note that we are adding ensemble-mean impacts (Δ_Pac + Δ_SO) to individual ensemble members of CESM2-LE. This synthetic large ensemble is meant to estimate the results for ensemble simulations that nudge the Southern Ocean and tropical Pacific SST simultaneously. It assumes the ensemble-mean impacts of pacemaker simulations are combined with the forced response and internal variability in the CESM2-LE simulations and that the impacts of pacemakers can be added linearly. A similar approach was taken in Kang et al.³³ to create synthetic ensemble simulations using SOPACE simulations.

The SUM ensemble provides valuable insights since it captures the observed SST trend in the ensemble mean (compare Fig. 3b, c). This is due to the remote impacts of the pacemaker simulations on SST trends outside the nudged area (see dashed lines in Fig. 3d, e). The SOPACE simulations affect the SST trend in the Southeast Pacific and around Antarctica³³ while the PacPACE simulations reverse the trend in the tropical Pacific and enhance the warming in the Southwest Pacific (Fig. 3e).

The SUM ensemble shows storminess trends are larger than individual pacemaker simulations both in the zonal mean and South Pacific (Fig. 4a, b) and exhibit significant storminess trends across the SH (Fig. 4e). More importantly, the trends in the SUM ensemble are not statistically different from the prescribed-SST GOGA simulations in both zonal mean (p value = 0.41) and South Pacific (p value = 0.26).