The future North Atlantic jet stream and storm track: relative contributions from sea ice and sea surface temperature changes

Using a novel set of coordinated simulations from four different models, the response of the wintertime (December–February) North Atlantic jet stream and storm track to prescribed sea surface temperature increases and sea ice loss is analysed and the underlying physical mechanisms investigated. Three out of the four models show a southward shift of the upper-level jet stream with an increase in jet speed over Europe, where the contribution of sea surface temperatures dominates over the effects of sea ice loss. However, the remaining model lacks the increase in jet speed over Europe, which originates from opposite responses of similar magnitude due to the future sea surface temperatures and sea ice cover. The jet stream responses are primarily driven by the change in the meridional temperature gradient and, as a consequence, baroclinicity. At the same time, momentum flux convergence acts as a secondary amplifying and dampening factor. The same three models see a significant eastward shift of the extratropical cyclone track density, which is equally driven by changes to sea surface temperatures and sea ice cover. A consistent feature across all models is a decrease in the frequency of extratropical cyclones in the Mediterranean. The responses of extratropical cyclones to future sea ice cover and sea surface temperatures do not exceed the inter-model climatological differences. Notable differences in the future response of the jet stream and storm track occur, and thus considerable uncertainty remains in how the European climate will respond to a warmer climate.

1 Introduction

Global warming due to rising greenhouse gas concentrations is accompanied by a warming of the average sea surface temperatures (SSTs) and loss of sea ice cover (SIC). In particular, the Arctic is severely impacted as it is currently warming up to 4 times faster than the global average, a phenomenon known as Arctic amplification (Rantanen et al., 2022), and is experiencing dramatic sea ice loss (Simmonds and Li, 2021). The sea ice loss strongly contributes to a locally enhanced warming and moistening of the lower troposphere (Screen and Simmonds, 2010). However, the effects of sea ice loss on the mid-latitude circulation is a topic of intense debate (Smith et al., 2022; Screen et al., 2022; Ye et al., 2023).

Mid-latitude weather is characterised by extratropical cyclones (ETCs), which are organised into storm tracks on climatological timescales. Another prominent atmospheric feature in the mid-latitudes is the eddy-driven jet stream, which acts as a guide for ETCs. It is critical to study the storm track and the jet stream together, as they are closely connected (Athanasiadis et al., 2010; Ronalds and Barnes, 2019; Ye et al., 2023). The response of the North Atlantic jet stream to a warmer climate shows considerable inter-model variability over Europe (Zappa et al., 2018; Oudar et al., 2020). Moreover, future changes to SSTs and SIC have been shown to have opposing influences on the jet stream. A poleward shift of the jet stream is associated with rising SSTs; meanwhile, SIC loss leads to an equatorward shift (Barnes and Screen, 2015; Screen et al., 2018; Yu et al., 2024). A similar pattern has been reported for the storm tracks (Yu et al., 2023). In particular, the North Atlantic jet stream response to sea ice is highly uncertain (Screen et al., 2018), especially in the Northern Hemisphere winter (Simpson et al., 2014; Hay et al., 2022). Furthermore, the changes to the storm tracks related to sea ice loss are tightly linked to the jet stream changes (Ye et al., 2023; Ronalds and Barnes, 2019), necessitating more studies simultaneously investigating jet stream and storm track changes.

A major source of uncertainty in future projections from fully coupled climate models is the differing amount of sea ice loss that different models predict (Notz and Community, 2020). Moreover, the region and magnitude of sea ice loss also have been shown to have a substantial impact on the atmospheric response. For example, links between sea ice loss and an intensification of the negative phase of the North Atlantic Oscillation on a climatological timescale have been identified by Screen (2017). Using an intermediate complexity model, McKenna et al. (2018) showed that moderate sea ice loss in the Atlantic (Pacific) sector leads to a negative (positive) Arctic Oscillation response; meanwhile, extensive sea ice loss in either sector leads to a negative Arctic Oscillation response.

The response of the mid-latitude tropospheric zonal wind, which is closely related to the jet stream, is proportional to the eddy momentum feedback (Smith et al., 2022; Screen et al., 2022). Smith et al. (2022) showed that models tend to underestimate responses of mid-latitude tropospheric zonal wind due to changes in SIC as the eddy momentum feedback is too weak in climate models compared to observations. It is important to note that the magnitude of the underestimation is strongly model-dependent. This limits the attribution of mid-latitude changes to differences in sea ice loss or model representation of the atmospheric interactions between high latitudes and mid-latitudes. To reduce the impact of different models' biases on the future response from future projections, efforts are made to simultaneously analyse outputs from numerous models (Eyring et al., 2016; O'Neill et al., 2016; Smith et al., 2019). However, this often limits studies to the use of multimodel means, which reduces the physical interpretability of future climate responses. On the other hand, studies focusing on the physical mechanism tend to employ one model, which have reduced accountability for model-related uncertainties (Levine et al., 2021; Dai and Song, 2020; Chemke et al., 2019).

This study is part of the “Climate Relevant interactions and feedbacks: the key role of sea ice and Snow in the polar and global climate system” (CRiceS, 2021) project. CRiceS aims to understand the role of the polar processes, such as feedback loops, in polar and global climates. This includes quantifying processes that drive interactions and teleconnections between the higher and lower latitudes. For this purpose, coordinated model experiments investigating the contributions of SST and SIC changes to the future climate response are performed. The CRiceS simulations are following previous studies on the effect of sea ice loss, which have performed model simulations with prescribed sea ice cover using coupled models (McCusker et al., 2017; Oudar et al., 2017) and atmospheric general circulation models (AGCMs) (Deser et al., 2010; Ogawa et al., 2018; Smith et al., 2019). These simulations, each encompassing a continuous 40-year simulation, are performed using four AGCMs (OpenIFS-43r3, EC-Earth3, CESM2, NorESM2) which are all forced with the same prescribed SST and SIC for historical and multiple future climate conditions (Naakka et al., 2024). The full set of simulations includes a baseline simulation with historical SSTs and SIC, future simulations where both SSTs and SIC are changed simultaneously according to different emission scenarios, and simulations where the SSTs and SIC are changed independently. Thus, the individual contributions of SIC and SST are obtained by leveraging the power of the full set of experiments – discussed further in Sect. 2.1. Prescribing either future SIC or SST while keeping the other at historical levels allows us to study the contributions of changes in SIC and SST in isolation, which is a limitation of fully-coupled climate simulations like CMIP6 (Coupled Model Intercomparison Project Phase 6) (Eyring et al., 2016). An additional benefit of prescribed SIC and SST is reduced internal variability, which improves the detection of future climate signals. To further increase the likelihood of detecting statistically significant responses, the experiments representing future climate conditions use stronger warming scenarios for SST and SIC compared to previous studies (Smith et al., 2019; Peings and Magnusdottir, 2014; Screen et al., 2012; Deser et al., 2010). Moreover, analysing carefully designed simulations from a relatively limited number of models, rather than all models from the CMIP6 archive, leads to a reduction in the uncertainties of the responses due to differences in the projected SIC and SST. Enabled by high output frequency of a wide selection of atmospheric variables, the CRiceS simulations permit one to examine the underlying physical mechanisms and to identify structural differences in physical mechanisms of the response to projected SIC and SST across the selected models.

This work aims to determine the wintertime response of the North Atlantic jet stream and storm track to changes in sea ice cover and sea surface temperatures and to quantify their relative contributions to future climate conditions with simultaneous changes to SST and SIC. Specifically, the following questions are investigated:

• How well do the four models agree on the climatologies of the historical climate simulation when SSTs and sea ice cover are identical? The focus is on the winter North Atlantic jet stream, storm track, and individual extratropical cyclones.

• Are there differences across models in the total future climate response when SSTs and sea ice cover are changed simultaneously?

• What are the relative contributions of SST and sea ice cover changes to the total future climate response for each individual model? In particular, do the contributions oppose or amplify each other?

• What are the physical mechanisms leading to the responses due to changed SSTs and sea ice cover?

The study is structured as follows. Section 2 presents the conducted model simulations and introduces the metrics used to study the North Atlantic jet stream and storm track. Section 3 compares the North Atlantic jet stream and storm track across models in current climate conditions. Section 4 addresses the response of the North Atlantic jet stream in future climate and the contributions of SST and SIC. The physical mechanisms underlying the responses are explored in Sect. 5. Section 6 focusses on the responses of ETCs, including the contributions of SSTs and SIC. Specifically, it investigates the storm track density and multiple important ETC-specific metrics like lifetime, maximum intensity, and more. The study is concluded in Sect. 7.

