Future changes in compound explosive cyclones and atmospheric rivers in the North Atlantic

The explosive development of extratropical cyclones and atmospheric rivers plays a crucial role in driving extreme weather in the mid-latitudes, such as compound windstorm–flood events. Although both explosive cyclones and atmospheric rivers are well understood and their relationship has been studied previously, there is still a gap in our understanding of how a warmer climate may affect their concurrence. Here, we focus on evaluating the current climatology and assessing changes in the future concurrence between atmospheric rivers and explosive cyclones in the North Atlantic. To accomplish this, we independently detect and track atmospheric rivers and extratropical cyclones and study their concurrence in both ERA5 reanalysis and CMIP6 historical and future climate simulations. In agreement with the literature, atmospheric rivers are more often detected in the vicinity of explosive cyclones than non-explosive cyclones in all datasets, and the atmospheric river intensity increases in all the future scenarios analysed. Furthermore, we find that explosive cyclones associated with atmospheric rivers tend to be longer lasting and deeper than those without. Notably, we identify a significant and systematic future increase in the cyclones and atmospheric river concurrences. Finally, under the high-emission scenario, the explosive cyclone and atmospheric river concurrences show an increase and model agreement over western Europe. As such, our work provides a novel statistical relation between explosive cyclones and atmospheric rivers in CMIP6 climate projections and a characterization of their joint changes in intensity and location.

1 Introduction

Atmospheric rivers (ARs) are narrow and elongated corridors of horizontal moisture transport usually associated with the cold front of an extratropical cyclone (Bao et al., 2006; Ralph et al., 2004, 2017). They play an essential role in the atmospheric water vapour cycle in the mid-latitudes, accounting for 90 % of the poleward moisture transport (Zhu and Newell, 1998; Guan and Waliser, 2015). Moreover, ARs drive wet and windy extreme weather events, particularly in western continental coasts, such as western Europe (Zhu and Newell, 1998; Guan and Waliser, 2015; Lavers and Villarini, 2015, 2013; Gimeno et al., 2016). In some coastal areas, intense ARs are associated up to 95 % of the time with extreme precipitation and up to 75 % with wind extremes (Waliser and Guan, 2017). As a result, ARs have widespread socioeconomic impacts; for example, the majority of European storms causing insurance losses of a billion US dollars or more are linked to ARs (de Vries, 2021).

Another meteorological feature that can lead to wet and windy extremes is explosive cyclones (ECs), which are also known as weather bombs (Roebber, 1984; Reale et al., 2019). ECs are rapidly intensifying mid-latitude cyclones. Historically, ECs have been identified as those with a deepening rate of more than 24 hPa in 24 h and scaled by the latitude (Sanders and Gyakum, 1980), although this definition has been repeatedly challenged, especially for the Southern Hemisphere (Allen et al., 2010). Events such as the Presidents' Day Snowstorm of 1979 (Schultz, 2022), characterized by poor forecast accuracy, sparked extensive research into this phenomenon. Upon landfall during their intensification phase, ECs can produce widespread damage and impacts associated with strong winds, heavy precipitation, and storm surges (Fink et al., 2012; Liberato et al., 2013; Ludwig et al., 2015; Seiler and Zwiers, 2016b; Reale et al., 2019; Ginesta et al., 2023).

ECs and the associated ARs thus play a crucial role in driving extreme weather hazards in the mid-latitudes (Liberato et al., 2013; Davolio et al., 2020). An illustrative example is Storm Alex, an EC associated with an AR that at first produced extreme winds in France and the UK and then led to record-breaking precipitation in Italy in October 2020 (Davolio et al., 2023; Ginesta et al., 2023). The climatological relationship between ECs and ARs has been previously studied and the literature shows that ARs are more often found in the surroundings of EC than non-ECs (Eiras-Barca et al., 2018; Zhang et al., 2019; Guo et al., 2020). ARs are important sources of moisture for cyclonic systems, and it has been suggested that they can enhance cyclone deepening through moist diabatic processes (Zhu and Newell, 1994; Ferreira et al., 2016), such as cloud condensation (Pinto et al., 2009). In addition, ECs with ARs show larger moisture inflow and deepen more rapidly than ECs without an AR but do not show significant differences in low-level baroclinicity nor upper-level potential vorticity, suggesting that diabatic processes are important contributors to their intensification (Zhang and Ralph, 2021). However, the extent to which these moist diabatic processes, compared to other factors such as upper-level forcing, influence cyclone intensification can vary from case to case (Pfahl and Sprenger, 2016; Ginesta et al., 2024).

