Northern Hemisphere Stratosphere‐Troposphere Circulation Change in CMIP6 Models: 2. Mechanisms and Sources of the Spread

We analyze the sources for spread in the response of the Northern Hemisphere wintertime stratospheric polar vortex (SPV) to global warming in Climate Model Intercomparison Project Phase 5 (CMIP5) and Phase 6 (CMIP6) model projections. About half of the intermodel spread in SPV projections by CMIP6 models, but less than a third in CMIP5 models, can be attributed to the intermodel spread in stationary planetary wave driving. In CMIP6, SPV weakening is mostly driven by increased upward wave flux from the troposphere, while SPV strengthening is associated with increased equatorward wave propagation away from the polar stratosphere. We test hypothesized

• About half of the projected stratospheric polar vortex (SPV) uncertainty in Climate Model Intercomparison Project Phase 6 (CMIP6) can be attributed to stationary planetary wave driving • Projected polar vortex weakening and strengthening are linked to increased upward and equatorward wave propagation respectively • A relationship is found between past lower stratospheric wind biases and SPV projections across CMIP6 models

Supporting Information:
Supporting Information may be found in the online version of this article.

Introduction
The large uncertainty in the projected wintertime Arctic stratospheric polar vortex (SPV) response to increased greenhouse gas concentrations has been shown to contribute to uncertainty in regional surface climate (Karpechko et al., 2022;Manzini et al., 2014;Oudar et al., 2020;Simpson et al., 2018;Zappa & Shepherd, 2017).This SPV uncertainty has not been reduced in the simulations by Climate Model Intercomparison Project Phase 6 (CMIP6) models compared to earlier multi-model intercomparisons (Ayarzagüena et al., 2020;Karpechko et al., 2022).While part of this uncertainty is irreducible due to the inherent internal variability of the climate system (Shepherd, 2014), the first part of this study (Karpechko et al., 2022, hereafter Part 1) showed that half of the spread in SPV projections by the end of the 21st century can be attributed to model uncertainty which could, in principle, be reduced by improving models.
Several mechanisms have been proposed to affect projected SPV changes.Shindell et al. (2001) and Lu et al. (2008) suggested that the strengthened meridional temperature gradient between the warming upper troposphere and cooling lower stratosphere and the tendency for the tropopause height to decrease with latitude will lead to an increase in the westerly thermal wind and an SPV strengthening.This mechanism is partly driven by thermodynamical and radiative processes independent of wave driving.Other studies found a projected strengthening of the equator-to-pole Brewer-Dobson circulation (BDC) that leads to a warmer lower polar stratosphere and a weakening SPV (Karpechko & Manzini, 2012;McLandress & Shepherd, 2009;Scaife et al., 2012;Sigmond et al., 2004).The BDC strengthening was linked to an increased upward planetary wave activity flux, in particular, from stationary planetary waves (Karpechko & Manzini, 2017;McLandress & Shepherd, 2009).Karpechko and Manzini (2017) further found that the magnitude of the increased planetary wave flux into the stratosphere depends on the response of the North Pacific low in the troposphere.However, other studies found unchanged or decreased planetary wave activity flux into the stratosphere (Austin et al., 2003;Butchart et al., 2010;Wu et al., 2019) and no link between changes in the North Pacific low and SPV projections (Wu et al., 2019, hereafter W19), suggesting that the response of the planetary waves to global warming is not robust across climate models.
While driving of SPV variability by planetary wave activity flux is well known (Polvani & Waugh, 2004), there is little understanding of how planetary waves will respond to global warming (Wills et al., 2019).Karpechko and Manzini (2017) proposed that a lengthening of the planetary waves due to the strengthening of the subtropical winds, an argument similar to that by Simpson et al. (2016) regarding lengthening of the intermediate-scale waves, could explain the eastward shift of the climatological low in the North Pacific and increase of the wave activity flux into the stratosphere projected by climate models.Other studies attributed the projected eastward shift and intensification of the North Pacific low to an eastward shift and intensification of the convection in the tropical Pacific (Müller & Roeckner, 2008;Wills et al., 2019;Zhou et al., 2014).
Barents-Kara sea ice loss has also been proposed as a factor increasing planetary wave activity flux into the stratosphere, weakening the stratospheric vortex (Kim et al., 2014;Sun et al., 2015) and pushing it toward Eurasia (Zhang et al., 2016).Recently, Kretschmer et al. (2020) linked the different-signed changes in SPV projections across CMIP5 models to the different timing of the projected sea ice loss in the Barents-Kara Seas region.The projected autumn and early winter sea ice loss is non-linear in time due to projections becoming ice-free at some point, with the rate of loss decaying proportionally to the initial amount of sea ice.Consequently, its influence on the atmospheric circulation is expected to be non-linear, causing a stronger SPV weakening as long as Barents-Kara sea ice is decreasing.Such a nonlinearity in the SPV evolution has indeed been reported in climate model projections (Kretschmer et al., 2020;Manzini et al., 2018).
Another mechanism whereby wave-mean flow interactions can affect the future SPV is changes in the wave propagation within the stratosphere (Sigmond & Scinocca, 2010;Walz et al., 2023).Sigmond and Scinocca (2010) proposed that strengthening of the subtropical jet in response to global warming will lead to a weakening of the subtropical barrier for equatorward planetary wave propagation through reducing (or eliminating) the region with negative refractive index.They further indicated that the refractive properties for planetary waves in the basic state climatology in the lower stratospheric region between the subtropical jet and the polar night jet (so called neck region) is crucial, with models characterized by a weak barrier in the present climate being more likely to lose the barrier in a warmer climate, leading to increased equatorward planetary wave propagation away from the SPV, thereby leading to an SPV strengthening.Walz et al. (2023) propose a similar mechanism based on idealized model results, but in addition to the changes in the subtropical jet they propose that an equatorward shift of the polar vortex in the lower stratosphere leads to a more upright waveguide and permits propagation to high latitudes, which favors planetary wave reflection and SPV strengthening.
Overall, the above discussion demonstrates that various mechanisms are potentially involved in the SPV response to global warming and might contribute to intermodel spread of SPV projections.W19 tested several of these mechanisms using CMIP5 simulations.They found that only 27% of the spread in SPV projections could be related to planetary wave driving and that neither the response of the North Pacific low nor the refractive properties of the background wind climatology were significantly related to the SPV projections.Given the importance of the SPV response for regional climate change and for future stratospheric ozone projections, in this study we evaluate the mechanisms involved in the SPV response using both CMIP5 and CMIP6 simulations.In addition, we attempt to explain parts of the projected SPV spread using a causal network analysis (see Section 2).
The study is structured as follows.The methods and the data are presented in Section 2. Section 3.1 discusses the changes in the stratospheric wave driving and its influence on the polar vortex.Section 3.2 evaluates planetary wave response in the troposphere with Section 3.2.1 analyzing the link to the upward propagation to the stratosphere and Section 3.2.2looking at how tropospheric wave changes are affected by the global warming.Section 3.3 analyzes changes in the refractive properties of the upper troposphere and stratosphere and to what degree they are related to the background climatology.Section 3.4 attempts to explain the spread in the SPV projections using a causal network analysis.Section 4 summarizes and discusses the results.This study is built on the results published in Part 1 and the reader is referred to that paper regarding the uncertainty in SPV response and its implications for regional climate change.