2 Data and methods

2.1 Models and simulations

This study utilises a set of simulations by the CRiceS consortium, consisting of atmosphere-only simulations by three Earth system models (EC-Earth3, CESM2, and NorESM2) and one general circulation model (OpenIFS-43r3). The conducted simulations were designed with the aim of studying the impacts of changes in sea surface temperatures (SSTs) and sea ice cover (SIC). This is achieved by running experiments prescribing different combinations of historical and future SIC and SSTs. While previous studies used multimodel means of SST and SIC to force their simulations (Smith et al., 2019; Deser et al., 2010), in the CRiceS simulations the SST and SIC boundary conditions are obtained from the Australian Earth system model ACCESS-ESM1.5 from the CMIP6 archive (Eyring et al., 2016). ACCESS-ESM1.5 produces an Arctic sea ice cover evolution for the historical period that is in reasonable agreement with observations and provides the best guess estimate for future SIC (Notz and Community, 2020). Monthly climatological means from 1950–1969 are used as historical boundary conditions for the seasonal cycle, which is annually repeated. Meanwhile, the future boundary conditions use output from 2080–2099 under either the shared socioeconomic pathway SSP 1-2.6 or SSP 5-8.5 scenario. The full simulation set consists of a baseline (BL) simulation (historical SST and SIC), two future simulations (FT_x) (future SST and SIC), two SST simulations (ftSST_x) (future SST and historical SIC), and two SIC simulations (ftSIC_x) (historical SST and future SIC), where x is either SSP126 or SSP585 and denotes which socioeconomic pathway the boundary conditions are taken from. The simulation set is summarised in Table 1. Each simulation was run for 40 years plus 1 year of spin-up, starting on the 1 January. The variables related to dynamics and thermodynamics are saved as 6-hourly output, enabling more diagnostics to be computed offline. A more detailed description of the models and simulations, including a basic meteorological analysis of temperature, mean sea level pressure, and precipitation, is found in Naakka et al. (2024).

Table 1 The experiment names of the CRiceS coordinated simulation set performed by OpenIFS-43r3, EC-Earth3, CESM2, and NorESM2. The rows correspond to the prescribed sea surface temperature (SST), while the columns correspond to the prescribed sea ice cover (SIC). SSP refers to the shared socioeconomic pathways from O'Neill et al. (2016). This study uses the BL and SSP 5-8.5 simulations (shown in bold).

The CRiceS simulation set bears close resemblance to the set-up used by the Polar Amplification Model Intercomparison Project (PAMIP) (Smith et al., 2019), yet a few notable differences exist. While PAMIP uses a 100-member ensemble of 1-year simulations, the CRiceS simulations are comprised of 40-year continuous simulations to enable analysis using a continuous time series. However, using PAMIP simulations, Peings et al. (2021) have shown that while the local thermal response to Arctic sea ice loss is very consistent across the different 100-member ensembles, the mid-latitude circulation response differs significantly. Therefore, results need to be carefully interpreted when isolating the response of mid-latitude circulation to sea ice loss at +2 K global warming. In the CRiceS simulations, the SST and SIC boundary conditions for SSP 1-2.6 correspond to +1.82 K global warming and for SSP 5-8.5 correspond +4.4 K warming compared to the baseline simulation. The SSP 5-8.5 simulations provide a stronger SIC and SST forcing (Fig. A1) compared to PAMIP simulations, which improves the signal strength compared to internal variability and therefore facilitating a more robust detection of responses due to SIC and/or SST changes.

This work analyses the 39 available complete Northern Hemisphere (NH) winters (DJF) using the SSP 5-8.5 (FT_SSP585, ftSST_SSP585, ftSIC_SSP585, and BL) simulation. The analysis is limited to the North Atlantic (NA) sector, defined as 95° W–45° E and 20–80° N, with a focus on Europe (15° W–35° E and 30–70° N), given the high population density in this region and the potential for large societal impacts if the jet stream and ETC characteristics change in the future. The model output is vertically interpolated from model levels to isobaric surfaces from 1000 to 50 hPa in 50 hPa intervals. The horizontal resolution is kept at the native resolution of each model. OpenIFS-43r3 and EC-Earth3 are run with TL255 horizontal resolution (0.7° × 0.7° at the Equator) and 91 vertical model levels. NorESM2 has a longitude–latitude resolution of 2.5° × 1.9°, while for CESM2 the resolution is 1.25° × 0.9°. NorESM2 and CESM2 are run with 32 model levels. The selected models can be grouped into two families based on their atmospheric component. The Integrated Forecast System (IFS), developed by the European Centre for Medium-Range Weather Forecasts, is the basis for EC-Earth3 (IFS cycle 36r4) and OpenIFS-43r3 (cycle 43r3). Similarly, the atmospheric component of NorESM2 was developed from the atmospheric component of CESM2, namely the Community Atmospheric Model version 6 (CAM6).

The responses are calculated as the difference in the December to February climatological means between a perturbed simulation and the baseline simulation (BL). The mathematical formulation is as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathrm{\Delta }\mathrm{FT}={\mathrm{FT}}_{\mathrm{SSP}\mathrm{585}}-\mathrm{BL},\\ & \mathrm{\Delta }\mathrm{SST}={\mathrm{ftSST}}_{\mathrm{SSP}\mathrm{585}}-\mathrm{BL},\\ & \mathrm{\Delta }\mathrm{SIC}={\mathrm{ftSIC}}_{\mathrm{SSP}\mathrm{585}}-\mathrm{BL}.\end{array}\end{array}$$

Note the scenario subscripts are omitted in the name of the response as this study only uses the SSP 5-8.5 scenario. The difference between the FT_SSP585 and baseline climatology is referred to as future response, denoted by ΔFT. The response to changes in sea surface temperature, ΔSST, is calculated from the difference between the SST simulation ftSST_SSP585 and baseline simulation BL, and the sea ice response, ΔSIC, is calculated likewise. It is important to highlight that the summation of ΔSIC and ΔSST does not result in ΔFT. The differences arise from the non-linear interactions when SIC and SSTs are changed simultaneously. Furthermore, in the SIC simulations, when sea ice is removed SSTs are not changed. The lacking change in SSTs in these locations contributes to the differences, but it has been shown to be a minor effect (Singarayer et al., 2006; Deser et al., 2010). Additionally, the internal variability also contributes to the differences between ΔFT and ΔSST+ΔSIC. For a more in-depth analysis of the non-linear interactions in the CRiceS simulation set, consult Naakka et al. (2024).

To test if the climatological means of perturbed experiments are statistically different from the baseline mean climatology for gridded data, a two-tailed t test with the 39 seasonal means for DJF as input was performed. Subsequently, the significance was controlled using false discovery rate (FDR) according to Wilks (2016). Furthermore, the consistency discovery rate (CDR), proposed by Peings et al. (2021), was applied to test for a significant consistent sign in the response. CDR testing in this study takes the following steps. We took 39 responses (ΔFT, ΔSST, ΔSIC) from the seasonal DJF means. A subsample of 20 from the 39 responses were drawn without repetition. The climatological average was calculated, and the sign was recorded. The subsampling was repeated 1000 times. If 900 iterations agree on the sign, the sign of the response was considered significantly consistent. In the figures, shading indicates the responses which show a CDR consistent response, and areas without stippling show statistically significant responses with the FDR corrected t test. Significant differences in ETC quantities described in Sect. 2.3 are detected using a Mann–Whitney U test with a threshold p value < 0.05.

2.2 Baroclinicity and momentum flux convergence

The strength of the NA jet stream is predominantly governed by two physical mechanisms. One is the thermal wind law, which describes vertical changes in the wind speed due to a horizontal temperature gradient. Relevant for the jet stream is the negative meridional temperature gradient (e.g. temperature decreases toward the poles) in the troposphere, which leads to an increasing zonal wind with height. The other mechanism originates from eddy–mean-flow interaction, whereby the momentum of atmospheric waves (eddies) is fed back to the jet stream (Eliassen and Palm, 1961; Hoskins et al., 1983).

The necessary metrics to study the mechanisms driving the NA jet stream require eddy–mean-flow separation, which is a common technique to study synoptic-scale atmospheric dynamics. This is achieved by low-pass filtering using a 21-weight Lanczos filter with a 10 d cut-off period. The longer timescale of 10 d compared to the more common 6 d is motivated by including breaking synoptic waves as part of the eddy flow (Rivière et al., 2018). For the filtering, the shoulder months (November and March) are added to the dataset and then subsequently discarded for analysis. The low-pass filter is applied to the potential temperature θ and wind components u and v. The eddy part of the wind components u^′ and v^′ is determined by subtracting the low-pass filter field from the full field.