A range of studies have investigated the impact of anthropogenic climate change on extratropical cyclones and ARs individually (Lavers et al., 2015; Zappa et al., 2013). The thermodynamic response of ARs to climate change is characterized by an increase in integrated water vapour transport (IVT). This increase is driven by the Clausius–Clapeyron relation, which implies a rise in moisture content in a warmer atmosphere. However, integrated water vapour (IWV) is expected to experience a larger increase than surface water vapour under climate change (Payne et al., 2020). This thermodynamic signal would act to increase the number of ARs detected in a warmer climate (Thandlam et al., 2022; Espinoza et al., 2018; Zhang et al., 2024; O’Brien et al., 2022; Wang et al., 2023). Similarly, the thermodynamic response acts to increase the precipitation within extratropical cyclones (Yettella and Kay, 2017). The dynamical response to climate change, such as changes in atmospheric circulation patterns, is more uncertain (Shepherd, 2014). In the North Atlantic, dynamic changes are mostly driven by changes in the eddy-driven jet, which serves as a guide for extratropical cyclones. The tug of war between the upper tropospheric warming and the Arctic amplification leads to a high uncertainty regarding the changes in the jet over the North Atlantic and western Europe (Shaw et al., 2016). However, climate models indicate a decline in the number of extreme extratropical cyclones in the North Atlantic and a local increase over the North Sea in EC frequency (Priestley and Catto, 2022; Zappa et al., 2013; Seiler and Zwiers, 2016a). Regarding ARs, studies also point to an increase in frequency, intensity, and size in western Europe and a northward shift in the AR location and landfall (Lavers et al., 2013; Ramos et al., 2016; Gao et al., 2016; Zhang et al., 2024).

The hazards associated with the joint occurrence of explosive cyclones and atmospheric rivers, especially along the western coast of Europe, underscores the importance of evaluating the projected changes in their concurrence in future climates. In this study, we assess future projections of the interplay between ECs and ARs using state-of-the-art models from the Coupled Model Intercomparison Project phase 6 (CMIP6; Eyring et al., 2016b). We specifically evaluate the frequency of EC and AR concurrence in the present climate for the ERA5 reanalysis and CMIP6 models, and in three end-of-century scenarios for the CMIP6 models. Moreover, we also assess the future changes in the intensity and location of such compound events.

This paper is structured as follows. Section 2 describes the datasets employed in this study, while Sect. 3 explains the methodologies for tracking ARs and cyclones and the calculation of their concurrences. In Sect. 4 we evaluate and discuss the performance of the CMIP6 models with ERA5 reanalysis. Section 5 shows the results and discussion of the future changes in frequency, intensity, and location of this compound event. Finally, Sect. 6 summarizes the key findings and provides the conclusions of this study.

Figure 1 EC track density climatology during the extended winter (ONDJFM) (1980–2009) for ERA5 (a), historical simulations of CMIP6 models (b), and the difference between ERA5 and CMIP6 (c). Panels (d, e, f) are the same as (a, b, c) but for non-ECs. The values are given in units representing the number of cyclone time steps per 3° spherical cap per month.

2 Data

We use the ECMWF reanalysis ERA5 (Hersbach et al., 2020) with a horizontal resolution of 0.25° × 0.25° as an observationally constrained reference for the current climate and validation of the global climate models (GCMs). For the AR detection, the variables used are specific humidity and meridional and zonal wind components at 1000, 925, 850, 700, 600, 500, 400, and 300 hPa (Sect. 3.2). For cyclone detection, the variable used is the sea level pressure (SLP). In both cases, the variables are at 6-hourly resolution during the extended winter period (October to March) between 1980 and 2009 for the North Atlantic region (15–75° N, 90° W–20° E) (hereafter referred to as the “data domain”). To mitigate issues caused by the cyclone tracking algorithm, which tends to generate stationary “artefact” cyclones along the western boundary, we apply a buffer zone of 10° to this data domain. This adjustment ensures that artefacts at the western boundaries are excluded. As a result, the domain used for all analyses presented in the results section is (25–65° N, 80° W–10° E), referred to simply as the “domain” from here on.

Further, we use one member each of six different GCMs from the CMIP6 dataset (Eyring et al., 2016a): MPI-ESM1-2-LR, MPI-ESM1-2-HR, NorESM2-MM, EC-Earth3, CMCC-ESM2, and MIROC6. More detailed information about the GCMs used can be found in the Appendix A1. We evaluate the listed GCMs for two periods: the current climate (1980–2009), using historical simulations, and the future climate at the end of the 21st century (2070–2099), where we used simulations following the SSP1-2.6, SSP2-4.5, and SSP5-8.5 forcing scenarios (Riahi et al., 2017). These 30-year datasets consist of 29 full winters and two partial winters (January to March for 1980 and 2070 and October to December for 2009 and 2099). Variables used for tracking ARs and cyclones in the GCMs are the same as used for ERA5. The current list of GCMs used in this study is limited by the availability of 6-hourly instantaneous variables for the historical and the three scenarios experiments in hybrid-sigma pressure model levels, which are required to interpolate to the necessary pressure levels for the IVT calculation when detecting ARs. This limitation in the analysed data prevents a complete assessment of model performance, as only a single member from each model is used, and a multi-member ensemble for each model would be necessary for a more robust evaluation of the individual model uncertainty. For this reason, we evaluate the results using the multi-model mean of the ensemble.

Figure 2 AR frequency climatology during the extended winter (ONDJFM) (1980–2009) for ERA5 (a), historical simulations of CMIP6 models (b), and the difference between ERA5 and CMIP6 (c). The values are given in units representing the percentage of time steps detected as AR.

3 Methods

3.1 Extratropical cyclone tracking

There are several detection and tracking algorithms for extratropical cyclones, most of them using either SLP or lower-tropospheric vorticity (Neu et al., 2013). Here, we detect and track cyclones based on SLP using the TempestExtremes algorithm developed by Ullrich et al. (2021); Ullrich (2020). This command-line software facilitates adaptable and rapid feature detection and tracking for extratropical cyclones and for ARs.