Materials and Methods
We analyze climate simulations by 45 CMIP6 models (Table S1 in Supporting Information S1).Where relevant, the analysis is complemented using simulations from 35 CMIP5 models (Table S2 in Supporting Information S1).For the CMIP6 models, 51% are high-top models (model lid above 1 hPa, Table S1 in Supporting Information S1) and 40% are high-top in CMIP5 (Table S2 in Supporting Information S1).Through most of the analysis we consider CMIP5 and CMIP6 ensembles separately.While combining the ensembles would increase the sample size, we find that in some aspects, such as the contribution of the resolved waves to the SPV spread and the role of the basic state for the projections, the two ensembles behave differently, a distinction which we emphasize.Further, the CMIP5 and CMIP6 models are not necessarily independent, and the slightly different forcings are used between CMIP5 and CMIP6 for the historical simulations and projection scenarios.
For consistency with Part 1 and W19, changes in the SPV are diagnosed by analyzing changes in the zonal wind speed at 10 hPa and averaged between 60º-75ºN (U 10 ).Unless otherwise stated, we focus on the change from 1980 to 2009 (past) to 2070-2099 (future) following high emission scenarios SSP5-8.5 (CMIP6) and RCP8.5 (CMIP5) during the boreal winter (DJF) season.The exception is the sea ice concentration over the Barents-Kara seas region which, following Kretschmer et al. (2020), is averaged over the October-December (OND) season.Note that the sea ice concentration is only available for 42 CMIP6 models (Table S1 in Supporting Information S1).
In Section 3.1.1we diagnose changes in the terms of the zonal momentum equation written in the Transformed Eulerian Mean formulation in the log-pressure vertical coordinate system (Andrews et al., 1987): The left hand side in the equation is the tendency of the zonal mean zonal winds (u), where the bar indicates zonal averaging and t denotes time.The first term on the right hand side is the Coriolis torque due to the meridional residual velocity v * with f being the Coriolis parameter.The second term is the divergence of the Eliassen-Palm (EP) flux (F → ), where ρ 0 is the background air density, a is the Earth radius and ϕ is the latitude.The third term represents the action of the nonconservative forces which, in our case, is represented by the parameterized drag due to breaking gravity waves (GWD) including both orographic and non-orographic gravity waves.Note that the contributions of the meridional and vertical wind advections are not considered here because few models have provided them and their contribution is assumed to be small.While Equation 1 is written for the wind tendency as opposed to the wind itself, a close link between wave driving and the zonal wind variability itself has been established on sufficiently long timescales, such as interannual (Newman et al., 2001;Polvani & Waugh, 2004) and centennial scales (Karpechko & Manzini, 2017;W19).
For a subset of 9 CMIP6 models (see Table S1 in Supporting Information S1), v * , F → , and GWD diagnostics are available as a part of DynVarMIP (Gerber & Manzini, 2016), and only this subset is analyzed in Section 3.1.1.The subset includes five high-top models (55%); thus, it has a similar ratio of high top models as the full set of models.All 9 models report GWD due to orographic gravity waves but only seven models report GWD due to non-orographic gravity waves (Table S1 in Supporting Information S1).
In Section 3.1.2we use F → diagnostics for the quasi-stationary waves calculated using monthly mean data which are available for all models.For consistency with W19, F → in Section 3.1.2is calculated in the (linear) pressure vertical coordinate system.The meridional (F ϕ ) and vertical (F p ) components of F → are given by: here, v is the meridional wind, θ is potential temperature, p is air pressure, and prime indicates deviation from the zonal mean.Note that in a stably stratified atmosphere, F p is opposite in sign to the eddy heat flux v′θ′.
In Section 3.1.2we consider the stationary wave F → budget of the polar stratosphere following Sigmond and Scinocca (2010) and W19 by looking separately at the F ϕ and F p components integrated along the boundaries of a polar box.The polar box is defined similarly to that in W19 and extends from 54°N to 90°N in latitude and from 100 hPa to 10 hPa in pressure.The net flux into the box is defined as: here, ϕ 1 = 54°N, ϕ 2 = 90°N, p 1 = 100 hPa, and p 2 = 10 hPa.The integrated fluxes are defined to be positive when directed upward (F 100 and F 10 ) and equatorward (F 54N ) so that F net is positive when F → is converging toward the box.Thus, F 100 and F 10 are the first and second integrals on the right hand side of Equation 3 respectively taken with the opposite sign.Note that F ϕ vanishes at 90°N.Sensitivity of the results to the box definition is tested in Section 3.1.2.
The global warming level (∆T2m) is estimated as a change in annual mean global-mean 2-m temperature.The change in the strength of the subtropical jet stream (∆U 100 ) is defined, following Karpechko and Manzini (2017), using DJF zonal mean zonal winds at 100 hPa averaged over 20°N-40°N.To diagnose the influence of the eastern equatorial Pacific warming on the atmospheric circulation we follow methodology by Oudar et al. (2020)  In Section 3.3 to test the link between the properties of the background flow and changes in SPV we use the squared refractive index n 2 for zonal wavenumber zero and zero phase speed: here, q ϕ is the meridional gradient of the quasi-geostrophic potential vorticity, N is the Brunt-Väisälä frequency and H = 7 km is the height scale.
We define intermodel consistency in projections when 90% of the models agree on the sign of change.For CMIP6 this means that 41 models (out of 45) are required to agree on the sign, while for CMIP5 32 models (out of 35) are required to agree on the sign.The exception is Section 3.1.1,where all 9 models are required to agree on the sign.
In Section 3.4 a causal inference approach (Kretschmer et al., 2020(Kretschmer et al., , 2021;;Pearl, 2013) is used to quantify the hypothesized pathways between SPV changes and the factors discussed in the introduction.The purpose of this approach is not to prove causality, but to quantify the hypothesized causal pathways and to test their consistency with the data (for an introduction of this approach see Kretschmer et al. (2021).For this, the reported factors and pathways affecting SPV change are summarized in the form of a causal network in Figure 1, with links between variables (the boxes) representing assumed causal relationships.We assume that the SPV response (blue box 9) is directly affected by global warming (red box 1), for example, through thermodynamically driven changes in the thermal wind.Moreover, there are different indirect pathways in which global warming affects the SPV.These include mediation through changes in the vertical propagation of the planetary waves from the troposphere (green box 7) and meridional propagation of the waves within the stratosphere (green box 6).Vertical wave propagation to the stratosphere, in turn, is hypothesized to be driven by changes in the tropospheric planetary waves (green box 5), which are affected by strengthening of the subtropical jet stream (green box 2), equatorial Pacific warming (green box 3) and Barents-Kara Seas ice loss (green box 4).These last factors are assumed to be directly caused by the global warming.Strengthening of the subtropical jet stream is also assumed to affect the meridional propagation of the stratospheric waves.Finally, model basic state zonal wind biases in the neck region (purple box 8) is assumed to affect the stratospheric planetary wave response, and thereby the SPV.Here, we restrict the analysis to the factors and links shown in Figure 1, keeping in mind that this network is only an attempt to provide a full picture, as there are potentially other factors and links at play.In particular, while we analyze changes in the gravity waves in Section 3.1.1we do not include them explicitely to the causal network, partly because these are not available for all models.Also changes in the factors such as BDC and the Quasi-Biennial Oscillation interact with planetary wave propagation but are not considered separately.

Results
Motivated by previous literature, we here analyze how the intermodel spread in SPV projections can be explained by processes depicted in Figure 1.We start by looking at the changes in the stratospheric waves in more detail, because their coupling with SPV is direct and well understood (boxes 6 and 7 in Figure 1).As a next step we demonstrate which of the tropospheric factors (boxes 1-5) are linked to the projected SPV changes, assumed indirectly through the stratospheric waves.Then we show how the model basic state (box 8) affects the SPV projections.Finally, we use a causal inference-guided (Kretschmer et al., 2020(Kretschmer et al., , 2021) ) multiple linear regression analysis of the processes in Figure 1 to quantify the overall contributions.

Changes in the Terms of the Momentum Budget
The multi-model mean changes of the three terms on the right hand side of Equation 1 are shown in Figure 2 together with the changes in the zonal mean winds for the 9 CMIP6 models for which these diagnostics are available.The lack of stippling in the SPV region of Figure 2a shows that there is an uncertainty in the sign of the stratospheric wind response across CMIP6 models, consistent with Part 1 and other studies.The changes in the zonal mean winds for the full set of 45 CMIP6 models are shown in Figure S1 in Supporting Information S1 and they are similar to those shown in Figure 2a.
The multi-model mean across 9 models shows a negative wind tendency due to ∇ • F → (which corresponds to EPflux convergence) throughout the mid-and high-latitude stratosphere above 50 hPa, indicating an increased deceleration of SPV (Figure 2b).The increased deceleration is balanced by a positive wind tendency from Coriolis torque (Figure 2d), both in the multi-model mean and across individual models (not shown).There is also an intermodel spread in ∇ • F → with one model (CanESM5) out of 9 projecting a positive change through most of the polar stratosphere, and consistent with this, a strengthened SPV.Such a strengthened SPV is also projected by several other models (see Table S1 in Supporting Information S1) despite a generally negative polar stratosphere ∇ • F → response in these models.This implies that factors other than changes in ∇ • F → affect SPV changes, as discussed below.
The impact of gravity waves on the zonal mean flow is largest in the mid-latitudes where orographic gravity waves are generated in the mountainous regions (Tibet, Rocky Mountains).All 9 models agree on enhanced easterly momentum tendencies from GWD between 30ºN and 50ºN and 30-10 hPa, and enhanced westerly momentum tendencies just below (Figure 2c).This change in GWD is attributable to the strengthening and upward shift of the subtropical jet which increases the altitude where the orographic gravity waves break, thereby allowing more gravity waves to propagate into the stratosphere (Garcia & Randel, 2008;Shepherd & McLandress, 2011;Sigmond & Scinocca, 2010).There is also a positive change in GWD in high latitudes above 10 hPa in the multi-model mean, although it is not consistent across models.
Figure 3 shows how the changes in the momentum budget terms are correlated with ∆U 10 .As expected, ∇ • F → is positively correlated with ∆U 10 : a stronger convergence leads to weaker winds (Figure 3a).According to Equation 1, the zonal wind tendency is affected by wave driving locally; however, the maximum of the correlation does not coincide with the region used to define U 10 .This is likely a consequence of the non-local action of the wave driving on the residual circulation (Haynes et al., 1991).While zonal winds adjust to the changes in the wave driving, so does the residual circulation.As a result, the residual circulation changes caused by the non-local wave driving also affect the winds, reducing the direct link between the winds and the local wave driving.Also, a portion of the zonal wind anomalies during the DJF season may have resulted from wave driving changes in prior months.
W19 found that less than a third of the intermodel spread in the CMIP5 polar vortex projections can be attributed to intermodel spread in ∇ • F → and suggested that accounting for the GWD forcing, not available in their study, might help to explain the spread.Figure 3b shows that the strongest correlation between GWD changes and ∆U 10 Journal of Geophysical Research: Atmospheres is in the subtropics between 10 and 50 hPa.Here, increased gravity wave forcing ( GWD), mostly due to orographic gravity waves (not shown), contributes to the strengthening of the shallow and deep branches of the BD circulation (Abalos et al., 2021;Butchart et al., 2010).The positive correlation between subtropical GWD driving and ∆U 10 can be understood by considering the influence of the parameterized GWD on resolved wave drag.Sigmond and Shepherd (2014) showed that strengthened gravity wave drag in the subtropics in their model ) in the high latitudes, which slows down the SPV.In agreement with this, we find a positive correlation across models (r = 0.6) between changes in subtropical GWD (10-50 hPa and 20º-40ºN) and polar ∇ • F → (10-50 hPa and 50º-85ºN), which likely explains the positive correlation between the subtropical GWD and ∆U 10 .Furthermore, combined forcing by GWD and ∇ • F → (Figure 3c) does not increase the correlation with ∆U 10 compared to ∇ • F → alone, suggesting that either GWD is small, or its influence on ∆U 10 is mediated by, and already included in, ∇ • F → changes.Thus, for this small subset of CMIP6 models we do not find support for the suggestion by W19 that accounting for GWD forcing would help to explain the spread in ∆U 10 not explained by ∇ • F → .Changes in f v * negatively correlate with ∆U 10 (Figure 3d), which is expected from the momentum budget because f v * largely balances ∇ • F → (Figure 2).

EP-Flux Budget for Stationary Waves in the Polar Stratosphere
While the previous section demonstrated the importance of the changes in the resolved waves for the intermodel spread in SPV, the results cannot be used to analyze the spread across all models.In this section we analyze F → calculated using monthly mean data (Equation 2), which are available for all CMIP models.While F → only represents quasi-stationary wave activity, we find a tight relationship (r = 0.93) between changes in ∇ • F → for the quasi-stationary waves and changes in the full ∇ • F → available from the CMIP6 database (Figure S2a in Supporting Information S1), consistent with W19 for CMIP5 data.Also, the pattern on changes in ∇ • F → for the quasistationary waves (Figure S1a in Supporting Information S1) and its correlation with U 10 (Figure S2b in Supporting Information S1) look broadly similar to those for the full ∇ • F → (Figures 2b and 3a respectively).
Therefore, we discuss the results assuming that most of the wave influence on the stratospheric zonal mean circulation changes is represented by the changes in the quasi-stationary F → .
Following W19 the F → budget over the polar box (see Section 2) is calculated separately for models that project SPV weakening (19 models, see Table S1 in Supporting Information S1 for model names) and strengthening (26 models) by the end of the 21st century using Equation 3.There are substantial differences between these two groups in both upward and equatorward fluxes.In the SPV weakening group (Figure 4a) there is an increase in F 100 , which, together with reduced F 54N , contributes to an increased convergence of the wave activity in the polar stratosphere.While part of this anomalous wave activity flux escapes to the upper stratosphere, the remaining part drives SPV weakening, consistent with expectations based on wave-mean flow interaction theory.
On the contrary, in the SPV strengthening group (Figure 4b), the F 100 does not change much while F 54N increases.
Together with increased F 10 , this results in an anomalous divergence of the wave activity away from the polar stratosphere leading to SPV strengthening.
Comparing our results with those from CMIP5 (W19, as well as Figure 5 and Table 1 below) we find that the contribution of F 100 to differences in SPV changes between the SPV weakening and strengthening models is qualitatively consistent, although it is more pronounced in CMIP6.The difference in F 54N between the weakening and strengthening groups in CMIP6 is large (and of opposing sign) and has a magnitude greater than the difference in F 100 .On the contrary, in CMIP5, the difference in F 54N is small.Note that W19 looked at changes in the F → budget for JFM winds as opposed to DJF, but we find essentially the same results as those in Figure 4 for the JFM winds (Figure S3 in Supporting Information S1).
Figure 5 shows intermodel spread in the F → budget changes and its components and their relation to ∆U 10 for both CMIP5 and CMIP6.Results of the statistical analysis are shown in Table 1.For the combined CMIP5 + CMIP6   models, the convergence of the EP-flux calculated over the polar box explains 37% of the intermodel variance.A larger fraction is explained for CMIP6 models (46%) than for CMIP5 (25%).To assess significance of the differences between CMIP5 and CMIP6 we randomly resampled, with replacement, from the CMIP6 models the same number of models as in CMIP5 and recomputed the correlations between the changes in the F → diagnostics and ∆U 10 1,000 times.The results show that the correlation coefficient between DJF ∆U 10 and ∆F net in CMIP5 is at the edge of the 5th-95th percentile range of the CMIP6 bootstrap sample (Table 1) while the correlations between ∆U 10 and ∆F 54N is outside the 5th-95th percentile range suggesting significant differences between the CMIP5 and CMIP6 ensembles.Note that our results for CMIP5 differ slightly from those published in W19 because of small differences in model selection.
Most models with weakened polar vortex (17 out of 19 for CMIP6 and 13 out of 18 for CMIP5; depicted by the blue dots) have increased F 100 but less than a half of the models with strengthened SPV (11 out of 26 for CMIP6 and 8 out of 17 for CMIP5; depicted by the red dots) have decreased F 100 .For CMIP6, the SPV strengthening is more related to increased F 54N (18 out of 26 models) but the dependence is less pronounced for CMIP5 (10 out of 17) consistent with a weak influence of ∆F 54N on SPV changes in CMIP5 (W19).Moving the southern boundary of the polar box to 45°N, to the region where the equatorward F ϕ strengthens in most models, as was done in Sigmond and Scinocca (2010), weakens the correlations between ∆U 10 and ∆F ϕ , and consequently ∆F net in both CMIP6 and CMIP5, suggesting that not all the EP-flux within this enlarged box affects ∆U 10 .
Figure 4 shows that the anomalous F → convergence extends up to 5 hPa, so that the convergence above 10 hPa would also affect SPV via downward control.The location of the upper boundary of the polar box at 10 hPa was chosen due to the availability of data in CMIP5.For CMIP6 we find that raising the upper boundary from 10 hPa to 5 hPa improves the correlation between DJF ∆U 10 and the convergence from 0.68 to 0.74, while raising to 1 hPa improves the correlation to 0.77.In CMIP6, ∆F 10hPa correlates significantly with ∆U 10 , and the correlation remains significant (r ∼ 0.5) also when raising the boundary to 1 hPa.A part of this flux is likely dissipated in the upper atmosphere below the model's lid, contributing to the SPV change.Furthermore, ∆F 10 significantly correlates with ∆F 100 (r ∼ 0.7); therefore, a part of the correlation with ∆U 10 is because they have a common driver.Importantly, the intermodel spread in ∆F 10 in CMIP6 is relatively small (Table 1) and it decreases further as the boundary is raised (not shown), making the relative contribution of ∆F 10 to the SPV change and its spread smaller than that of ∆F 100 and ∆F 54N .
In CMIP5, ∆F 10hPa does not correlate with ∆U 10 , but its intermodel spread is comparable to that of ∆F 100hPa and ∆F 54N .Changes in the sponge layer, used in the low-top models to avoid spurious reflection of the planetary waves from a model's lid, may affect planetary wave propagation.However, we obtain results similar to those in Table 1 when the analysis is restricted to CMIP5 high-top models only (not shown).Thus, changes in the upper stratospheric EP-flux represent a considerable source of SPV uncertainty in CMIP5 but its physical interpretation is unclear.In summary, we conclude that about half of the spread in ∆U 10 across CMIP6 models is linked to the spread in the changes in stationary wave driving, mostly ∆F 100 and ∆F 54N .We will return to this point in Section 3.4.

Changes in Z′ 500 and Their Coupling With the Stratospheric Changes
We now consider factors affecting changes in the tropospheric waves and their link to the changes in the upward wave fluxes.
We diagnose changes in the tropospheric planetary waves by looking at Z′ 500 .In the mid-latitudes the Z′ 500 change (∆Z′ 500 ) in both CMIP5 and CMIP6 projects onto a wavenumber 2 pattern with decreasing geopotential heights over the central and eastern Pacific, and Atlantic and increasing geopotential heights over North America and Eurasia (Figures 6a and 6c).These projected changes suggest an eastward shift of the climatological planetary waves and there is a remarkable similarity in the changes between CMIP5 and CMIP6 ensembles.Potential causes of these changes include lengthening of the planetary waves due to increase in the upper-level zonal winds, and eastward shift of tropical convection; however, full understanding of the mechanisms involved is still lacking (Wills et al., 2019).The spread across models maximizes over the North Pacific (Figures 6b and 6d) where models disagree on the sign of the changes.
To investigate the association between the changes in the tropospheric waves and the stratospheric circulation we regress DJF mean ∆Z′ 500 onto standardized DJF ∆F 100 (Figures 7a and 7c) and ∆U 10 (Figures 7b and 7d).While ∆Z′ 500 affects upward planetary wave propagation (i.e., upward wave coupling), it can also be affected by wave reflection (i.e., downward wave coupling) or other stratospheric impacts on tropospheric waves.Therefore, a fully causal interpretation of the regression coefficients is not possible, and the following discussion relies on additional information available from nudging experiments (e.g., Simpson et al., 2018) and lagged analysis (e.g., Karpechko & Manzini, 2017).
A significant statistical relationship between ∆F 100 and ∆Z′ 500 is seen across a broad range of latitudes and longitudes (Figures 7a and 7c).The patterns resemble a regression pattern found by Karpechko and Manzini (2017) in a lagged analysis of reanalysis and historical simulations when the troposphere leads, suggesting that upward coupling contributes more strongly to the resulting pattern than downward coupling.The link between ∆F 100 and ∆Z′ 500 is similar in both CMIP5 and CMIP6, with an anomaly correlation coefficient between the two regression patterns north of 20ºN of 0.74.In particular, models with a stronger upward ∆F 100 exhibit a negative change of Z′ 500 over the North Pacific (∆Z′ 500NP , indicated by the green box in Figure 7) in both CMIP5 and CMIP6 (see Table 2).This negative ∆Z′ 500 change represents a deepening of the climatological low there, which constructively interferes with climatological wavenumbers 1 and 2 (Garfinkel et al., 2010), amplifying their amplitude and associated wave propagation into the stratosphere; however, it is mostly changes in wavenumber 1 amplitude (∆Z′ 500W1 ) that contribute to ∆F 100 (Table 2).
Across-model variations in ∆Z′ 500NA and ∆Z′ 500NE also significantly contribute to across-model variations in ∆F 100 in both CMIP5 and CMIP6 (Figure 7).∆Z′ 500 in these regions are correlated with those in ∆Z′ 500NP (Figure S4 in Supporting Information S1).For example, the across model correlation between ∆Z′ 500NP and ∆Z′ 500NA is 0.85 in CMIP6 and -0.63 in CMIP5, suggesting that these regions are picking out different parts of the same anomalous planetary wave response.The correlation between ∆Z′ 500NP and ∆Z′ 500NE is also significant although less pronounced ( 0.50 in CMIP6 and -0.36 in CMIP5), with ∆Z ′ 500NE contributing more strongly to changes in wavenumber 2 amplitude (∆Z′ 500W2 ) than ∆Z′ 500NP (Table 2).
In CMIP6, ∆U 10 is correlated to ∆Z′ 500 in an opposite way than to ∆F 100 , as expected (Figure 7b).The two regression patterns are anticorrelated with an anomaly correlation coefficient of 0.78.Models with a weakening SPV are associated with a more negative ∆Z′ 500NP , in agreement with Karpechko and Manzini (2017), and also with a more positive ∆Z′ 500NA and ∆Z′ 500NE .Consistent with the lack of agreement in the sign of ∆U 10 the signs of ∆Z′ 500NP , ∆Z′ 500NA and ∆Z′ 500NE vary across models (Figure 6).
In CMIP5, the relation between ∆Z′ 500 and ∆U 10 appears different (Figure 7d).In particular, there is no strong correlation between ∆Z′ 500NP and ∆U 10 (Table 2), consistent with the results of W19.On the other hand, a relationship between ∆Z′ 500 and ∆U 10 for CMIP5 is found over the North Atlantic.Here, models with strengthening U 10 exhibit decreasing Z′ 500 over the North Atlantic and increasing Z′ 500 over the subtropical Atlantic contributing to a positive North Atlantic Oscillation (NAO) change.A similar relationship, but much weaker, can also be distinguished in CMIP6.Such a pattern was found earlier in stratosphere-nudging experiments (e.g., Hitchcock & Simpson, 2014;Simpson et al., 2018) indicating that a downward influence of SPV changes on the NAO may dominate in this signal over the Atlantic in CMIP5.

Potential Factors Affecting ∆Z′ 500
We here assess the role of different proposed drivers of tropospheric planetary wave changes (Figure 1) in CMIP models.In this section we use linear regression to investigate spatial structures of the relationships.Once more we  notice that a causal interpretation is limited, given two way troposphere-stratosphere coupling as well as other potentially confounding factors.
First, we look at the association between ∆Z′ 500 and global-mean near-surface warming (∆T2m) across CMIP6 models (Figure 8a), which is strongest south of 40°N.Global warming rates affect the drivers of the planetary waves, but it also may affect SPV by strengthening the thermal wind, that is, a thermodynamical mechanism (see Figure 1).There is a good resemblance between this pattern and the multi-model mean ∆Z′ 500 pattern (Figure 6a).In other words, the projected ∆Z′ 500 pattern scales with ∆T2m and larger warming levels lead to stronger ∆Z′ 500 .
The low correlation between ∆T2m and the circulation change indices (Table 3) does not necessarily mean that they are insensitive to warming levels because different effects may suppress each other, as will be further addressed in Section 3.4.
It was proposed (Cattiaux & Cassou, 2013) that the influence of ∆T2m on ∆Z′ 500 is mediated by the warming of the equatorial Pacific which causes Rossby wave generation by strengthened divergent winds in the upper troposphere, and north-eastward propagation of the waves across the Pacific Ocean and North America toward the  Atlantic Ocean.Equatorial Pacific warming is strongly correlated with ∆T2m (r ∼ 0.9), therefore it is not a surprise that the patterns of the correlations between various Niño indices and ∆Z′ 500 resembles that of ∆T2m (not shown).Instead, following Oudar et al. (2020), we assess how the differential temperature anomalies of the eastern Pacific with respect to the western Pacific (∆N3-N4) is linked to ∆Z′ 500 .Figure 8b suggests that while the differential heating correlates with an eastward shift of the Z′ 500 pattern, it does not significantly affect the areas relevant to explaining the uncertainty in the SPV response (Figure 7b), such as Z′ 500NP , nor it is linked to ∆F 100 (Table 3).
Another potential factor affecting ∆Z′ 500 is the strength of the subtropical jet stream response (∆U 100 , see Section 2).∆U 100 is significantly correlated with ∆T2m (r = ∼0.7 in both CMIP5 and 6), however, unlike ∆T2m, it also correlates significantly with ∆Z′ 500NP , ∆Z′ 500W1 (Figures 8c and Table 3) as well as with the upward wave activity flux (Table 3), consistent with Karpechko and Manzini (2017).Models with more positive ∆U 100 show more negative ∆Z′ 500NP .A similar relationship between ∆U 100 and ∆Z′ 500 is found in CMIP5 models (Figure S5 in Supporting Information S1), although the link between ∆U 100 and the upward wave activity flux is small (Table 3).Note that ∆U 100 can also affect equatorward planetary wave propagation in the stratosphere (Sigmond & Scinocca, 2010); thus there is more than one possible causal link between ∆U 100 and ∆U 10 as also indicated in Figure 1.
Finally, the decline of the Barents-Kara Sea ice (-∆SIC BKS ), coinciding with a strongly warming Arctic, is associated with an increase of Z′ 500 in the Eastern Arctic region and a lowering of Z′ 500 in the Bering Strait-Arctic region (Figure 8d).In the Arctic, this differential heating projects onto wavenumber 1 with a maximum near 30ºE 80ºE and a minimum near 110ºW 180ºW.In mid-latitudes, models with greater Barents-Kara Sea ice loss (-∆SIC BKS ) are associated with a more negative ∆Z′ 500NP , and therefore may contribute to increased F 100 (Table 3).Note, that the across model correlation between -∆SIC BKS and the global warming level (∆T2m), is low (r = 0.22), which is the result of most models having become ice-free already before the end of the 21st century, while global warming levels continue to rise.This nonlinearity of -∆SIC BKS in response to global warming would also imply a nonlinearity in its effect on upward planetary wave propagation (Kretschmer et al., 2021).This makes it overall harder to understand the role of sea ice.Accounting for potential nonlinearities and different timing of models becoming ice-free was outside the scope of this work.

Influence of the Basic State on Changes in the Stratospheric Planetary Waves
Strengthening of the subtropical jet affects the refractive properties of the zonal mean stratospheric flow which in turn affects wave propagation (Sigmond & Scinocca, 2010), a process represented by box 8 in Figure 1.It was shown in Section 3.1 that models projecting strengthening of SPV (positive ∆U 10 ) are associated with an increase in the equatorward flux of planetary wave activity in the stratosphere (∆F ϕ ). Figure 9 shows changes in n 2 (Equation 4) separately for models simulating weakening and strengthening of SPV.While both groups of models show a general increase of n 2 in the subtropical lower stratosphere, the key difference, as we argue following Sigmond and Scinocca (2010), is that the strengthening SPV group projects a near disappearance of the region with a negative n 2 (black region in Figure 9).The region with a negative n 2 represents a barrier for equatorward wave propagation, refracting the waves toward higher latitudes instead.Consequently, an erosion of the barrier allows waves to propagate further equatorward contributing to a positive ∆F ϕ and strengthening of the SPV.Note that n 2 in Figure 9 was calculated for wavenumber zero to allow for a like-for-like comparison with previous literature (e.g., Sigmond & Scinocca, 2010;W19), and that the actual n 2 for planetary wavenumber 1, and even more wavenumber 2, would be more negative.Hence, the barrier for higher wavenumbers might still be present, at least in a climatological sense.Sigmond and Scinocca (2010) proposed that future changes in the barrier depend on the basic state refractive index.Models that have a smaller barrier in the basic state climate will more likely lose the barrier in a warmer climate, leading to increased equatorward ∆F ϕ and strengthened SPV.Sigmond and Scinocca (2010) proposed this link in the context of a single model study.W19 subsequently searched for this linkage across the CMIP5 models but found no significant relationship between the basic state n 2 or zonal wind structure in that region and SPV change.In contrast, for CMIP6, we do find a difference in basic state (past) n 2 between the models that exhibit a weakening SPV and the models that exhibit a strengthening SPV that is aligned with the Sigmond and Scinocca (2010) hypothesis, as shown in Figure 9.The models that project a positive ∆U 10 also have a weaker barrier in the past climate, while the models that project a negative ∆U 10 have a stronger barrier in the past climate which remains relatively large also in the warmer climate.
This contrast between the CMIP5 and the CMIP6 model behavior is further demonstrated in Figure 10 which shows the link between the past zonal wind climatology and the change in the SPV.The refractive index is closely linked to the structure of the zonal mean zonal wind, which is a more observable quantity.Sigmond and Scinocca (2010) argued that the models with stronger climatological winds in the "neck" region between the subtropical jet stream and the SPV also have a smaller region of climatological negative basic state n 2 and a weaker barrier to wave propagation.We would, therefore, expect a positive correlation between the basic state climatological winds in that region and ∆U 10 which is, indeed, found in CMIP6 (Figure 10c) while it is entirely absent in CMIP5 with a weak negative correlation with climatological winds found in that region (Figure 10a).
The correlation between ∆U 10 and basic state winds maximizes in the neck region and the correlation between UNECK (see Section 2) and ∆U 10 is 0.57, which is statistically significant at p = 0.05 (Figure 10d) and is also robust to different scenarios since a similar connection is found when comparing the abrupt4xCO2 simulations with the pre-industrial control (Figure S6 in Supporting Information S1).In contrast, weak and insignificant negative correlations between UNECK and ∆U 10 are found in CMIP5 (Figure 10b).
Any linkage between a past quantity and a future projected change offers the potential for an emergent constraint to narrow down future projections.While we must be cautious as to the robustness of the linkage between the past zonal winds in the neck region between the tropospheric westerlies and the SPV, given that it is not present in CMIP5 as shown in Figure 10a (Hall et al., 2019), it is worth considering whether this relationship has potential to constrain future projections.
We show that the observation-based estimates of UNECK from several reanalyzes (Figure 10d) sits within the model distribution but there is a tendency for more models to exhibit a positive bias than a negative bias in CMIP6.Using the emergent constraints method of Simpson et al. ( 2021), the right hand panel of Figure 10d shows an estimate of the constrained distribution based on the observed values of UNECK and the CMIP6-derived relationship between UNECK and ∆U 10 (we show the same for CMIP5 for comparison but such a constraint should not be used since the expected relationship between UNECK and ∆U 10 is not found there).The most likely adjustment from the constraint is slightly more negative ∆U 10 than the CMIP6 ensemble mean, but since the uncertainty is very large and the observed UNECK does not lie near the edge of the model distribution, it does not provide a particularly strong constraint on the future projections of ∆U 10 .Nevertheless, if this relationship between UNECK and ∆U 10 remains in future multi-model intercomparisons, and if uncertainties can be reduced for example, through the use of larger ensembles, then there is scope for further considering this relationship as a method of narrowing down model uncertainty in ∆U 10 .

Quantifying the Contributions to Projected SPV Changes
How much of the intermodel spread in projected SPV changes can be explained using the factors represented in Figure 1 and considered in more detail in the previous sections?Table 1 showed that in CMIP6 nearly half of the spread (r = 0.68) can be attributed to ∆F net , that is, stationary wave driving in the polar stratosphere.Here we attempt to quantify the different pathways affecting SPV changes as summarized in the form of a causal network in Figure 1.
We first build a multiple linear regression model (MLR) to explain the projected SPV changes (∆U 10 , box 9 in Figure 1) using upward (∆F 100 , box 7 in Figure 1) and meridional (∆F 54N , box 6 in Figure 1) wave propagation as predictors.The results of the calculations for CMIP6 are shown in Figure 11a and summarized in Table 4.For CMIP6, the MLR predicted and actual SPV changes are strongly correlated (r = 0.72, MLR #1) suggesting that ∆F 100 and ∆F 54N together can indeed explain roughly half (r 2 ∼ 0.5) of the projected SPV spread and confirming that they both contribute to ∆F net .Adding ∆F 10 does not improve the model (see Table 4, MLR #2), suggesting that the upper stratospheric flux change is not independent from ∆F 100 and ∆F 54N .For CMIP5 (Table S3 in Supporting Information S1), the correlation between actual ∆U 10 and those predicted by MLR with ∆F 100 and ∆F 54N as predictors is only r = 0.39, thus explaining only 15% of the ∆U 10 variance.The correlation between ∆U 10 and ∆F net (r = 0.50) in CMIP5 can only be recovered when all the three fluxes (∆F 100 , ∆F 10 , and ∆F 54N ) are included in the MLR, confirming that ∆F 10 contributes significantly to SPV changes in CMIP5, but not in CMIP6.Note that while the wave fluxes explain a considerable fraction of ∆U 10 spread, they are not completely independent of ∆U 10 because the zonal mean flow also affects wave propagation (as indicated by the bidirectional links in Figure 1).
We will next quantify how much the assumed drivers of stratospheric wave propagation, that is, boxes 1-5 in Figure 1, as discussed in Sections 3.2 and 3.3, can explain the projected SPV changes.We focus on discussing the  results for the effects of the subtropical wind response (∆U 100 , box 2 in Figure 1) and the past climatology of neck region winds (UNECK, box 8 in Figure 1).Moreover, we include the global-mean near-surface warming (∆T2m, box 1 in Figure 1) as a predictor for the effects not mediated through ∆U 100 and UNECK.Note that we tested all the factors discussed in Sections 3.2 and 3.3 and found that differential warming in the Equatorial Pacific (∆N3-N4, box 3 in Figure 1) and the Barents-Kara Sea ice loss (-∆SIC BKS , box 4 in Figure 1) do not provide additional information.Results of the MLR calculations for CMIP5 are given in Table S3 in Supporting Information S1, but because of the lack of correlations between these drivers and ∆U 10 in CMIP5, the following discussion is focused on CMIP6.
A regression model including ∆T2m, ∆U 100 and UNECK (MLR #3, Table 4) has all three regression coefficients significant, which shows that a larger global warming response stronger climatological UNECK winds lead to a strengthening SPV response, while a stronger subtropical jet stream response leads to a weakening SPV response.This model (r = 0.67) explains nearly as much of the spread in ∆U 10 (Figure 11b) as when using the EPfluxes as predictors (MLR#1, Figure 11a).These results are robust to overfitting, which we tested using "leaveone-out" cross-validation.More precisely, we iterated through all CMIP models and calculated regression coefficients using all but this one left-out model, and used the regression coefficients to reconstruct ∆U 10 for the left-out model.The correlation coefficients between these reconstructed ∆U 10 and CMIP6-predicted ∆U 10 are r = 0.66 for MLR#1 (with the 95% confidence intervals ranging from 0.45 to 0.80) and r = 0.56 for MLR#3 (the 95% confidence intervals 0.31 to 0.73), which are only slightly smaller than the correlation coefficients obtained by fitting all models.
To understand the effect of ∆U 100 (box 2) on ∆U 10 , we regress out ∆T2m (box 1), which is assumed to act as a common driver (see Figure 1).The regression coefficients for ∆T2m and ∆U 100 (MLR #4, Table 4) are both significant and opposite in sign.In other words, if interpreted causally based on Figure 1, this means that strengthening of the subtropical jet (U 100 , box 2) due to global warming leads to a weakening of SPV (U 10 , box 9), while the remaining effects of global warming lead to a strengthening of SPV.Importantly, the opposite sign relationships can explain the overall low correlation between ∆T2m and ∆U 10 (r ∆T2m; ∆U 10 = 0.07).This correlation can be recovered using the regression coefficients from MLR#4 and following the path-tracing rule (Kretschmer et al., 2021;Wright, 1921), that is, adding the products of the regression coefficients along the different pathways: Here the first summand is the estimated strength of the pathway from ∆T2m to ∆U 10 which is not mediated through ∆U 100 , and the second summand is that mediated via ∆U 100 (for an overview, see Kretschmer et al., 2021).Put differently, if from ∆T2m one subtracts its (standardized) contribution to SPV via the subtropical jet response, ∆T2m -(-0.74)× ∆U 100 , one gets a correlation of these residuals with ∆U 10 as high as r = 0.58 (in contrast to the low correlation of ∆T2m and ∆U 10 of only 0.07).The regression coefficients for ∆T2m and ∆U 100 in MLR#4 are smaller compared to those in MLR#3 because we exclude UNECK which is assumed to affect ∆U 10 through both ∆F 100 and ∆F 54N .
As a way to evaluate the assumed network in Figure 1, we also test if it is consistent with the data.For example, when regressing SPV changes (∆U 10 , box 9) on ∆F 54N , ∆F 100 and ∆U 100 (boxes 6, 7, 2) we obtain a small insignificant regression coefficient for ∆U 100 (not shown), which is consistent with its causal effect to be fully mediated by ∆F 54N and ∆F 100 (i.e., there is no direct link from box 2 to box 9 in Figure 1).Similar results are obtained when combining ∆T2m, UNECK, ∆N3 -N4 and -∆SIC BKS (boxes 1, 8, 3, 4) with the wave flux component responses (boxes 6 and 7) as predictors of ∆U 10 (not shown).Note again, that both ∆F 54N and ∆F 100 depend on ∆U 10 , therefore interpreting their correlation in a causal sense should be done with caution.
Overall, the results show that nearly half of the spread of the projected SPV changes in CMIP6 can be explained from different assumed pathways as summarized in Figure 1.This includes (a) a weakening of the SPV due to subtropical jet stream strengthening and increases in wave activity fluxes into the stratosphere, as well as (b) strengthening of the SPV due to equatorward wave activity refraction, which is sensitive to the background flow in the neck region and also depends on the level of the global warming.Barents-Kara sea ice loss in autumn might contribute to a weakening of the SPV, but the end of the century change is likely dominated by other factors given that the sea ice is gone in autumn in most models (e.g., Årthun et al., 2021).Differential sea surface temperature responses across the Pacific basin can drive anomalous extratropical atmospheric wave responses (Figure 8b), but these patterns do not align with the wave patterns associated with the SPV spread (Figures 7b and Table 3).

Discussion and Conclusions
Constraining future changes in SPV and reducing the large uncertainty projected by climate models remains challenging.The present paper makes several important steps toward better understanding the mechanisms contributing to future SPV changes and sources of the uncertainty.
First, we show that changes in stationary planetary wave fluxes are the dominant factor contributing to the spread in SPV response in CMIP6.We find that a half of the intermodel spread in SPV change is explained by the changes in the convergence of the stationary planetary waves in the polar stratosphere.This result contrasts to that found for CMIP5 where only a quarter of the spread in SPV response could be attributed to the response in the EPflux convergence (also found in W19).
Further, we find that both meridional and vertical propagation of the planetary waves is important for explaining the spread in SPV changes in CMIP6.Weakening of the SPV is mostly driven by increased upward propagation of the stationary planetary waves from below.Statistical analysis shows that a key factor in explaining the enhanced upward wave propagation and weakening SPV response is a strengthening of the subtropical jet stream, which increases the amplitude of the tropospheric planetary wavenumber 1 mostly by shifting and strengthening the North Pacific low.
On the other hand, we find that, in CMIP6, an increased equatorward wave propagation related to a weakened subtropical wave propagation barrier in the future (Sigmond & Scinocca, 2010) leads to a strengthening of the SPV.The changes in the equatorward wave propagation partly depend on the model's basic state climatology.Specifically, aligned with the arguments of Sigmond and Scinocca (2010), CMIP6 exhibits a significant link between the climatological strength of the zonal wind in the neck region and the SPV response as a result of climatological differences in refractive index.Idealized dynamical core model experiments demonstrated that tropical upper tropospheric warming beyond a certain threshold leads to SPV strengthening, and that this threshold is dependent on model resolution and background climatology (Walz et al., 2023;Wang et al., 2012).The mechanistic pathway in the idealized model of initially more equatorward wave propagation due to subtropical jet strengthening is consistent with the results found here; however, the degree to which enhanced wave reflection, which leads to a regime transition to a strong SPV in the idealized models (Walz et al., 2023) plays a role for comprehensive models can only be verified from daily EP fluxes, which are not available for the CMIP models.
By employing the causal network analysis, we also uncover a positive correlation between an SPV strengthening response and the global-mean near-surface warming response, which appears after controlling for the assumed pathway through the subtropical jet.This association might at least partly represent a thermodynamic mechanism, for example, a thermal wind contribution due to the strengthening meridional temperature gradient near the tropopause.An alternative explanation is that it represents a weakening of the wave propagation barrier and so contributes to increased equatorward wave refraction and thus strengthening of the SPV.The causal network analysis used here is a helpful framework in unpacking and interpreting correlations of projected changes (e.g., between ∆T2m and ∆U 100 ).However, it requires expert knowledge regarding the choice of considered mechanisms, thus, conclusions are conditional on the network quality and the data used to describe the involved mechanisms (Kretschmer et al., 2021).
The combined influence of the three above factors (subtropical jet response, global warming response and past climatological wind biases) explains nearly the same fraction of the SPV uncertainty as the stationary planetary wave fluxes.While this result is an important step toward understanding the SPV response, it does not answer the important question of why models disagree even on the sign of the SPV change.As the competing processes push the change to opposite directions, akin to a tug of war, the winning side might be determined by factors that are difficult to detect.In CMIP6, the prediction at least partly depends on the model's basic state climatology; however, used as an emergent constraint this relationship is not strong enough to significantly reduce the projected uncertainty in the SPV response and questions remain over why such a relationship is not present in CMIP5.
A major outstanding question is why the SPV projections in CMIP5 models are substantially less constrained by the changes in the planetary waves.While earlier studies proposed an additional role of gravity waves, our results show that their direct contribution is small, but their influence may be manifest indirectly through the stationary planetary wave response.Specifically, changes in the gravity wave drag lead to changes in the planetary wave drag, so that most of the gravity wave influence on zonal mean flow is mediated by the stationary planetary wave response (Sigmond & Shepherd, 2014).Although we find little evidence for a significant direct role of gravity waves for the SPV spread in CMIP6, we can not completely rule out their significant role in CMIP5.The other possible candidate is the influence of the sponge layer implemented in low-top models, which constitute about a half of CMIP5 models (Charlton-Perez et al., 2013).However, restricting the analysis to CMIP5 high-top models does not change the result.Further, we find evidence for a stronger downward influence of the circulation changes in CMIP5 than in CMIP6 (Figure 7), which might mask the upward influence.The reason for this difference between the ensembles is not clear.
Another potential source of uncertainty is the projected rate of the sea ice loss which varies strongly across both CMIP5 and CMIP6 models (Notz and SIMIP Community., 2020).The uncertainty in the atmospheric response to sea ice loss can also play a role as the sea ice loss not only contributes to generation of planetary waves but also affects refractive properties of the atmosphere, and so the propagation of the atmospheric waves (Smith et al., 2022).
It is also important to point out that the statistical approach adopted in the present paper has its limitations.Analysis of the advanced dynamical diagnostics requested as a part of DynVarMIP has helped, in particular, to better understand the role of the gravity waves.However, since only a limited number of the modeling centers provided DynVarMIP diagnostics, they could not be applied to examine the full spectrum of the model response.
We thus encourage the modeling centers to produce the DynVarMIP diagnostics for the next CMIP efforts.Coordinated multi-model experiments designed to test the mechanisms discussed in this paper are required to improve understanding.
and use the DJF differential temperature anomalies of the eastern Pacific with respect to the western Pacific calculated as the change in Niño 3 index minus Niño 4 index (∆N3-N4).Note that Oudar et al. (2020) used the reversed index, Niño 4 minus Niño 3. Finally, to assess the influence of the basic state winds on SPV change we follow W19 and use DJF basic state zonal mean winds calculated here at 70 hPa, 45°N-55°N and over the past 1980-2009 period (UNECK).

Figure 1 .
Figure 1.Assumed causal pathways (black arrows) and factors (boxed) on how global warming (red box 1) affects the stratospheric polar vortex response (blue box 9) considered in this study.Processes triggered by global warming are shown in green boxes.The assumed impact of model basic state zonal wind biases in the neck region (UNECK) is represented in purple (box 8).Variables used as a diagnostic for the specific factor are given in parentheses in each box and are defined in Section 2.

Figure 2 .
Figure 2. (a) Latitude-pressure cross sections of DJF multi-model mean historical climatology (contours, 1980-2009) and change following the SSP5-8.5 scenario (shading, 2070-2099 minus 1980-2009) in the zonal mean zonal winds (u); Panel (b) same as panel (a) but for the zonal wind tendencies due to ∇ • F → , Panel (c) same as panel (a) but for the zonal wind tendencies due to GWD; Panel (d) same as panel (a) but for the zonal wind tendencies due to the Coriolis torque.Only 9 CMIP6 models with the momentum budget data (Table S1 in Supporting Information S1) are used.The gray bar at 10 hPa indicates the region used for defining U 10 .Line contours mark multi-model mean climatologies of the respective variables: (a) 10, 20, 30, 40 m/s; (b-c) ±1, 2, 4, 6 m/s/day.Dashed contours are for negative values.Stippling shows where all 9 models agree on the sign of the response.

Figure 3 .
Figure 3. Correlation between changes in (a) U 10 and the tendency due to ∇ • F → , (b) U 10 and GWD, (c) U 10 and [ 1 ρ 0 a cos ϕ ∇ • F → + GWD] and (d) U 10 and the Coriolis torque across the subset of 9 models.Line contours mark multi-model mean climatologies of the respective variables (±1,2,4,6 m/s/day).Dashed contours are for negative values.Stippling shows correlation coefficients significant at p = 0.05 according to a two-sided t-test.

Figure 4 .
Figure 4.The change in DJF F → (arrows) and tendency of the zonal winds due to ∇ • F → (shading) together with changes in u (contours) in CMIP6 for (a) the stratospheric polar vortex (SPV) weakening group, (b) the SPV strengthening group and (c) their difference.The contour interval is 1 m/s and solid and dashed contour denote positive and negative values respectively.Numbers next to the red arrows show changes in F 100 (number below the blue box), F 10 (number on top of the blue box), and F 54N (number left of the blue box) in units of 10 4 kg m s 4 .The changes in F 100 and F 10 are defined as positive for upward fluxes, and in F 54N for equatorward fluxes.Numbers inside the box are the change in F net in units of 10 4 kg m s 4 .

Figure 5 .
Figure 5. Analogy of Figure 4 in W19.Scatterplots between model mean ∆U 10 and (a) ∆F net ; (b) ∆F 100 ; (c) ∆F 10 ; (d) ∆F 54N .Blue (red) dots depict stratospheric polar vortex weakening (strengthening) models.Solid (open) dots depict CMIP6 (CMIP5) models.The EP-flux change values here are in units of 10 4 kg m s 4 .Wind changes are in m/s.Diamonds are the model mean values.The correlation coefficients marked with stars are significant at p = 0.05.The correlation coefficients weighted by the ratio of the standard deviations, as in W19, are shown within the parentheses.

Figure 6 .
Figure 6.(a, c) DJF multi-model mean ∆Z′ 500 and (b, d) standard deviation across model mean ∆Z′ 500 between 1980-2009 and 2070-2099 in (a-b) CMIP6 and (c-d) CMIP5.Areas where 90% of the models agree on the sign of change are marked in (a, c) with black dots.Black contours show multi-model mean Z′ 500 climatologies in the historical period (1980-2009).The contours are ±50,100,150 m, with negative contours shown with dashed lines.The North Pacific (NP), North American (NA) and Northern European (NE) areas defined in Section 2 are marked with green boxes and labeled in (d).

Figure 7 .
Figure 7. Across-model regressions of DJF mean ∆Z′ 500 on (a, c) standardized ∆F 100 and (b, d) standardized ∆U 10 in (a-b) CMIP6 and (c-d) CMIP5.Areas where regression coefficients are significant at p = 0.05 according to a field significance test (Wilks, 2016) are highlighted with black dots.Black contours show multi-model mean Z′ 500 climatology in the historical period (contours ±50,100,150 m, negative contours shown with dashed lines).The North Pacific (NP), North American (NA) and Northern European (NE) areas discussed in the text are marked with green boxes and labeled in (d).

Figure 9 .
Figure 9. Composites of squared refractive index (n 2 ) in CMIP6 models simulating (a, c) weakening stratospheric polar vortex (SPV) and (b, d) strengthening SPV in the (a-b) past (1980-2009) and (c-d) future (2070-2099).The negative n 2 is depicted in black.Contours represent the corresponding zonal mean zonal wind composites (contour interval 5 m/s with the thick black line representing the 20 m/s contour).Only models projecting ∆U 10 larger than 2 m/s in the absolute sense (17 strengthening SPV models and 13 SPV weakening models) are included in the composites.

Figure 10 .
Figure 10.(a, c) The correlation between ∆U 10 and past zonal mean zonal wind for (a) CMIP5 and (c) CMIP6.The black horizontal line depicts the level and latitudes used to define the "neck region winds" (UNECK).(b-d) ∆U 10 versus past UNECK for (b) CMIP5 and (d) CMIP6.The gray cross hairs show 95% uncertainty ranges for the individual models based on large ensembles following Simpson et al. (2021).UNECK values are shown for 4 reanalyzes along with an uncertainty range centered on the ERA5 estimate, also based on uncertainties from large ensembles.In the right portion of the panel, the actual ∆U 10 values are shown along with a constrained uncertainty range obtained using the emergent constraints procedure of Simpson et al. (2021) using their Bayesian Hierarchical Model regression method.The colored bar shows the 95% confidence interval, the white horizontal lines show the 66% confidence interval, and the black horizontal line shows the mean constrained value.The dashed colored lines mark the best-fit regressions.

Table 1 Correlation
Coefficients (r) Between the F → Diagnostic Changes Over the Polar Stratosphere Box and ∆U 10 Across CMIP5 and CMIP6 Models Note.The results are shown for DJF F → correlated with DJF ∆U 10 as well as with JFM ∆U 10 .Correlation coefficient in bold are statistically significant at p = 0.05.Correlation coefficients in italic are not significant.The 2.5th-97.5thpercentile range of bootstrap sample are shown for CMIP6 in parentheses.Also shown are across-model standard deviations (STD) in 10 4 kg m s 4 for each of the F → diagnostic.

Table 2
Across Model Correlation Coefficients in CMIP5 and CMIP6 Between ∆Z′ 500NP , ∆Z′ 500NA ∆Z′ 500NE , Amplitudes of Planetary Wave 1 and 2 at 500 hPa (∆Z′ 500W1 and ∆Z′ 500W2 ) and ∆F 100 Note.Values for CMIP5 are given in parentheses.Values in bold are significant at p = 0.05 according to a two-sided t-test.Values in italic are not significant.
All variables are standardized prior to use in MLR.Regression and correlation coefficients significant at p = 0.05 according to a two-sided t-test are shown in bold.Values in italic are not significant.