In isobaric coordinates, the horizontal temperature gradient is proportional to the potential temperature gradient, which is quantified by the meridional component of the baroclinicity vector B_y. The baroclinicity vector, B_s, is given by

$$\begin{array}{}\text{(2)}& {\mathit{B}}_{\mathrm{s}}=-{\displaystyle \frac{\mathrm{\nabla }\stackrel{\mathrm{‾}}{\mathit{\theta }}}{\sqrt{S}}},\end{array}$$

where

$$S={\displaystyle \frac{\mathrm{1}}{h}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\theta }}}{\partial p}},$$

and

$$h={\displaystyle \frac{R}{p}}{\left({\displaystyle \frac{p}{p_{\mathrm{0}}}}\right)}^{R/c_p},$$

as used by Cai and Mak (1990) and Schemm and Rivière (2019). B_s is calculated from the low-pass-filtered potential temperature $\stackrel{\mathrm{‾}}{\mathit{\theta }}$. Additionally, the temperature gradient $\mathrm{\nabla }\stackrel{\mathrm{‾}}{\mathit{\theta }}$ is normalised by the stability S, which uses the scale height h, which includes the gas constant for air R, specific heat capacity of air at constant pressure c_p, and reference pressure p₀=1000 hPa. The zonal wind speed and the Eady growth rate are proportional to the meridional component of the baroclinicity B_y.

Multiple approaches have been developed to characterise the eddy–mean-flow interaction and its effects on the jet speed, most notable being the Eliassen–Palm flux (Eliassen and Palm, 1961) and E vectors (Hoskins et al., 1983) – and many more exist (Trenberth, 1986; Plumb, 1985). The Eliassen–Palm flux uses the zonal averages and quantifies eddies as deviations from the zonal average. On the other hand, E vectors utilise all three spatial dimensions, where the eddies and mean flow are obtained by using temporal filters. Both methods are physically motivated by the transfer of momentum between eddies and mean flow. Additionally, the mathematical form is identical for zonal averages (Trenberth, 1986). The present article uses the meridional component F_ϕ of the E vectors, defined as

$$\begin{array}{}\text{(3)}& F_{\mathit{\varphi }}=a\cos \mathit{\varphi }\left(-u^{\prime }{v}^{\prime }\right),\end{array}$$

where a is the radius of Earth and ϕ is the latitude. The variables u^′ and v^′ symbolise the eddy in terms of the horizontal wind components. Specifically, the momentum flux convergence,

$$\begin{array}{}\text{(4)}& \mathrm{MFC}=〈{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial F_{\mathit{\varphi }}\cos \mathit{\varphi }}{\partial \mathit{\varphi }}}〉,\end{array}$$

is calculated and results in the zonal mean, denoted by 〈⋅〉. The momentum flux convergence (MFC) is proportional to the acceleration of the zonal wind speed u for frictionless motion, and it quantifies the feedback of eddies on the jet. A similar metric was used by Smith et al. (2022) to assess the zonal wind speed response to changes in Arctic sea ice.

While the temperature gradient and eddy–mean-flow interaction mechanisms are not directly quantitatively comparable, they are tightly connected in the jet stream area. By comparing the sign of changes in baroclinicity and momentum flux convergence, a qualitative interpretation of compounding or opposing mechanistic effects on the responses of the zonal wind due to sea ice and SST changes is possible.

2.3 Extratropical cyclone tracking

Extratropical cyclones (ETCs) are objectively identified by the TRACK algorithm (Hodges, 1994, 1999). The tracking is performed using mean sea level pressure, as it is available as direct model output for all four models. The mean sea level pressure is truncated to T63 resolution with wave numbers 5 and lower removed, from which local minima are identified as ETCs. Like Priestley et al. (2023), the individual ETCs are filtered according to multiple criteria: they have to be (1) mobile (travel at least 1000 km), (2) long-lasting (have a lifetime of at least 48 h), (3) affect Europe (at least 48 h within the European box as defined in Sect. 2.1), and (4) and occur in the NH winter (genesis date in DJF).

Multiple quantities are derived from the TRACK output for each entire ETC track satisfying conditions (1), (2), (3), and (4). Track count originates from the number of tracks, track duration is obtained from the number of time steps, genesis latitude corresponds to the latitude at the first time step, and latitudinal displacement is the difference in latitude between the last and first time step. The mean speed is obtained by calculating the speed at each time step from the change in coordinates and subsequently averaged over the track. The maximum intensity is the maximum of the negative pressure anomaly from the T63 mean sea level pressure within a single ETC track. For ETCs satisfying conditions (1), (2), and (4), the track density is computed using spherical kernels (Hodges, 1996), and the likelihood that a given point is affected by an ETC is estimated.

3 Intercomparison of the baseline simulation

This section compares the baseline simulation across models. For readability, a detailed description of the baseline simulation of OpenIFS-43r3 is given in each subsection. Subsequently, the key differences between the three remaining models and OpenIFS-43r3 are presented.

3.1 The horizontal perspective

The NA jet stream, identified using the 250 hPa zonal wind speed, originates at the North American east coast and ends over Europe. In OpenIFS-43r3 (Fig. 1a, colours), the NA jet stream is strongest (36 and 40 m s⁻¹) over the western NA between 30–40° N. Toward the east, the jet is narrower and weaker, and the jet core moves northward, commonly referred to as the tilt of the NA jet. The eddy-driven jet, identified via the 850 hPa zonal wind as in Woollings et al. (2010), is strongest (9 to 12 m s⁻¹) in the central Atlantic (Fig. 2a, contours), which is east of the upper-level maximum. Similar to the upper-level jet, the low-level jet shows a northward tilt from west to east. The structure of both the upper-level and lower-level NA jet stream in the OpenIFS-43r3 baseline simulation is in good agreement with what is found in reanalysis (Fig. A2) and historical simulations from coupled climate models (Harvey et al., 2020), despite the baseline simulation not being directly comparable to these as there is no interannual variability of SST and SIC.

Figure 1 The 39-year DJF mean of zonal wind speed u at 250 hPa (in colour shading, m s⁻¹) and baroclinicity B_y at 500 hPa (in black contours, 10⁻⁴ s⁻¹) for (a) OpenIFS-43r3, (b) EC-Earth3, (c) NorEMS2, and (d) CESM2.

Figure 2 The 39-year DJF mean of ETC track density (in colour shading) and zonal wind speed u at 850 hPa (in black contours, m s⁻¹) for (a) OpenIFS-43r3, (b) EC-Earth3, (c) NorEMS2, and (d) CESM2. The ETC track density unit is number of ETCs per 5° spherical cap per winter season (DJF).

The baroclinicity B_y is diagnosed at 500 hPa. The maxima and minima in B_y closely coincide with the maxima and minima in the 250 hPa zonal wind. In the NA (Fig. 1a, contours), the maximum in B_y is located slightly northward of the upper-level jet maximum, yet their respective shapes match closely. As was the case with the jet stream, B_y shows a northward tilt and an eastward decrease in magnitude.

The ETC track density is tightly linked to the presence of baroclinicity and the jet stream. Three strong maxima are identified (Fig. 2a, colours): one over the eastern USA, another at the North American east coast, and the last one southeast of Greenland. Generally, the ETC track density is highest in regions of high baroclinicity and on the northern flank of the lower-level and upper-level jet stream, thus following the northward tilt across the NA. Further prominent features are elevated ETC track density in the Mediterranean and Baltic Sea. The ETC densities are in good agreement (not shown) with reanalysis (Gramcianinov et al., 2020) and coupled climate models (Priestley et al., 2020).

Generally, the baseline simulation from EC-Earth3 (Figs. 1b and 2b) is similar to that from OpenIFS-43r3 when the zonal wind speed at 250 and at 850 hPa, baroclinicity, and ETC track density are considered. Minor differences are found in the ETC densities, which are 40 % higher in the northeastern USA and mainland Europe compared to OpenIFS-43r3. NorESM2 and CESM2 have comparable baseline simulations to OpenIFS-43r3 and EC-Earth3. However, there are multiple distinguishing features. Most notably, the low-level jet stream is considerably faster and oriented more zonally, and the maximum is shifted toward the east (Fig. 2c, d, contours) in both NorESM2 and CESM2 compared to OpenIFS-43r3 and EC-Earth3. Furthermore, the upper-level jet stream speed reduces more gradually toward the east, which results in higher jet speed over Europe (Fig. 1c, d, colours) in NorESM2 and CESM2. Moreover, the ETC track density is higher extending from the North American east coast over the southeast of Greenland to northern Europe (Fig. 2c, d, colours). Lastly, in NorESM2 and CESM2 the baroclinicity is lower outside the jet stream regions.