To identify extratropical cyclones, we recognize candidate “nodes” corresponding to local minima in the SLP field with the same set-up as in Ullrich et al. (2021). Nodes within a 6° great circle distance (GCD) of each other are merged. Next, to connect these candidate nodes into tracks, the distance between consecutive detections should not exceed 6° GCD. The tracks must persist for a minimum of 24 h, the maximum duration between two detections is set at 6 h, and the cyclones must have moved at least 12° GCD to filter out stationary lows, such as the Icelandic Low. The relevant code and the detailed set-up are available in the Appendix A2.

In addition, extratropical cyclones are classified as ECs if their normalized deepening rate (NDR), as defined by Sanders and Gyakum (1980), is equal to or higher than 1:

$$\begin{array}{}\text{(1)}& {\mathrm{NDR}}_c={\displaystyle \frac{{\mathrm{DR}}_{\mathrm{24}\mathrm{h}}}{\mathrm{24}\mathrm{hPa}}}{\displaystyle \frac{\sin \left(\mathrm{60}\mathit{°}\right)}{\sin \left(\mathit{\phi }\right)}},\end{array}$$

where DR_24 h is the pressure difference over 24 h measured at the storm centre and φ is the latitude at its second time step. Cyclones that do not fulfil this condition are classified as non-ECs. Table A2 in Appendix A4 lists the number of cyclone tracks of ECs and non-ECs in each model and scenario and in ERA5 within the North Atlantic domain. The number of cyclone tracks detected fits inside the range number of cyclones detected in Neu et al. (2013) using several cyclone tracking algorithms.

Figure 1 shows the number of cyclone time steps detected within a 3° spherical cap per month, hereafter referred to as cyclone track density, for ECs and non-ECs in both ERA5 and CMIP6 models. These results agree with Priestley and Catto (2022) for CMIP6 and Zappa et al. (2013) for CMIP5 track densities despite using different tracking algorithms. Our tracking method shows some differences between ERA5 and CMIP6 models (Fig. 1c, f). CMIP6 models show smaller values for cyclone track density along the North Atlantic storm track, particularly east of Newfoundland, south of Greenland, and in the North Sea. For non-ECs, CMIP6 show a northward shift in the storm track with lower cyclone track density in the south and higher density in the north of the domain.

3.2 AR detection and tracking

There are two primary approaches to detecting ARs: one involves using the integrated water vapour, commonly applied to satellite data, while the other more broadly used approach consists of computing the IVT (Gimeno et al., 2014; Shields et al., 2018). In this study, we calculate the IVT for both reanalysis and model data. The IVT is defined for each grid point as follows:

$$\begin{array}{}\text{(2)}& \mathrm{IVT}={\left[{\left({\displaystyle \frac{\mathrm{1}}{g}}\underset{\mathrm{1000}\mathrm{hPa}}{\overset{\mathrm{300}\mathrm{hPa}}{\int }}qudp\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{1}}{g}}\underset{\mathrm{1000}\mathrm{hPa}}{\overset{\mathrm{300}\mathrm{hPa}}{\int }}qvdp\right)}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where q is the specific humidity, u is the zonal wind component, v is the meridional wind component, and g is the gravitational acceleration. Moreover, we separately compute the eastward (IVT_E) and northward (IVT_N) components of IVT, which correspond to the two terms inside of the brackets in Eq. (2), respectively.

For the detection of ARs, we also use the TempestExtremes algorithm developed by Ullrich et al. (2021); Ullrich (2020). We find candidates for atmospheric rivers by detecting ridges in the IVT field. Ridges are defined as points where the Laplacian of the IVT is below $-\mathrm{4}\times {\mathrm{10}}^{\mathrm{4}}$ kg m⁻² s⁻¹ rad⁻², as this operator identifies elongated areas and regions of local minima. In addition, the IVT should be higher than 250 kg m⁻¹ s⁻¹. Each candidate should have an area larger than 4×10⁵ km⁻². These thresholds were defined and tested by Ullrich et al. (2021) and show strong agreement with other tracking algorithms (Collow et al., 2022). The detected candidates are concatenated if at least one grid point is identified as an AR in consecutive time steps, meaning that the AR area at a consecutive time step spatially overlaps with the previous AR. In addition, they should last 60 h. The successful ARs that fulfil all these requirements are labelled as AR grid points and will be used later to determine if cyclones are concurrent with an AR (Sect. 3.3). This AR tracking methodology is less sensitive to a generalized increase in IVT in the future climate due to the Clausius–Clapeyron relationship because we detect AR candidates using the Laplacian of the IVT instead of the IVT field. Thus, AR candidates are detected when having a pronounced gradient of IVT (not dependent on the background IVT), as discussed further in Sect. 5. This detection algorithm may influence the observed changes in the number of AR tracks in future climates. The relevant code is available in Appendix A3, and the detailed number of AR tracks detected is reported in Table A3 in Appendix A4. Figure 2a shows the percentage of time steps with a detected AR for ERA5, hereafter referred to as AR frequency. Our tracking methodology accurately reproduces the AR frequency when compared with Guan and Waliser (2015). The AR frequency in the historical simulations of CMIP6 models is higher compared to ERA5 and exhibits a southward shift, with more AR time steps detected primarily in the lower mid-latitudes (Fig. 2b, c).

3.3 Concurrences

For each extratropical cyclone (both ECs and non-ECs) we compute the maximum deepening point (MDP), which is the maximum difference in SLP between two consecutive 6-hourly time steps. This metric allows us to evaluate the influence of ARs on the development of the cyclone before and after its maximum intensification, which is when the potential impact of ARs on cyclone deepening is expected, following the approach used by Eiras-Barca et al. (2018).

Subsequently, we determine whether a specific time step of an extratropical cyclone (EC or non-EC) is linked to an AR if at least one grid point within 1500 km from the centre of the cyclone is part of an AR (detected and tracked independently; see Sect. 3.2). This AR search radius is consistently applied across all time steps of the EC and non-EC tracks. By selecting a 1500 km radius, our methods align with those of (Eiras-Barca et al., 2018), with the primary difference between the two methods being the AR and cyclone tracking algorithms used. Almost 50 % of the identified ARs are located in the southeastern quadrant of the cyclone (Supplement Figs. S1 and S2). The rationale for using a 1500 km radius is based on the consideration that an AR can influence the cyclone by delivering moisture to the warm conveyor belt or feeder airstream (Dacre et al., 2019), which are typically located in the southeastern quadrant and within this distance from the cyclone centre.

Figure 3 shows an example of our detection methodology applied to storm Xynthia. EC Xynthia underwent rapid intensification before making landfall on 27 February 2010. It caused widespread damage across western European countries, especially France and Spain. In Fig. 3, the black dots represent the cyclone's path, while the red crosses indicate the time steps concurrent with the presence of an AR. The shaded areas depict the regions identified as ARs at each concurrence time step during the cyclone trajectory. Xynthia was associated with an AR during its intensification phase until landfall, suggesting that the AR may have contributed to its intensification (Liberato et al., 2013).

Figure 3 Example of the detection and tracking of an EC (Xynthia, 26 February 2010 at 18:00 UTC to 1 March 2010 at 00:00 UTC) concurrent to an AR in ERA5. Red crosses mark the location of the cyclone at times when it concurred with an AR. Black dots indicate no concurrence. Shaded areas depict the regions identified as ARs at each concurrence time step during the cyclone trajectory.

We further evaluate the internal variability of the CMIP6 concurrences by analysing the spread within the multi-model ensemble. Additionally, we assess the inter-annual variability of the concurrences by calculating the standard deviation at each time step of the rate of coincidence between cyclones and ARs over the 29 full extended winter seasons (ONDJFM).

Figure 4 Rate of coincidence between ECs and ARs (a) and non-EC and ARs (b) for ERA5 and the historical runs of CMIP6 models (1980–2009) during the extended winter (ONDJFM). Inter-annual variability as the standard deviation of the rate of coincidence over the 30 winter seasons between ECs and ARs (c) and non-ECs and ARs (d).

4 Concurrences in present climate

In this section, we evaluate the concurrence of ECs and ARs in the ERA5 reanalyses and compare it with those obtained in climate models. Figure 4a shows the rate of coincidence between ECs and ARs, indicating the fraction of ECs that are associated with an AR (see Sect. 3.3) as a function of time from the MDP. ERA5 shows a maximum rate of coincidence at around the MDP. This maximum is about 0.72, meaning that 72 % of the ECs are associated with ARs. The rate of coincidence is minimum 36 h before and after the MDP, with 55 % of the ECs associated with an AR. These results agree with the ones obtained by Eiras-Barca et al. (2018), despite changes of up to 0.05 in the rate of coincidence that were probably due to differences in the tracking method and data used. The curves of the CMIP6 models show a similar shape to that of ERA5, with a maximum rate of coincidence centred at the MDP and two minimums 36 h before and after the MDP. Some of the models, specifically CMCC, MPI-HR, MPI-LR, and MIROC6, have higher rates of coincidence for almost the whole lifetime of the ECs. CMCC and MPI-LR show a maximum of around 80 % of ECs associated with an AR, which is 10 % more than ERA5. In contrast, EC-Earth3 and NorESM2-MM show very similar rates of coincidence along the EC lifetimes, with differences up to 0.02 with respect to ERA5. Figure 4c shows the evolution of inter-annual variability as the standard deviation of the rate of coincidence between ECs and ARs over the 29 winter seasons at each time from −36 to +36 h from the MDP. CMIP6 models and ERA5 show an inter-annual variability between 0.06 and 0.18 that is higher at the first and last time steps from the MDP. CMIP6 models reproduce similar concurrence variability to ERA5, and differences in concurrence rate between them and ERA5 are within the internal variability, computed as the spread of the multi-model ensemble. However, MPI-LR has a higher rate of coincidence with respect to ERA5 and the rest of the models 36 h before the MDP.

Figure 4b shows the rate of coincidence of non-ECs and ARs. ERA5 has a much smaller variation in the rate of coincidence for non-ECs than for the ECs. There is a maximum in the coincidence rate for ERA5 at the MDP, around 50 %, with two minimums if around 0.45 % at 36 h before and after the MDP of the cyclones. Thus, there is a relatively small variation in the coincidence rate throughout the lifetime of the non-ECs. These results for non-ECs also agree with the ones obtained by Eiras-Barca et al. (2018). All models resemble ERA5 regarding the shape of the curve. However, CMCC has a higher rate of coincidence throughout the lifetime of the cyclones of up to 0.02 with respect to ERA5. In contrast, EC-Earth3 and MIROC6 have lower rates of coincidence of up to almost 0.1 with respect to ERA5. The models that have the most similar rates with respect to ERA5 are MPI-HR, MPI-LR, and NorESM2-MM. We emphasize that caution is needed when interpreting these model performances, as they may be influenced by internal variability due to the use of only a single member. Additional sensitivity tests of the impact of resolution on ERA5 are available in the Supplement. Figure 4d shows the inter-annual variability for the rate of coincidence between non-ECs and ARs. In this case, CMIP6 models and ERA5 show an inter-annual variability between 0.08 and 0.12, showing larger values at the first and last time steps for some models, reaching 0.18 for MPI-LR. CMIP6 models reproduce similar concurrence variability to ERA5, and the differences in concurrence rate between them and ERA5 are within their inter-annual variability.