While NorESM2 and CESM2 share a similar pattern, distinguishing them from OpenIFS-43r3 and EC-Earth3, differences also exist between NorESM2 and CESM2. NorESM2 shows a 10 % stronger upper-level jet over Europe (Fig. 1c, colours) compared to CESM2 (Fig. 1d, colours). Another feature is the ETC densities over northern Europe, where CESM2 (Fig. 2d, colours) has higher values than NorESM2 (Fig. 2c, colours). This is related to higher ETC counts in CESM2 compared to NorESM2, which is discussed in Sect. 3.3.

3.2 The zonal cross-section perspective

This study's primary focus is on Europe (30–70° N, 15° W–35° E) for multiple reasons. First, Europe is a densely populated region where the jet stream and ETCs may lead to societal impact. Second, this is where the exit region of the NA jet is located, and the left-hand jet exit is an area of upper-level forcing that induces ascent, which is favourable for extratropical cyclone intensification. Third, the four models have shown differences in the zonal wind over Europe at 250 hPa (Fig. 1) and 850 hPa (Fig. 2). The cross-sections of the zonal mean are utilised in Fig. 3, where the zonal wind speed u (in colours), the baroclinicity B_y (in blue contours), and the eddy momentum flux convergence (MFC) (in black contours) are displayed. This allows for insight into the European jet properties across models.

Figure 3 Zonal mean (15° W–35° E) of the 39-year DJF mean of zonal wind speed u (in colour shading, m s⁻¹), eddy momentum flux convergence (MFC) (in black contours, m² s⁻²), and baroclinicity B_y (in blue contours, 10⁻⁴ s⁻¹) for (a) OpenIFS-43r3, (b) EC-Earth3, (c) NorEMS2, and (d) CESM2.

The OpenIFS-43r3 baseline cross-section is presented in Fig. 3a. The edge of the subtropical jet stream located over northern Africa (around 30° N and 200 hPa) is visible as high values of u. A second maximum (values between 16 and 18 m s⁻¹) in zonal wind speed at 250 hPa and between 50 and 60° N is associated with the jet exit region of the NA jet stream. An area of increased baroclinicity B_y (600 to 400 hPa) is found below and north of the NA jet maximum. In addition, this NA jet stream maximum is co-located with a maximum in MFC, which contributes to the formation of the NA jet stream.

EC-Earth3 (Fig. 3b) closely matches OpenIFS-43r3 in patterns and magnitudes of u, B_y, and MFC. On the other hand, in CESM2 (Fig. 3d), the NA jet stream is a vertically deep feature, with a strong meridional gradient in wind speed and has a maximum speed between 22 and 24 m s⁻¹. The MFC, with values between 100 and 125 m² s⁻², is considerably stronger compared to OpenIFS-43r3 and EC-Earth3. This indicates that the eddy–mean-flow interaction contributes substantially to the jet structure in CESM2. Meanwhile, the area of strong B_y below the jet maximum is vertically thicker compared to OpenIFS-43r3 and EC-Earth3, which contributes to a stronger jet maximum.

Lastly, NorESM2 (Fig. 3c) shares the vertically deep structure and a similar maximum NA jet speed as the CESM2 baseline simulation. The B_y shows small differences to CESM2, mainly a thicker layer exceeding 12 × 10⁻⁴ s⁻¹ at the latitude of the NA jet. Most notably, NorESM2 has the highest values in MFC (between 125 and 150 m² s⁻²). Of all four models, the jet stream in NorESM2 has the highest contribution from eddy momentum flux convergence.

3.3 The ETC properties

Looking at the properties of ETCs affecting Europe (defined as spending at least 48 h within 30–70° N, 15° W–35° E), Fig. 4 provides insight into the ETC count, track duration, mean speed, maximum intensity, the genesis latitude, and latitudinal displacement. All models agree well on the distribution's overall shape and the range of values across all six ETC quantities. However, multiple statistically significant key differences in the mean values are addressed in the following. The detailed differences in the mean values are reported in Table A1 in Appendix A.

Figure 4 The distribution of six ETC quantities and their mean values shown by individual markers for OpenIFS-43r3 (red, cross), EC-Earth3 (blue, square), NorESM2 (green, inverted triangle), and CESM2 (purple, circle). Given in brackets in the legend is the total number of ETCs during 39 DJF periods. The six ETC quantities include ETC count per winter season (DJF) (a), ETC lifetime (b), mean speed (c), maximum intensity (d), genesis latitude (e), and latitudinal displacement (f). The y axis gives the probability density. Note that the y axes are different between panels.

The most apparent differences are present in the ETC count. NorESM2 has a significantly lower number of ETCs affecting Europe (1550 total, 39.7 per winter season) than the three other models, while CESM2 has the highest number (1814 total, 46.5 per winter season). The mean values for OpenIFS-43r3 (1694 total, 43.4 per winter season) and EC-Earth3 (1720 total, 44.1 per winter season) are not statistically different from each other. Referring back to the track densities in Fig. 2, the higher number of ETCs in CESM2 is the origin of the higher ETC track density compared to OpenIFS-43r3 and EC-Earth3. However, regardless of the lower ETC count in NorESM2, NorESM2 presents a higher ETC track density compared to OpenIFS-43r3 and EC-Earth3. This results from the combination of low mean speed and long lifetime of ETCs in NorESM2, which amplifies the ETC track density in the storm track.

The mean lifetime of ETCs is similar in EC-Earth3 (120.6 h) and OpenIFS-43r3 (124.9 h), while it is significantly higher in CESM2 (132.5 h) and NorESM2 (148.8 h). Next, in NorESM2, the mean speed of ETCs with a value of 11.0 m s⁻¹ is significantly slower than in the other three models. Meanwhile, EC-Earth3 has the fastest ETCs on average (11.9 m s⁻¹). There is no statistical difference in mean speed between CESM2 (11.3 m s⁻¹) and OpenIFS-43r3 (11.2 m s⁻¹).

Furthermore, the mean maximum intensity is very similar in OpenIFS-43r3 (33.3 hPa), CESM2 (32.8 hPa), and NorESM2 (32.5 hPa) with no significant differences. EC-Earth3 is the exception, with a significantly higher mean maximum intensity of 35.3 hPa. Moreover, out of the four models, CESM2 has the most equatorward mean latitude of genesis of ETCs affecting Europe (39.4° N), and EC-Earth3 has the most poleward genesis location (41.8° N). These models differ significantly from OpenIFS-43r3 (40.1° N) and NorESM2 (40.6° N), with the difference between the latter two models being insignificant. Lastly, latitudinal displacement is significantly different in all models, in order from least to most poleward: EC-Earth3 (12.8° N), OpenIFS-43r3 (14.0° N), CESM2 (15.1° N), and NorESM2 (15.8° N).

4 The jet response to climate change

4.1 The total future response of the jet

The response of the upper-level jet stream to changed SSTs and SIC in the SSP5-8.5 scenario is shown in Fig. 5. All models show an apparent deceleration on the poleward side of the subtropical jet stream over northern Africa. This is due to an upward shift of the subtropical jet maximum (Fig. 6), driven by higher SSTs in the tropics and extratropics which leads to a warmer and deeper troposphere. Another common feature across all models is the reduced 250 hPa zonal wind speed on the poleward side of the NA jet stream southeast of Greenland.

Figure 5 Future response ΔFT in DJF mean of zonal wind speed u at 250 hPa (in colour shading) and the baseline climatology at 250 hPa (in black contours) for (a) OpenIFS-43r3, (b) EC-Earth3, (c) NorEMS2, and (d) CESM2. Shading indicates a consistent sign according to CDR testing, and stippling indicates statistically insignificant changes (FDR-corrected t test).

Figure 6 Future response ΔFT of zonal mean (15° W–35° E) of the 39-year DJF mean of zonal wind speed u (in colour shading) and the baseline climatology (in black contours) for (a) OpenIFS-43r3, (b) EC-Earth3, (c) NorEMS2, and (d) CESM2. Shading indicates a consistent sign according to CDR testing, and stippling indicates statistically insignificant changes (FDR-corrected t test).

The most notable and distinguishing aspects across the models are the changes in wind speed in the NA jet stream at 250 hPa (Fig. 5). In particular, OpenIFS-43r3 only shows a significant increase on the equatorward side of the jet located over the central Atlantic (Fig. 5a). The remaining three models form a group where the jet speed increases on the jet's southern side and in the jet exit region over Europe. To gain more insight into the disagreement between OpenIFS-43r3 and the remaining three models, Fig. 6 depicts the zonal mean cross-section over Europe (15° W–35° E), where the responses differ the most. Evidently, the non-existent response in the NA jet exit over Europe in OpenIFS-43r3 is also visible in the cross-section (Fig. 6a). The only statistically significant features are a deceleration of the polar vortex (60–70° N, 150–50 hPa) and an upward shift of the subtropical jet stream. In contrast, EC-Earth3 responds with a barotropic increase of the NA jet stream between 850–350 hPa and 50–60° N (Fig. 6b). Around the jet maximum at 200 hPa, there is a strong vertical shear in ΔFT, suggesting an additional baroclinic contribution. A similar response in the SSP 5-8.5 scenario to EC-Earth3 is found in NorESM2 (Fig. 6c) and CESM2 (Fig. 6d). However, the maximum increase in zonal wind speed u is located equatorward of the NA jet maximum in the respective baseline simulation. Furthermore, in NorESM2 and CESM2, the response increases in magnitude with height throughout the troposphere (850–200 hPa), indicating a more baroclinically driven response.

4.2 The contribution of SST and SIC

A key goal of the present paper is to investigate the contributions of SST and SIC changes to the combined climate response. The SST response, ΔSST, and SIC response, ΔSIC, of the 250 hPa zonal wind are presented in Fig. 7. EC-Earth3, NorESM2, and CESM2 show very similar structures in both contributions. ΔSST (Fig. 7c, e, f) closely resembles the future response (Fig. 5b, c, d), with an increase on the southern side and in the eastern exit region of the NA jet stream. Hence, the future response is largely dominated by changes in SSTs. The ΔSST (Fig. 7d, f, h) spatial pattern exhibits a tripole structure. There is a deceleration south of Greenland and around 20–30° N in the North Atlantic and an increase in zonal wind speed at 40–50° N. The exact geographical location and extent of the tripole structure vary between EC-Earth3, NorESM2, and CESM2. Generally, the increases and decreases in ΔSST and ΔSIC are co-located, resulting in an amplifying effect in the future response.

Figure 7 SST response ΔSST (a, c, e, g) and SIC response ΔSIC (b, d, f, h) of DJF mean of zonal wind speed u at 250 hPa (in colour shading) and the baseline climatology at 250 hPa (in black contours) for (a, b) OpenIFS-43r3, (c, d) EC-Earth3, (e, f) NorEMS2, and (g, h) CESM2. Shading indicates a consistent sign according to CDR testing, and stippling indicates statistically insignificant changes (FDR-corrected t test).

OpenIFS-43r3 is the outlier in ΔSST and ΔSIC regarding the NA jet stream at 250 hPa. The response to SST (Fig. 7a) has broadly the same spatial pattern as in the other models, but it is lower in magnitude. The maximum increase is between 4.0 and 4.5 m s⁻¹, while it exceeds 6.0 m s⁻¹ in the other models. However, the key difference setting OpenIFS-43r3 apart from the other models is identifiable in the SIC response (Fig. 7b). There is an absence of a tripole structure and a consistent decrease (determined using CDR) relative to the baseline simulation in the wind speed in the NA jet exit region over Europe. The SST-related increase and SIC-related decrease over Europe compensate each other, resulting in the lack of the response in ΔFT (Fig. 5a).

Figure 8 SST response ΔSST (a, c, e, g) and SIC response ΔSIC (b, d, f, h) of the zonal mean (15° W–35° E) of the 39-year DJF mean of zonal wind speed u (in colour shading) and the baseline climatology (in black contours) for (a, b) OpenIFS-43r3, (c, d) EC-Earth3, (e, f) NorEMS2, and (g, h) CESM2. Shading indicates a consistent sign according to CDR testing, and stippling indicates statistically insignificant changes (FDR-corrected t test).

The ΔSST and ΔSIC of zonal wind speed cross-section from the OpenIFS-43r3 simulations are shown in Fig. 8a and b. This reaffirms the compensating contributions of ΔSST and ΔSIC. Due to increased SSTs, the subtropical jet shifts upward, and the NA jet stream increases in speed between 300 hPa and 100 hPa. Meanwhile, ΔSIC shows a reduction in the NA jet speed throughout the troposphere (Fig. 8b). The change in SIC is also the origin of the deceleration of the polar vortex in the future response (Fig. 8a).

EC-Earth3, NorESM2, and CESM2 present structurally similar ΔSST for the cross-section of zonal wind speed. The SST responses (Fig. 8c, e, g) show an intensification of the NA jet stream equatorward of the jet core in the baseline simulation. All three models exhibit a positive response in the zonal wind speed through the vertical column, with the strongest response at upper levels (300–100 hPa). However, in CESM2 the response is stronger compared to ΔSST in NorESM2 and EC-Earth. For NorESM2 and CESM2, ΔSIC (Fig. 8f, h) is characterised by a reduction on the poleward side and an increase on the equatorward side. Lastly, EC-Earth3 shows no changes to the u cross-section in response to the changes in SIC (Fig. 8d).

5 The mechanisms of the jet response

5.1 The mechanisms for SST response

A feature common to all models is the impact of SSTs on baroclinicity driving the change of the subtropical jet (Fig. 9a, c, e, g, 30–35° N). The baroclinicity between 850 hPa and 300 hPa decreases, while an increase is visible from 300 to 100 hPa. This corresponds well with the upward shift of the subtropical jet described in Sect. 4.2. The origin of the increased baroclinicity above 300 hPa is an increased upper-level warming in the tropics (not shown) and an upper-level cooling in the potential temperature in the subtropics/mid-latitudes (Fig. 10a, c, e, g, 30–35° N, 200–100 hPa). The deepening of the troposphere due to warming SSTs and the subsequent lifting of the thermal tropopause explains the cooling. The tropospheric deepening, and subsequent increase in jet height, also explains the increase in momentum flux convergence located above the maxima in the respective baseline simulation for all models (Fig. 11a, c, e, g, 40–60° N, 300–100 hPa).

Figure 9 SST response ΔSST (a, c, e, g) and SIC response ΔSIC (b, d, f, h) of the zonal mean (15° W–35° E) of the 39-year DJF mean of baroclinicity B_y (in colour shading) and the baseline climatology (in black contours) for (a, b) OpenIFS-43r3, (c, d) EC-Earth3, (e, f) NorEMS2, and (g, h) CESM2. Shading indicates a consistent sign according to CDR testing, and stippling indicates statistically insignificant changes (FDR-corrected t test).

In OpenIFS-43r3 (Fig. 9a), the baroclinicity associated with the NA jet exit (45–65° N) increases at upper levels (400–200 hPa); this results in the shallow response of the zonal wind speed seen in Fig. 8a. The increase in baroclinicity is caused by a larger increase in potential temperature at lower latitudes (30–50° N) compared to higher latitudes (50–70° N). This meridional difference in potential temperature increase is the strongest between 400–300 hPa (Fig. 10a). MFC increases at the location of the zonal wind speed increase (45–65° N, 200–100 hPa). Like the zonal wind speed response, OpenIFS-43r3 shows the weakest response in MFC among all models.

Figure 10 SST response ΔSST (a, c, e, g) and SIC response ΔSIC (b, d, f, h) of the zonal mean (15° W–35° E) of the 39-year DJF mean of potential temperature θ (in colour shading) and the baseline climatology (in black contours) for (a, b) OpenIFS-43r3, (c, d) EC-Earth3, (e, f) NorEMS2, and (g, h) CESM2. Shading indicates a consistent sign according to CDR testing, and stippling indicates statistically insignificant changes (FDR-corrected t test). Note that different colour scales are used because the absolute values of ΔSST and larger than those of ΔSIC.