In summary, for ECs, the coincidence rate maximum model difference with respect to ERA5 occurs during the MDP, reaching around 0.1. Even though this represents almost 15 % of the ERA5 value, the inter-annual variability of the datasets is comparable. For non-ECs, the different models show a more constant difference with ERA5 with a maximum of 0.05 or around 11 % of the ERA5 value but smaller than the inter-annual variability. Overall, the models reproduce the qualitative features of the life cycle of the rate of coincidence between both ECs and non-ECs and ARs. The quantitative differences with ERA5 are generally smaller than the datasets' internal variability.

Figure 5 Rate of coincidence between ECs and ARs (a) and non-ECs and ARs (b) for the historical runs and ERA5 (1980–2009) and the future scenario runs (2070–2099) during the extended winter (ONDJFM). Solid lines show the multi-model mean, shape points are individual models, and shading shows the inter-model spread. The difference in the multi-model mean rate of coincidence between the forcing scenarios and historical runs for ECs and ARs (c) and non-EC and ARs (d) is also shown. Solid lines show the multi-model mean difference, the shading shows the inter-annual variability of the multi-model ensembles in future scenarios and the historical data, and the dashed line shows the inter-annual variability for ERA5.

5 Future projections

In this section, we analyse changes in the concurrence of cyclones and ARs in CMIP6 models between the historical simulations (1980–2009) and the three scenario simulations (2070–2099).

5.1 Changes in concurrence frequency

Figure 5a and b show a clear increase in the rate of coincidence of ARs with ECs and non-EC events across all warming scenarios. In the SSP5-8.5 scenario, the multi-model mean rate of coincidence increases by approximately 0.14 throughout the lifetime of ECs and 0.12 for non-ECs. For the SSP1-2.6 and SSP2-4.5 scenarios, the increases are around 0.07 and 0.08 for ECs, respectively, and 0.05 and 0.07 for non-ECs.

Table A2 in Appendix A4 shows the number of EC and non-EC tracks and Table A3 shows the number of AR tracks detected for ERA5 and the CMIP6 models in each scenario. Hence, these tables also show the changes (in parentheses) in the occurrence of the features individually in future scenarios compared to historical ones. The number of ECs decreases in all future scenarios with respect to the historical data and for all models, with the largest decrease for the high-emission scenario (SSP5-8.5). The decrease in the number of non-ECs is smaller when compared to the EC, although for SSP1-2.6 not all the models show a decrease in non-ECs, but for SSP2-4.5 and SSP5-8.5 all the models project a decrease. Regarding the number of AR tracks, in four models (EC-Earth3, NorESM2-MM, MIROC6, CMCC-ESM2) there is an increase in all emission scenarios. However, the changes are not linear, meaning some of the models (MPI-LR, MPI-HR, NorESM2-MM, MIROC6) depict a larger number of ARs in the SSP1-2.6 scenario than in the SSP5-8.5. The non-proportional increase in the number of AR tracks with warmer scenarios at the end of the century could be explained by the nature of the AR tracking, which detects ARs in the Laplacian of the IVT (Sect. 3.2). However, Table A3 only shows the number of AR tracks and does not account for the AR duration or extension, which might result in an increase in AR activity under a warmer climate (Zhang et al., 2024; O’Brien et al., 2022). Despite not finding a strong increase in ARs' tracks, the total number of times an AR is detected in the surroundings of both types of cyclones nonetheless shows a clear increase in all scenarios (Fig. 5a, b), with the increase becoming larger with higher warming.

The increase in the rate of coincidence between cyclones and ARs is partly the result of a decrease in the total number of cyclone time steps and a direct increase in the number of concurrent cyclone and AR time steps. To clarify this, we calculated the absolute number (including all CMIP6 models) of cyclones concurrent with ARs at the MDP and it increases across all scenarios: 4.6 %, 7.1 %, and 6.3 % for SSP1-2.6, SSP2-4.5, and SSP5-8.5, respectively. In addition, the absolute number of cyclone time steps also at the MDP decreases: −3.6 %, −5.2 %, and −12.3 % for SSP1-2.6, SSP2-4.5, and SSP5-8.5, respectively. Hence, the combined effect of an increase in AR occurrence and a decrease in the number of cyclones (with the latter being especially relevant for SSP5-8.5) explains the rise in the rate of coincidence as the level of warming increases (Fig. 5a, b). These findings suggest that changes in the characteristics of ARs, cyclones, or their interactions may be driving the observed changes, rather than the result being merely a statistical artefact of more cyclones and ARs occurring individually.