EC-Earth3 shows a very similar SST response to OpenIFS-43r3 with two distinct differences (Fig. 9c). First, the upper-level baroclinicity increase (45–65° N, 400–200 hPa) is stronger. This is the consequence of a stronger upper-level warming in EC-Earth3 (Fig. 10c, 30–40° N, 400–200 hPa). Furthermore, the enhanced baroclinicity response explains the stronger increase in the upper-level part of the NA jet stream exit in EC-Earth3 compared to OpenIFS-43r3, as shown in Sect. 4.2. Second in EC-Earth3, there is a weak increase in baroclinicity from 40 to 50° N and between 700 and 400 hPa, which results from a more pronounced warming at the subtropical middle troposphere (30–40° N, 400–700 hPa). However, this increase in baroclinicity does not explain the increase in zonal wind speed at lower levels (500–850 hPa) in Fig. 8c, as the maximum is located between 49–54° N. This requires the momentum flux convergence (Fig. 11c), which shows an increase at the lower levels in EC-Earth3 (47–53° N, 500–850 hPa). An increase in MFC leads to a barotropic increase in zonal wind speed, present in Fig. 8c. Furthermore, there is a decrease in MFC between 60 and 70° N at upper levels (400–200 hPa), which is co-located with an increase in baroclinicity. These effects are compensating, resulting in a reduced increase in zonal wind speed co-located with a MFC decrease.

Figure 11 SST response ΔSST (a, c, e, g) and SIC response ΔSIC (b, d, f, h) of the zonal mean (15° W–35° E) of the 39-year DJF mean of momentum flux convergence (MFC) (in colour shading) and the baseline climatology (in black contours) for (a, b) OpenIFS-43r3, (c, d) EC-Earth3, (e, f) NorEMS2, and (g, h) CESM2. Shading indicates a consistent sign according to CDR testing, and stippling indicates statistically insignificant changes (FDR-corrected t test).

The ΔSST of baroclinicity in NorESM2 (Fig. 9e) is structurally similar to EC-Earth3. The strongest increases in baroclinicity are found at 400 hPa to 200 hPa and from 40 to 60° N. Additionally, NorESM2 also features a weak but significant increase in baroclinicity at lower levels (500–850 hPa, 35–50° N), which is more pronounced relative to EC-Earth3. Contrary to EC-Earth3, the origin of the baroclinicity response in NorESM2 is slightly different. Figure 10e shows that the upper-level warming between 30 and 45° N is weaker in NorESM2 compared to EC-Earth3. Likewise, the upper-level potential temperature increase between 45 and 70° N is weaker, resulting in a comparable baroclinicity response. NorESM2, like the other models, shows the strongest increase in MFC co-located with the zonal wind speed increase (Fig. 11e). Moreover, like EC-Earth3, NorESM2 exhibits the compensating effect of baroclinicity and MFC, found at 60 to 70° N and 400 to 200 hPa, resulting in a non-significant change in zonal wind speed (Fig. 8e).

Lastly, CESM2 shows a baroclinicity response structurally similar to NorESM2, but with an amplified magnitude (Fig. 9g). This results from weaker warming of the atmospheric column at higher latitudes (55–70° N) in CESM2 compared to NorESM2, while the warming at lower latitudes (30–55° N) is comparable (Fig. 10g). Unlike the other models, the change in baroclinicity is largest south of the maximum in the baseline simulation. Together with the change in MFC, this increases the jet speed located southward of the maximum in the baseline simulation. Furthermore, MFC contributes a weak amplification of the zonal wind speed at lower levels (850–600 hPa, 52–58° N) in Fig. 11g.

5.2 The mechanisms for SIC response

OpenIFS-43r3 shows the geographically most extensive reduction in baroclinicity across all models (Fig. 9b), with large decreases spanning from 47 to 70° N and from 850 to 300 hPa. The baroclinicity response originates from deep vertical warming at 60 to 70° N (Fig. 10b) with the strongest warming located at the surface. This warming is directly induced by sea ice removal, exposing the sea surface to the atmosphere. A secondary effect is a cooling at lower latitudes (35–50° N), which has its maximum at 500 to 300 hPa. This dipole structure has been identified in previous studies (Screen et al., 2022; Labe et al., 2020). Furthermore, MFC shows a non-significant reduction co-located with its baseline maximum (Fig. 11b). However, the change in baroclinicity is the primary driver of the zonal wind speed response to sea ice loss, and MFC plays a secondary role.

The response of the zonal wind speed to changes in sea ice in EC-Earth3 is not significant. This is the result of opposing weak effects of MFC and B_y. EC-Earth3's baroclinicity shows the weakest ΔSIC of all models (Fig. 9d), which is due to the weaker vertical and broader horizontal extent of warming due to sea ice loss. Relevant to the NA jet stream exit is the reduction in B_y found at 50 to 60° N and 850 to 600 hPa, above which a consistent minor increase in MFC is present. A similar structure is found at 65 to 70° N, where a MFC decrease is associated with an increase in B_y. The dominant factor in NorESM2 is MFC (Fig. 11f). There is a dipole with an increase on the equatorward side and a decrease on the poleward side of the baseline maximum. This results in the equatorward shift of the NA jet stream exit in NorESM2. Supporting the equatorward shift in zonal wind speed is a decrease in baroclinicity on the poleward side (60–70° N, 850–700 hPa) and an increase on the equatorward side (40–50° N, 500–300 hPa), shown in Fig. 9f. The cause of the reduction in baroclinicity is shallow warming at high latitudes (60–70° N, 850–700 hPa), resulting directly from the removal of sea ice (Fig. 10f). Additionally, the increase in B_y results from a warming at 30 to 40° N and a cooling at 50 to 65° N.

CESM2 responds to the decreased SIC with a vertically extensive reduction in B_y, reaching from 850 hPa to 400 hPa located between 53 and 70° N (Fig. 9h). This reduction originates similarly to OpenIFS-43r3 from deep high-latitude vertical warming (60–70° N, 850–400 hPa) and an adjacent cooling equatorward of the warming (Fig. 10h). The MFC presents a dipole structure, which is not associated with a corresponding dipole in the zonal wind speed, indicating the baroclinicity change is the dominant factor for the SIC response in CESM2.

6 The extratropical cyclone response

6.1 The future response of ETC track density

Figure 12 shows the future response of the ETC densities for all four models. In OpenIFS-43r3, the majority of the ETC track density changes compared to the baseline simulation, including those over Europe, are non-significant. In EC-Earth3, there is an increase in ETC track density exceeding 4 ETCs per 5° spherical cap per winter season downstream of the baseline storm track (Fig. 12b). This is tied to the eastward extension of the NA jet stream and is located northward of the zonal wind speed increase. Additionally, there is a significant decrease east of Greenland and over the European continent. NorESM2 is set apart by a significant decrease in ETC track density in the storm track in the future simulation (Fig. 12c). The maximum decrease is found around Iceland with a magnitude > 12 ETCs per 5° spherical cap per winter season. Similarly to EC-Earth3, the increase in ETC track density over northern Europe is tied to the increase in the NA jet stream. CESM2 presents an increase downstream of the baseline storm tracks (Fig. 12d). As with other models, this is co-located with the increase the NA jet stream speed in the future simulation. Also, CESM2 shows its strongest reduction in ETC track density eastward of Greenland and over the western Atlantic. All models show a decrease in ETC density over the Mediterranean; however, this is only significant in three out of four models, with OpenIFS-43r3 showing non-significant changes.

Figure 12 Future response ΔFT of the 39-year DJF mean ETC track density (in colour shading) and the baseline climatology (in black contours) for (a) OpenIFS-43r3, (b) EC-Earth3, (c) NorEMS2, and (d) CESM2. The unit is number of ETCs per 5° spherical cap per winter season (DJF). Shading indicates a consistent sign according to CDR testing, and stippling indicates statistically insignificant changes (FDR-corrected t test).

6.2 The SST and SIC contributions

Due to the sparser nature and higher variability of the ETC dataset compared to, for example, the zonal wind speed, the contribution of SSTs and SIC do not sum up to the future response, resulting in a larger residual. Nonetheless, it is useful to look at the effects of changed SST and SIC individually.

Figure 13 SST response ΔSST (a, c, e, g) and SIC response ΔSIC (b, d, f, h) of the 39-year DJF mean ETC track density (in colour shading) and the baseline climatology (in black contours) for (a, b) OpenIFS-43r3, (c, d) EC-Earth3, (e, f) NorEMS2, and (g, h) CESM2. The unit is the number of ETCs per 5° spherical cap per winter season (DJF). Shading indicates a consistent sign according to CDR testing, and stippling indicates statistically insignificant changes (FDR-corrected t test).