The spread across CMIP6 models is the primary source of uncertainty when evaluating changes in the rate of coincidence and exhibits overlap between different scenarios (Fig. 5a, b). The choice of model thus has a significant impact on the results when quantifying the increase in cyclone–AR compound events. Nevertheless, the inter-model spread of the SSP5-8.5 scenario is well separated from that of the historical period for both ECs and non-ECs (Fig. 5a, b). The inter-model spread further exceeds inter-annual variability for the models we analyse (Fig. 5). Indeed, the inter-annual variability of the model ensemble is consistently smaller than the multi-model mean differences between any of the three warming scenarios and historical data for ECs (Fig. 5c) and systematically smaller than the multi-model mean differences for non-ECs in the SSP5-8.5 scenario and at most lead times for the SSP2-4.5 scenario (Fig. 5d). Our results thus point to anthropogenic climate change exerting a strong influence on the increase in coincidence rates of cyclones and ARs, which is most evident for higher warming levels and for the latter cannot be attributed to natural variability or inter-model spread.

Figure 6 Mean IVT-max for the historical runs and ERA5 (1980–2009) and the future scenario runs (2070–2099) of ARs associated with ECs (a) and ARs associated with non-EC (b) during the extended winter (ONDJFM). Solid lines show the multi-model mean, shape points are individual models, and shading shows the inter-model spread. The difference in the multi-model mean IVT-max between the forcing scenarios and historical run for ECs and ARs (c) and non-EC and ARs (d) is also shown. Solid lines show the multi-model mean difference, the shading shows the inter-annual variability of the multi-model ensembles in future scenarios and the historical data, and the dashed line shows the inter-annual variability for ERA5.

5.2 Changes in AR intensity

The intensity of ARs can be quantified in different ways. Here we use the maximum value of IVT (IVT-max) within the detected AR. We use this metric because it represents a good proxy for the total transport of moisture of an AR and is not dependent on the boundaries of the detected AR (Ralph et al., 2017). Moreover, the IVT-max is usually located at the core of the AR, and this makes this variable also a good proxy for AR extension and duration as ARs reaching higher IVT values at their core tend to be larger and last longer (Guan et al., 2023). Our results indicate that the IVT-max of ARs associated with ECs and non-ECs increases proportionally to the warming in future scenarios (Fig. 6). When comparing the IVT-max between scenarios, it is particularly relevant that under the SSP5-8.5 the multi-model mean remains above 1250 kg m s⁻¹ for more than 48 h. This extended duration suggests that on average ARs under these conditions will primarily be hazardous, as classified by the AR scale from Eiras-Barca et al. (2021). For ARs with non-EC events, the IVT-max is constant throughout the entire cyclone lifetime, suggesting a uniform inflow of atmospheric moisture transport. On the other hand, for ARs associated with ECs, the IVT-max reaches its peak around the MDP of the cyclone's life cycle, showing the maximum moisture transport during the cyclone's most active phase (Fig. 6a, b). The AR intensity for ERA5 is larger than any model for the historical period because ERA5's resolution is almost 4 times higher than the CMIP6 models and attains larger values of IVT-max. Additional sensitivity tests of the impact of resolution on the ERA5 results are available in the Supplement. Despite this, the ERA5 and CMIP6 curves show the same qualitative behaviour during the cyclones' lifetime, and the inter-annual variability in both datasets is of the same order of magnitude (Fig. 6).

When evaluating uncertainties for changes in the AR intensity, inter-annual variability is smaller than the inter-model spread. Both are dependent on model resolution, as the IVT-max magnitude is sensitive to resolution, and all datasets are treated on their native grids. For ARs with ECs and non-ECs, the inter-model spread of both the SSP2-4.5 and the SSP5-8.5 scenarios is well separated from that of the historical period (Fig. 6a, b). Moreover, the difference in the multi-model mean between the future scenarios and the historical period is more pronounced than the inter-annual variability for all three scenarios and for both ECs and non-ECs (Fig. 6c,d). Thus, the IVT-max increases in all scenarios, and our results point to this being larger than both inter-model spread and internal model variability in most scenarios, and thus it is ascribable to anthropogenic warming.

Figure 7 Mean SLP for the historical runs and ERA5 (1980–2009) and the future scenario runs (2070–2099) of ECs with ARs (a), non-EC with ARs (b), ECs without ARs (c), and non-ECs without ARs (d) during the extended winter (ONDJFM). Solid lines show the multi-model mean, shape points show individual models, and the shading shows the inter-model spread. The difference in the multi-model mean SLP between the ECs with AR and without AR (e) and non-EC with AR and without AR (f) is also shown.

5.3 Changes in cyclone intensity

In this section, we examine the cyclone core SLP to assess changes in the intensity of ECs and non-ECs with the presence of ARs or lack thereof under climate change. The presence of ARs influences the evolution of SLP through the life cycle of the two types of cyclones differently, but in both cases when an AR is present the cyclone is deeper (Fig. 7e, f). The concurrence of ECs with ARs makes the ECs substantially deeper, especially before the MDP, where ECs are between 2.5 hPa and almost 10 hPa deeper. After the MDP the influence of the AR on the cyclone intensity decreases. Thus, ARs play an important role before the EC MDP, making these storms deeper from an earlier stage and resulting in longer-lasting and more intense cyclones. The particular influence in the first half of the cyclone life span suggests that moisture brought by the AR plays a key role in the EC intensification. On the other hand, the non-ECs also get deeper when occurring together with an AR, with a deepening between 2 and 7 hPa, having its peak also just before the MDP. The results from ERA5 show similar behaviour for both types of cyclones when compared to the models (Fig. 7b, d). Before the MDP, the models tend to simulate lower SLP for ECs with ARs and higher SLP for ECs without ARs. After the MDP, the models generally simulate higher SLP for both ECs with and without ARs. For non-ECs, the models have higher SLP values after the MDP compared to ERA5. However, the ERA5 values fall within the ensemble spread of historical values, indicating that they are within the uncertainty range of the models.

The influence of climate change on cyclone intensity for any of the four types of compound events is very limited, as the historical and scenario lines, along with their respective spreads in Fig. 7, are close to each other or overlap. All three forcing scenarios show the same life cycle behaviour and intensity changes as the historical ensemble. Figure A1 shows that the difference in the multi-model mean between the forcing scenarios and the historical is negligible (less than 1 hPa). It is slightly larger only for the ECs without ARs, probably due to the reduced number of these events (but as a result their inter-annual variability is also larger). The < 1 hPa shift can be compared to the models' inter-annual variability of between 1 and 2 hPa (Fig. A1). Similarly, the multi-model spreads of all scenarios are about ± 4 hPa (Fig. 7). Thus, the models do not show a robust change in cyclone intensity when averaged across the North Atlantic basin under climate change scenarios.

Figure 8 Rate of coincidence between ECs and ARs for ERA5 (a) and the multimodel-mean of the CMIP6 historical runs (b). Differences between future projections and historical periods of the three scenarios: SSP1-2.6 (b), SSP2-4.5 (c), and SSP5-8.5 (d). Dots denote where all the CMIP6 models agree on the sign of change.

5.4 Changes in the location of concurrent ECs and ARs

To assess the spatial distribution of changes observed across the North Atlantic, we show the rate of EC-AR coincidence over a 3° spherical cap. In the historical period, the rate is highest along the climatological North Atlantic storm track, extending from the western to the northeastern part of the basin (Fig. 8a, b). CMIP6 models reproduce closely this spatial pattern. Differences between future projections and the historical period show an increase in the rate of coincidence, all the stronger for higher warming levels. The high-emission scenario shows the largest increase and a high model agreement across almost the entire domain. The southern part of the domain shows a noisier pattern and weaker model agreement due to the reduced number of events in the area. Only the high-emission scenario shows a clear increase in the rate of coincidence over Europe, showing model agreement over the British Isles, northern France, and the Iberian Peninsula. In other words, a larger proportion of landfalling ECs and ARs at their maximum deepening point is expected under the SSP5-8.5 scenario. This highlights the possibility of an increase in wet and windy extremes in western Europe in a high-emission future.

6 Conclusions

We have used six global climate models (GCMs) participating in CMIP6 to evaluate the change in the concurrences of cyclones and atmospheric rivers (ARs) in three different future scenarios under climate change in the North Atlantic. We have compared the performance of the models using the ERA5 reanalysis. Our main findings are summarized as follows.

• For the present period, in ERA5 nearly 72 % of the explosive cyclones (ECs) are associated with an AR at the maximum deepening point (MDP). This higher rate of coincidence around the MDP indicates that ARs are more likely to occur when the EC is at its peak deepening stage. Despite some biases in the magnitude of the coincidence rate, CMIP6 models exhibit qualitatively similar concurrence rates to ERA5. In contrast, the concurrence rate over the lifespan of non-explosive cyclones (non-ECs) shows much smaller variation compared to ECs.

• In future scenarios, there is an increase in the rate of coincidence between ECs and ARs, with the magnitude of the increase being proportional to the level of warming. All models agree on the sign of the change. In the high-emission scenario (SSP5-8.5), there is a maximum increase in the rate of coincidence from 0.72 to 0.87. Even in the best-case scenario (SSP1-2.6), there is an increase of around 0.07. For non-ECs, there is an increase in the rate of coincidence, ranging from 0.05 to 0.13 across scenarios. The increase in SSP5-8.5 and SSP2-4.5 scenarios is larger than the inter-annual model variability for both ECs and non-ECs, but the inter-model spread of the historical period is only well separated from that of the SSP5-8.5 scenario.

• In all warming scenarios there is an increase in AR intensity. This is larger than the inter-annual variability for ARs associated with both ECs and non-ECs. Moreover, the inter-model spread of both the SSP2-4.5 and the SSP5-8.5 scenarios is well separated from that of the historical period. Under the SSP5-8.5 scenario, the maximum integrated vapour transport (IVT) of ARs associated with ECs is projected to exceed 1250 kg m⁻¹ s⁻¹ for more than 48 h on average, indicating exceptional and hazardous AR conditions.

• ECs are deeper when associated with an AR. The presence of an AR is associated with ECs deepening earlier compared to those without, making them deeper for a longer period and potentially more hazardous. Non-ECs are also deeper when having an AR in their surroundings, but their intensification is lesser compared to the ECs with ARs.

• The concurrence of ARs with ECs will be more frequent in the North Atlantic basin in future climates. There is an increase in both agreement among models and the magnitude of the change with the degree of warming, with the SSP5-8.5 scenario showing the largest change. Under this most severe scenario, Europe is exposed to this increase, particularly the southern Iberian Peninsula, the British Isles, France, and Scandinavia.