OpenIFS-43r3 shows a reduction of ETC track density over the eastern USA and western Atlantic as a response to increased SSTs (Fig. 13a). There is an increase in ETC in the eastern Atlantic (48–60° N, 40–5° W), north of the zonal wind increase (Fig. 7a), which extends the storm track to the east. Meanwhile, in the SIC response (Fig. 13b), a clear feature is a dipole southeast of Greenland, effectively leading the southward shift. In the Mediterranean, ΔSST and ΔSIC oppose each other, with ΔSST showing a decrease in ETC track density and ΔSIC showing an increase.

EC-Earth3's ΔSST is characterised by a southeastward shift of the storm track (Fig. 13c) with the strongest decrease in ETC densities between Greenland and Iceland. Moreover, there is a significant decrease in ETC track density in the Mediterranean area and the Arctic Ocean and an extensive increase over the eastern Atlantic Ocean (40–58° N, 40–5° E). The SIC response presents a similar dipole (Fig. 13d), with a reduction southeast of Greenland and an increase over the central Atlantic (40–58° N, 50–20° W). The compounding effects of ΔSST and ΔSIC explain the clear increase over northern Europe and decrease over central Europe in the future response in EC-Earth3. Both ΔSST and ΔSIC show a strong connection with the jet stream changes (Fig. 7c, d). The widespread increases in ETC track density are located northward of the increase in jet speed.

NorESM2 shows significant reductions in ETC track density from the eastern USA along the western Atlantic toward Greenland and in the Mediterranean in response to increased SSTs (Fig. 13e). In these regions, ΔSST is the dominant contribution to ΔFT. The change in SIC (Fig. 13f) causes a reduction over Iceland and the Arctic Ocean and an increase over northern Europe and the British Isles. Hence, both SSTs and SIC contribute to the increase over northern Europe and the British Isles, but ΔSIC has a larger magnitude (8 to 9 ETCs per 5° spherical cap per winter season).

In CESM2, similar to NorESM2, the changed SSTs reduce the ETC track density along the baseline simulation storm track and the Mediterranean (Fig. 13g). However, the reduction in the Arctic Ocean due to SST changes is more pronounced in CESM2. The increases in ETC track density in the North Atlantic are statistically insignificant, and like the other models, they are located northward of the increases in the zonal wind. The SIC response has the largest impact in the Arctic Ocean (Fig. 13h, 60–80° N, 25° W–20° E). The future response of CESM2 shows a clearer signal compared to the individual contributions of SST and SIC. However, the individual contributions have similar spatial patterns, indicating that the future response in CESM2 is not dominated by either SST or SIC.

6.3 The impacts on ETC properties

Figure 14 presents the changes in the mean value of ETC count, ETC lifetime, mean speed, maximum intensity, genesis latitude, and latitudinal displacement for ETC tracks affecting Europe (at least 48 h within 30–70° N, 15° W–35° E). Similar to the track densities, ΔSIC and ΔSST do not necessarily add linearly to exactly equal the future response. However, they provide insight into the dominant contribution to the future response and help to compare the responses across models. The diamonds in Fig. 14 represent the mean value of the ETC properties in the baseline simulation, previously discussed in Sect. 3.3 (Fig. 4). The inclusion of the baseline simulation sets the responses of each model in the context of model differences. Considering ETC properties, the majority of future, SIC, and SST responses are smaller than the differences between the individual models. In particular, ETC count and genesis latitude are properties which show responses of a magnitude similar to the model differences.

Figure 14 The change in mean values of ETC count per winter season (denoted DJF⁻¹) (a), ETC lifetime (b), mean speed (c), maximum intensity (d), genesis latitude (e), and latitudinal displacement (f) for OpenIFS-43r3 (red), EC-Earth3 (blue), NorESM2 (green), and CESM2 (purple) with saturated colours indicating statistical significance (Mann–Whitney U test with p value < 0.05). Diamonds represent the baseline simulation, and crosses represent the corresponding experiment. Arrows are added to aid the interpretation of the direction from baseline to the corresponding experiment.

The ETC count per winter season (denoted DJF⁻¹) shows a decrease in the future response originating from ΔSST (Fig. 14a) across models, with the change being statistically significant in NorESM2 (ΔFT, −2.3 DJF⁻¹) and CESM2 (ΔFT, −4.5 DJF⁻¹). Models disagree on the effect of changed SIC on ETC count. However, OpenIFS-43r3 has the largest change amongst models by +2.4 per winter season (DJF), which corresponds to the increase in ETC track density in the Mediterranean (Fig. 13b).

The SIC and SST responses of ETC lifetime contribute similarly in magnitude to the future response (Fig. 14b). EC-Earth3 and CESM2 present an increase in lifetime for ΔSST and ΔSIC, which results in a significant increase in the future response (EC-Earth3: +5.0 h, CESM2: +7.5 h). Meanwhile, OpenIFS-43r3 shows a non-significant decrease in ETC lifetime for ΔSST and an increase for ΔSIC, which compensate each other in the future response. There are only minimal, non-significant changes to the ETC lifetime in NorESM2. The mean speed of ETCs does not show any significant changes (Fig. 14c), besides a deceleration in EC-Earth3's future response (−0.27 m s⁻¹). Additionally, the models disagree on the sign of the change for the future, SST, and SIC responses.

For the mean maximum intensity (quantified by the magnitude of the pressure anomaly), models agree on an increase in response to the SST changes, indicating a deepening of ETCs, and a decrease in response to the SIC changes, corresponding to a reduction of the pressure minimum (Fig. 14d). However, many of these changes are not significant. NorESM2 shows a significant decrease due to SIC changes (−1.8 hPa) and a minimal non-significant change due to ΔSST, which results in a reduced maximum intensity on average in the future response. Next, while EC-Earth3 presents a significant increase in maximum intensity for ΔSST, it shows a non-significant decrease in ΔFT and ΔSIC.

The future response of the genesis latitude is dominated by the change in SSTs with a significant increase across all models (Fig. 14e). Averaged across models, the genesis latitude increases by +1.1° N in ΔFT and +1.3° N in ΔSST. The SIC responses are non-significant across models. Figure 14f shows the latitudinal displacement of ETCs, which is significantly reduced in OpenIFS-43r3, EC-Earth3, and NorESM2 for the future response. For NorESM2, the change originates from the SST change and the interaction of the SST and SIC changes. Meanwhile, ΔSIC dominates for OpenIFS-43r3. The SST responses are not significant in any of the four models. Combining the changes in the genesis latitude and latitudinal displacement, ETCs are propagating further north in the SST change only simulation. However, in the future climate scenario, ETCs travel more zonally in OpenIFS-43r3, EC-Earth3, NorESM2, which indicates important interactions of SST and SIC changes.

7 Discussion

Multiple aspects of the results require a more detailed discussion. Foremost, it is important to consider the robustness of the results obtained here given that a previous study (Peings et al., 2021) suggested that longer simulations, or larger ensembles, are required to isolate consistent atmospheric responses to sea ice loss. The atmospheric response due to sea ice loss found in this study is weak in the mid-latitudes. This may be partially caused by the models underestimating the responses due to eddy feedback that is too weak, as shown previously by Smith et al. (2022) and Screen et al. (2022). Nonetheless, consistent signals in the sign of the response are identified using the CDR, leading to confidence in the sign of the response if not the magnitude of the response. In comparison to previous studies, the sea ice loss forcing applied here is stronger, making it easier to distinguish the response from internal variability. Furthermore, the SIC response combined with the responses due to changing SST are coherent with the future responses, where both SST and SIC are changed. While the results should not be overinterpreted, they provide evidence for the physical mechanisms contributing to the future change to the North Atlantic jet stream and storm track in addition to the largely descriptive results of how the jet and storm track change.

The future winter North Atlantic jet stream (identified at 250 hPa, Fig. 5) intensifying on its southern flank and extending further into Europe in EC-Earth3, CESM2, and NorESM2 is consistent with findings across CMIP3, CMIP5, and CMIP6 by Harvey et al. (2020), albeit for a weaker warming scenario (SSP 2-4.5) compared to this study (SSP 5-8.5). The remaining future response of OpenIFS-43r3 is a scientifically important outlier, particularly as it is the same atmospheric model component used in EC-Earth3, albeit a much newer version. The SST response is structurally consistent across all models (Fig. 7). The reduction of the jet speed induced by sea ice loss at 250 hPa over Europe in OpenIFS-43r3, CESM2, and NorESM2 is present in the PAMIP simulation (Ye et al., 2023). However, the spatial overlap of the SST and SIC responses in OpenIFS-43r3 is the important difference with the other models, which explains the lack of an extension of the jet stream over Europe. Using greenhouse gas and sea ice loss forcing, Oudar et al. (2020) have shown a similar compensating effect on zonal wind speed at 850 hPa over Europe.

There are striking differences in the vertical extent of the warming induced by sea ice loss at high latitudes across models (Fig. 10). The vertical extent of high-latitude warming has been shown to be important for mid-latitude atmospheric circulation. For example, the vertical extent of high-latitude warming has been investigated as a contributor to the Eurasian cooling pattern in future climate projections (He et al., 2020; Labe et al., 2020). Moreover, idealised simulations suggest an enhanced reduction in the zonal wind speed with high-latitude warming aloft compared to the same magnitude of warming near the surface (Kim et al., 2021). Results by Xu et al. (2023) indicate that sea ice loss in the Barents Sea and Kara Sea contributes to a stronger high-latitude warming aloft, which in turn impacts mid-latitude circulation. Similarly to Labe et al. (2020), this study shows enhanced upper-tropospheric cooling in the mid-latitudes with higher vertical extent of warming at high latitudes in OpenIFS-43r3 and in weaker form in CESM2 (Fig. 10b, h). This is associated with a decrease in baroclinicity in the mid-troposphere (Fig. 9b, h), which impacts the jet speed over Europe (Fig. 8b, h). In contrast, EC-Earth3 shows a moderate vertical extent of warming without an associated upper-tropospheric cooling in the mid-latitudes, which results in weak changes in baroclinicity and no significant impact on zonal wind speed. Lastly, despite a shallow vertical extent in high-latitude warming, NorESM2 shows a weak deceleration of zonal wind speed, which is explained by the changes driven by changes in the momentum flux convergence rather than in baroclinicity.

The future response of the storm tracks found in this study closely resemble results from CMIP6 (Priestley and Catto, 2022), where it was shown that under the SSP 5-8.5 scenario the ETC track density increases over the British Isles and northern Europe, while it decreases in the Mediterranean. Previous studies showed that the North Atlantic storm track shifts south in response to sea ice loss (Hay et al., 2023; Yu et al., 2023). This study supports these finding. However, this southward shift is model dependent, with EC-Earth3 and NorESM2 showing strong responses while CESM2 and OpenIFS-43r3 are comparatively weak (Fig. 13). In contrast to Yu et al. (2023), ΔSST and ΔSIC have contributions of similar magnitude to the future response of the North Atlantic storm track for each individual model. It is important to note that the storm tracks exhibit high variability; therefore, the precise decomposition of the contributions of SST and SIC, including the non-linear interactions between SST and SIC, requires further investigation with statistical methods like principal component analysis.

Furthermore, the sea ice response results of the ETC properties agree with Hay et al. (2023), who show that ETCs tend to have a longer lifetime and a reduced intensity. Similarly, the number of ETCs reduces in most models (EC-Earth3, CESM2, NorESM2), while it increases in others (OpenIFS-43r3). In agreement with CMIP6 results, the future responses of all four models agree on a reduction in the number of ETCs (Priestley et al., 2023) and a more zonal propagation (Crawford et al., 2023).

8 Conclusions

This study presents the responses of the North Atlantic jet stream and storm track to future sea surface temperatures and sea ice cover and quantifies their contributions using the set of atmosphere-only model simulations from the CRiceS project (CRiceS, 2021). Furthermore, the potential driving mechanism of the responses were examined by using the baroclinicity and momentum flux convergence. The key findings of this study are as follows:

• The baseline simulations of the four models agree on the general shape of the North Atlantic jet stream, which is in reasonable agreement with results from the ERA5 reanalysis and CMIP6 (Harvey et al., 2020). A key difference is that the zonal wind speed (u) in the upper troposphere is higher over Europe in NorESM2 and CESM2 compared to EC-Earth3 and OpenIFS-43r3. While the baroclinicity is comparable across models, the difference in the jet strength over Europe is primarily explained by the magnitude of momentum flux convergence in the baseline simulations, which is larger in NorESM2 and CESM2. The ETC track density is similar across models, and the patterns of the ETC track density in the baseline simulation closely match previous findings based on CMIP6 simulations (Priestley et al., 2020), as well as PAMIP (Yu et al., 2023).

• In the future response, three out of the four models show an increase in zonal wind speed at 250 hPa on the southeastern side of the baseline jet core over the eastern North Atlantic and reaching far into Europe, which has been shown in CMIP6 simulations (Harvey et al., 2020). The ΔSST contribution dominates the future response by an increase in the upper tropospheric baroclinicity. The momentum flux convergence contributes to an equatorward shift of the NA jet stream in the SIC response.

• The exception is OpenIFS-43r3 with no significant change in zonal wind speed over Europe. The SST response is similar to the other three models; meanwhile, the SIC response shows a decrease in zonal wind speed over Europe, which is driven by a decrease of more than 10 % in baroclinicity at lower levels (850–400 hPa). These overlapping opposite responses result in no significant changes to the jet stream over Europe in the future scenario.

• The ETC track density reduces along the east coast of North America and Greenland in the future scenario. SST and SIC changes both significantly contribute to the reduction. A consistent feature is a significant reduction of ETC track density in the Mediterranean basin, which is driven by the changes in SSTs. The increases in ETC track density over the eastern North Atlantic, leading to an eastward extension of the storm track, are associated with the extension of the NA jet stream. This eastward extension (and reduction) in the Mediterranean basin is also present in CMIP projections (Priestley and Catto, 2022). Similar to the jet stream changes, OpenIFS-43r3 shows no significant increases over Europe for the ETC track density in the future scenario.

• The future changes to the mean values of the ETC-specific properties are superseded by the inter-model differences. For example, mean speed and maximum intensity show no consistent statistically significant changes for two or more models. However, ETCs affecting Europe tend to become less frequent, which is linked to the reduced ETC track density in the Mediterranean. They also originate more poleward in the future scenario, which is driven by changes in SSTs, and travel more zonally, driven by the interaction of SST and SIC changes.

• The total future response of the North Atlantic at 250 hPa is dominated by the changes in sea surface temperature. On the other hand, the contributions of changes in sea ice cover and sea surface temperature to the total future response of the North Atlantic storm track are of similar magnitude.

This study is based on atmosphere-only simulations with annually repeating sea ice cover and sea surface temperatures, which limits the direct comparability between the studied simulations and results from reanalysis and CMIP. Furthermore, the changes presented are the mean responses of the jet stream and storm track, leaving the investigation of changes to the variability and extreme events of the jet stream and storm track to future research. Moreover, one aim of the coordinated simulations is to constrain the variability of the climate system. However, the internal atmospheric climate variability is a major source of uncertainty in this study. This affects both the study between the differences between simulations within one model and the comparison between the four different models. In particular, the metrics related to extratropical cyclones are strongly affected by atmospheric internal variability. In addition, motivated by the model-specific physical interpretability, the number of models is limited. Thus, it is difficult to assess whether the disagreement of the OpenIFS-43r3 response is an outlier or representative of a different part of the spread of state-of-the-art models.

This study encourages further investigation into the interannual variability of the baroclinicity and eddy-driven mechanisms of the jet stream and how this connects to the storm track activity. Moreover, as this study investigates the mean changes to extratropical cyclones, further work on extreme cyclones and associated extreme events of wind and precipitation is necessary.

Notable differences in the future responses occur, despite the models having similar baseline climates and being forced with the same changes to SSTs and sea ice cover. Furthermore, although the SST response dominates, the sea ice cover (SIC) response is significant in some areas. For example, in OpenIFS-43r3 the SIC response is large enough to counteract the impact of changing SSTs. Overall, substantial uncertainties remain in how the jet stream and storm track in the North Atlantic and Europe will change in the future.

Appendix A

Table A1 Additional information for Fig. 4. For the six ETC quantities, the mean value for each model is reported on the diagonal (in parentheses), and the differences between the models (row – column) are given in the off-diagonal cells. A Mann–Whitney U test with p value < 0.05 is highlighted in bold font.

Figure A1 DJF mean of the boundary conditions of sea surface temperature (SST) (shading) and sea ice cover (SIC) (white where SIC exceeds 0.5) for the baseline simulations (a) and future SSP 5-8.5 simulation (b).

Figure A2 DJF mean of zonal wind speed u in ERA5 reanalysis (years: 1950–1969). Shown are the 250 hPa level in colour shading and the 850 hPa level in black contours (m s⁻¹).