Our results for concurrence rates of cyclones and ARs in the present climate are broadly consistent with Eiras-Barca et al. (2018). One key difference is that they calculated the MDP using a 24 h time window, while we used a 6 h window. These differences in time intervals, as well as detection and tracking configurations for cyclones and ARs, may account for the minor differences in concurrence rates. Furthermore, our study uses the latest ECMWF reanalysis, ERA5 (Hersbach et al., 2020), whereas their study used the earlier ERA-Interim dataset.

Previous studies found an increase in the IVT of ARs under climate change (Payne et al., 2020), as well as an increase in their frequency (Espinoza et al., 2018; Wang et al., 2023; Ramos et al., 2016). These results align with the increase in IVT-max detected for ARs associated with ECs and non-ECs and the increase in concurrences between cyclones and ARs, partly driven by an enhanced AR frequency. Moreover, other studies have found an increase in the frequency and severity of extratropical cyclones under climate change, mainly over the British Isles, and an eastward extension of the storm track activity over Europe (Priestley and Catto, 2022; Zappa et al., 2013; Seiler and Zwiers, 2016a). Our results show a generalized increase in compound events of ECs with ARs in most of the North Atlantic basin. A robust increase in concurrence over the British Isles, Iberia, and northern France is only observed under the most severe climate change scenario. However, we did not detect a clear upward trend in the individual frequency of EC or AR tracks across the entire North Atlantic basin. This apparent contradiction suggests that changes in the characteristics or dynamics of ECs and ARs, rather than their frequency, may be driving the observed increase in concurrence. This is a significant finding that needs further investigation, as the underlying physical mechanisms for this increase remain unclear.

Our analysis has limitations that should be acknowledged. The main constraint is the reduced number of CMIP6 models and members used. The number of models or ensemble members used is limited by the availability of data on multiple vertical levels at the 6-hourly resolution that are necessary to compute IVT. We acknowledge that using only one member per model does not facilitate a comprehensive model intercomparison; more members for each model would be needed to adequately assess model uncertainty, specifically the biases of the models relative to reanalysis data. We chose all models and members from CMIP6 where these variables were available for the historical, SSP1-2.6, SSP2-4.5, and SSP5-8.5 experiments. In particular, we deemed it important to go beyond the high-emission scenario (SSP5-8.5) and also look at the implications of lower warming levels. While most of the results presented indicate a stronger signal for the highest-emission scenario (SSP5-8.5), this scenario has been deemed unrealistic (Hausfather and Peters, 2020). Therefore, we emphasize that our results should be interpreted with consideration of various scenarios. A further caveat is that we used a single tracking algorithm, namely the TempestExtremes software (Ullrich et al., 2021). While we have compared our results for the present period with a previous study that uses different tracking algorithms (Eiras-Barca et al., 2018), this does not detract from the fact that our results depend on the detection and tracking method (Neu et al., 2013). Furthermore, future studies should aim to delve deeper into isolating the dynamic signal from the thermodynamics of the climate change response. Finally, future work should explore future changes in the wet and windy extremes associated with the compound meteorological events investigated here.

Appendix A

A1 Information about the CMIP6 models

Mauritsen et al. (2019)Müller et al. (2018)Seland et al. (2020)Döscher et al. (2022)Cherchi et al. (2019)Tatebe et al. (2019)

Table A1 Description of the CMIP6 models and members used for the historical, SSP1-2.6, SSP2-4.5, and SSP5-8.5 scenarios.

A2 TempestExtremes code for detecting and tracking extratropical cyclones

To identify extratropical cyclones, we use the executable DetectNodes from the TempestExtremes tracking algorithm (Ullrich et al., 2021; Ullrich, 2020), which recognizes candidate “nodes” corresponding to local minima in the SLP field. Subsequently, we employ StitchNodes also from the TempestExtremes tracking algorithm (Ullrich et al., 2021; Ullrich, 2020) to connect these candidate nodes into tracks. The tracking codes set up with the parameters used are as follows:

A3 TempestExtremes code for detecting and tracking ARs

For the detection and tracking of ARs, we use the executables from the TempestExtremes tracking algorithm (Ullrich et al., 2021; Ullrich, 2020) DetectBlobs to detect ARs and to connect ARs or “blobs” we use the executable StitchBlobs with the following parameters:

A4 Number of cyclone tracks and AR tracks detected in ERA5 and CMIP6

Section 3.1 describes how cyclones are detected, tracked, and classified as ECs or non-ECs. The result of this process is the track of each cyclone. In the following table we summarize the number of individual EC and non-EC tracks for each dataset and the percentage difference from the historical period.

Table A2 Number of EC and non-EC tracks detected in each dataset. In parentheses we provide the percentage difference between the future scenario and the historical period for each model.

Section 3.2 describes how ARs are detected and tracked. The result of this process is the track of each AR. In the following table we summarize the number of individual AR tracks for each dataset and the percentage difference from the historical period.

Table A3 Number of AR tracks detected in each dataset. In parentheses we provide the percentage difference between the future scenario and the historical period for each model.

A5 Cyclone intensity differences between future scenarios and historical simulations

Figure A1 The difference in the multi-model mean SLP between the forcing scenarios and historical runs is shown as solid lines, the inter-annual variability of the multi-model ensembles of mean SLP is shown with shading, and the inter-annual variability for ERA5 is shown with the dashed line for ECs with ARs (a), non-ECs with ARs (b), ECs without ARs (c), and non-ECs without ARs (d).