Extratropical storms shape midlatitude weather and vary due to the slowly evolving climate and the rapid changes in synoptic conditions. While the influence of each factor has been studied extensively, their relative importance remains unclear. Here, we quantify the climate's relative importance in mean storm activity and individual storm development using 84 years of ERA-5 data and convolutional neural networks. We find that the constructed model predicts over 90% of the variability in the mean storm activity. However, a similar model predicts about a third of the variability in individual storm properties, such as maximum intensity, showing their variability is dominated by synoptic conditions. Isolating the impact of present-day climate change on individual storms shows it contributes to about 0.1% for storm-intensity variability, whereas their contribution to storms' heat-anomaly variability is over three times greater, highlighting that focusing on variables directly tied to global warming offers a clearer attribution pathway. When attempting to attribute extreme extratropical weather events, often associated with cyclones and anticyclones, to changes in Climatic Forcing (CF) (Philip et al., 2022;Ginesta et al., 2023), a major uncertainty arises from the Synoptic Variability (SV) of storms. Over extended periods, the climate is expected to be influenced primarily by CF, such as solar insolation, sea surface temperature anomalies, and the greenhouse effect. However, due to the chaotic nature of weather, the dynamics of individual storms exhibit large SV, which results from the fine structure of the meteorological conditions during their formation. Although SV and CF play significant roles in determining the intensity and trajectory of individual storms, their relative importance has yet to be quantified. Bridging this knowledge gap can uncover the importance of the dynamical response of storms to changes in CF when studying individual weather events. Therefore, in this study, using a statistical approach, we attempt to quantify the importance of CF to the growth of individual storms and the climatology of storm activity. The effect of CF on midlatitude climate and weather has traditionally been studied by decomposing the flow into a mean (spatially or temporally averaged) and an eddy component. In this framework, the mean flow largely reflects climatological features such as the jet stream, while eddies describe synoptic-scale weather systems. Hence, given a particular CF, serving as a boundary condition, one expects a characteristic midlatitude flow pattern, consisting of a mean flow and a spectrum of eddies. Our central goal is to quantify how strongly these synoptic eddies (storms) depend on the mean flow, both in the time-averaged sense (climatological storm activity) and at the scale of individual storm growth. The interaction between the mean flow and midlatitude storms is well studied, both from theoretical and idealized perspectives (Charney, 1947;Eady, 1949;Phillips, 1954;Harnik & Chang, 2004), climatological alternation in the intensity of the storm tracks (Nakamura, 1992;Schemm & Riviere, 2019), and the effect on individual storms (Orlanski & Katzfey, 1991;Orlanski & Chang, 1993;Riviere & Joly, 2006;Schemm et al., 2020;Hadas et al., 2023;Hadas & Kaspi, 2025). Despite this extensive literature, the fundamental question-namely, the degree to which the mean flow dictates the characteristics of individual storms, has not been addressed quantitatively. In this study, we use 84 years of ERA-5 reanalysis (Sec. 2.1) and tracks of about 100,000 cyclones and 50,000 anticyclones (Sec. 2.2) to quantify the extent to which storm activity is determined by the mean flow from a statistical perspective. To investigate this relationship, we use two modeling approaches. First, the climatological model (Sec. 2.3) tests whether the link between the mean flow and the climatological storm activity is unique (Sec. 3.1). Second, the single storm model (Sec. 2.4) evaluates how strongly the mean flow influences the properties of individual storms (Sec. 3.2). Finally, we discuss the implications of our findings for the attribution of weather events and climatic changes to current anthropogenic climate change (Sec. 3.3). Data from the ERA-5 reanalysis by the European Center for Medium-Range Weather Forecasts (Hersbach et al., 2020;Soci et al., 2024), spanning from 1940 to 2023, is employed for current climate assessment. The ERA-5 reanalysis provides estimates of atmospheric variables with a horizontal resolution of 31 km and 137 vertical levels. (SLP) data is utilized to track cyclones and anticyclones, while three-hourly temperature (T ), specific humidity (Q), zonal wind (U ), and meridional wind (V ) at pressure levels 300, 500 and 850 hPa are utilized to estimate the atmospheric state and its climatology. A feature point tracking algorithm (Hodges, 1995;Tamarin & Kaspi, 2016, 2017) is applied to SLP data to identify and characterize extratropical cyclones and anticyclones. To reduce noise, the data is smoothed to a T63 resolution. The background, defined as zonal wavenumbers 0-4, is removed to isolate the synoptic scale dynamics. Only SLP anomalies deeper than 8 hPa are tracked to focus on significant events. After identification, systems are tracked, and their location and intensity are recorded if they persist for more than 48 hours and move more than 500 km westward. Additionally, systems that reach peak intensity over topography higher than 1 km are filtered out. The tracking algorithm returns the positions and depth (in hPa) of the storms from genesis to lysis. The output of the tracking algorithm is used to quantify the growth of storms using four metrics: the maximum intensity, defined as the maximum depth of the storm normalized by the sinus of latitude (Sanders & Gyakum, 1980), the growth time, defined as the time between identification and maximum intensity, and the meridional and zonal rate of propagation of the storm from genesis to maximum intensity. In addition, the maximum 850 hPa temperature anomaly associated with a cyclone is defined as the maximum temperature anomaly around the cyclone center in a composite of 15  +/- in the zonal direction and 10  +/- in the meridional direction. An anomaly is defined as the deviation from the 84-year grid and time of the year average.
The connection between the mean flow and climatological storm activity is treated as an image-to-image regression problem , solved using Convolutional Neural Network (CNN), due to its ability to efficiently learn spatial hierarchies and patterns in structured data. The model input is the three-dimensional distribution of U , V , Q, and T over the midlatitudes (plus and minus 20  -80  in for the NH and SH, respectively), where bar represent 90-day moving mean (45 days in each direction). The output is the distribution of vertically integrated Eddy Kinetic Energy: where the bar is defined the same as above, the prime is defined as a deviation from a 10-day moving mean to focus on synoptic activity, and the integration is performed between 300 and 850 hPa. The EKE is calculated over the same domain. divided into trainvalidation-test set in ratios of 0.5  0.2  0.3. An in-depth discussion of the data appears in Sec. S1. Averaging the data over time substantially simplifies the regression task by reducing the effective image size in the frequency domain. Because atmospheric waves typically follow a dispersion relation-where slower waves tend to be larger-this temporal averaging further diminishes the impact of small-scale waves on variability (Fig. S1 a-c). This simplification achieved can be quantified by the number of principal components required to capture 90% of the variability in the dataset. Analyzing this reveals that filtering reduces the number of required components by 80% (Fig. S1d). In addition, time averaging significantly simplifies the dynamics, as it nullifies (or at least significantly reduces) the contribution of many physical processes to the overall system (Vallis, 2017). Consequently, much simpler models can simulate the relation between mean and eddy flow, relative to numerical weather prediction models (Lam et al., 2023;Bi et al., 2023;Kochkov et al., 2024). The model is built out of a series of convolutional layers (Krizhevsky et al., 2012;Guan et al., 2023). Given an input, it predicts the shape and scale parameters of the underlying Gamma distribution, which are optimized using negative log-likelihood loss (Eq. 2 in the supporting information), similar to Guillaumin & Zanna (2021). An in-depth discussion of the model, including the loss function, structure, hyperparameters, and steps taken to avoid overfitting, is provided in Sec. S2. The extent to which the mean flow sets the growth of individual storms is tested by constructing a model that predicts the average outcome of a storm's growth based on the mean flow it experiences at genesis. The model input is the three-dimensional composites of U , V , Q, and T around the storm at genesis. The model output is the outcome of growth, namely the maximum intensity, the growth time, and the meridional and zonal rate of propagation of the storm from genesis to maximum intensity (Sec. 2.2). Tracking all cyclones and anticyclones between the years 1940 and 2023 results in approximately 100,000 cyclones and 50,000 anticyclones. An in-depth discussion of the data appears in Sec. S3. The model consists of multiple convolutional layers, each followed by max-pooling to reduce dimensionality, thus performing feature extraction. The resulting features then pass through two dense layers to regress on the storm properties. To assess the model's robustness, cross-validation is performed by reshuffling the train-validation-test splits and training/testing the model 50 times. Additionally, we compare the CNN to a dense neural network (no feature extraction) and a random forest (Breiman, 2001), implemented via RandomForestRegressor and MiniBatchKMeans from scikit-learn (Pedregosa et al., 2011). Further details on the model architecture, hyperparameters, and this comparison can be found in Sec. S4. Climatological storm activity becomes largely independent of synoptic variability (SV) when averaged over sufficiently long timescales. Consequently, if the mean atmospheric state captures most of the climate forcing (CF) information, it should explain much of the variability in climatological storm activity. To test this, we evaluate how well a climatological model can predict EKE maps using only mean-flow fields (Sec. 2.3). Comparing the model's predictions with the climatological EKE averaged from September 28 to December 27, 2017 shows that the model successfully reproduces major features such as the EKE "hot spots" along the Northern Hemisphere Atlantic and Pacific storm tracks, as well as many smaller-scale seasonal and interannual variations. To quantify performance, we use the coefficient of determination: where y i and y pred are the true and predicted EKE at each grid point over time (test set), and y is the overall mean EKE. As shown in Fig. 1d, the model achieves R 2 values above 90% in most regions (shading), especially in regions of high climatological storm activity (contours). Overall, about 93% of the total EKE variability is explained on average. These results confirm that the mean flow contains the core CF information essential for determining the climatology of storm activity. Shortening the averaging window for EKE reduces the model's predictive skill. However, as long as the window remains longer than the synoptic scale (roughly 10 days), the model still explains above 50% of the variability (Fig. S3). Hence, the mean flow dominates storm-activity variability up to the scale of a single synoptic cycle. Moreover, model performance is relatively insensitive to which atmospheric variables are provided (Fig. S4), with a single variable alone accounting for over 90% of the variance. This finding reflects that temporal averaging enforces consistency among mean flow variables. For the dry variables, the thermal wind balance quantitatively demonstrates this, as it allows the inference of the mean atmospheric wind profile based on the mean temperature profile and vice versa. The consistency with the mean atmospheric water content can be understood qualitatively since atmospheric dynamics influence the injection, removal, and distribution of water vapor, which couples back through diabatic effects, and therefore, their climatologies are tightly connected. While the mean flow largely sets the climatology of storm activity, individual storms can also be influenced by SV. To quantify the relative importance of the mean flow, we construct a model that predicts storm growth based solely on the mean flow X (Sec. 2.4).
Because X encapsulates most of the relevant CF information for storm climatology (Sec. 3.1), it is expected to provide an effective representation of CF's influence on individual storm dynamics. Meanwhile, SV is captured by the remainder of the atmospheric initial state, X  = X  X, where X is the full initial state. Since storm dynamics also depend on X  , there is a range of possible growth outcomes for any given X (Mana & Zanna, 2014). Consequently, a model that relies only on X, cannot distinguish among these different outcomes. Instead, with a mean square error loss, it learns the single best-fit prediction for the observed distribution of outcomes, i.e., the mean of that distribution (DelSole & Tippett, 2022). If X strongly drives storm growth variability, each X i leads to a relatively narrow distribution of possible outcomes, separated from other X j =i , yielding a lower model error. Conversely, if SV dominates, each X is associated with a broad distribution of outcomes, which diminishes the predictive skill of the mean-flow-only model, resulting in higher model error. To illustrate this concept, we approximate the distribution of maximum storm intensities under different mean flow conditions by clustering storms into four groups using K-means applied to X . The resulting distributions are distinct, indicating that even coarse sorting by mean flow yields differences in growth outcomes. However, the within-cluster spread, which is approximately the SV, is large relative to the differences between the cluster means, which is due to the variability of CF. To go beyond this qualitative demonstration, we train a neural network to learn the mapping between X and the expected storm growth . The network's predictions closely match the average maximum intensity observed within each cluster . The model's mean square error, normalized by the total variance in storm intensity, then provides a quantitative estimate of the contribution of SV to storm growth variability. One minus this value gives the contribution of climatic variability, and also corresponds to the model's R 2 . Fig. 3 compares the single-storm model's residual distribution against the overall distribution to illustrate the relative contributions of mean flow and SV. As indicated by the similar widths of the residual (dashed lines) and the original distribution (solid lines), most of the variability in individual storm properties cannot be explained by the mean flow. Quantitatively, for cyclones (anticyclones), the climatic variability accounts for only about 30% (38%) of the the maximum intensity variability , 7% (7%) of the growth-time variability , 34%(25%) of the zonal propagation rate variability, and 25% (23%) of the meridional propagation rate variability . The low climatic variability values show that most variability among storms on Earth arises from SV rather than from CF. A series of sensitivity tests confirms the robustness of our single-storm model results. First, using only a Dense Neural Network (DNN) or a nonparametric Random For- est regressor yields nearly identical performance (Table S1). Second, training the model on at least 20,000 samples appears sufficient for convergence (Fig. S5). Finally, the model's predictive skill increases significantly only when the averaging window X is shorter than about 10 days-an interval that begins to capture more synoptic variability (Fig. S6). Due to anthropogenic forcing, the mean climate has changed in the last decades Vallis et al. (2015); Forster et al. (2021); Woollings et al. (2023). In this section, the portion of single-storm variability attributable to climate change is isolated. This is achieved by comparing the predictability found in Sec. 3.2 to a modified test set in which the climate trend is "scrambled," enabling a distinction between the long-term climate signal and other climate-driven variability. Following prior studies (e.g., Woollings et al., 2023), the climate trend at each grid box is determined by linearly regressing mean flow on the year for each date. A modified test set is then created by removing the correct trend and adding the trend from a random year, thus preserving seasonal and interannual variability while eliminating the actual long-term signal. Finally, the storm model's performance (Sec. 2.4) on this modified dataset is compared with that on the original test set, and the drop in skill is attributed to the climate change trend. In this section, only the years 1980-2023 were used, in order to focus on a period where the climate trend is large (as can be assessed qualitatively from the global mean surface temperature trend, Forster et al. (2021)), and to reduce biases related to data availability (Soci et al., 2024). To estimate the uncertainty, we train and test the model 50 times, where each time we assign different storms for training, validation, and test- ing. From this comparison, the climate trend in the mean flow explains approximately 0.10(+/-0.01)% of cyclone intensity variability . Thus, accurately discerning such a weak signal poses a substantial challenge for attribution efforts, especially given the dominance of synoptic variability found in Sec. 3.2. Nevertheless, these results do not imply that climate change is undetectable in the midlatitudes. As shown by many studies and supported here , temporally averaging storm activity diminishes the role of synoptic variability. Once the long-term signal rises above natural variability, its detection becomes feasible (Priestley & Catto, 2022;Chemke et al., 2022;Chemke & Coumou, 2024;Shaw et al., 2024). Furthermore, certain event attributes respond more directly to climate change, making them more amenable to attribution (Trenberth et al., 2015;Catto et al., 2019). For example, temperature anomalies accompanying cyclones (see Sec. 2.2) are expected to show a stronger climate-change signal due to robust warming. Indeed, re-training the model to predict temperature anomalies and scrambling their trend shows that climate change currently contributes to about 0.35(+/-0.017)% of the variability , more than three times higher than for cyclone intensity. These findings underscore that, while attributing individual storm intensity changes to climate forcing remains difficult, examining variables more directly linked to a warming climate can offer a clearer attribution pathway. No new data sets were generated during the current study.
The three-dimensional atmospheric variables on pressure levels (e.g., zonal and meridional wind) are available from Copernicus through the dataset "ERA-5 hourly data on pressure levels from 1940 to present" (Hersbach et al., 2023a). Sea level pressure data is available from Copernicus through the dataset "ERA-5 hourly data on single levels from 1940 to present" (Hersbach et al., 2023b). where  (shape) and  (scale) parameters define the distribution and  is the Gamma function. The mean and standard deviation of this distribution are  and  2 , respectively. Following Guillaumin & Zanna (2021)-but substituting the Gamma PDF-the negative log-likelihood (loss) function we optimize is: where y i is the true EKE in pixel i and  i ,  i are predicted from the input. Hence, the model outputs two images per input image: one for the shape parameters and one for the scale parameters, enabling joint optimization of both moments. To compute the R 2 score, we set the model's EKE prediction to y pred = . Empirically, this approach achieves about 7% higher R 2 (92% vs. 85%) than a model trained using mean square error. Our model is built from convolutional layers, each followed by a ReLU activation function. To address the anisotropic input dimensions, we use kernels of size 3 x 9 (meridional x zonal). We apply zero-padding in the meridional direction and circular padding in the zonal direction to preserve image dimensions. Each layer has eight filters, effectively capturing spatial features. The final layer has two filters-one for the shape parameter and one for the scale parameter-and uses a softplus activation (log 1 + e x ) to ensure positivity. Systematic tuning of the number of layers and filters (Fig. S2) shows that performance increases with the total number of parameters up to around 10 4 . This improvement behaves similarly whether we add more layers or more filters. Consequently, we use five layers with eight filters each, yielding about 10 4 parameters and an R 2 above 90%. To avoid overfitting the model parameters, we monitor the validation loss after each epoch (i.e., after each pass over the training set), selecting the model checkpoint with the lowest validation loss. To prevent overfitting to hyperparameters, all final results reported in the paper are evaluated on a separate test set, which is never used during training or hyperparameter tuning. Using default settings in Scikit-Learn, LinearRegression, MLPRegressor, and Ran-domForestRegressor achieve R 2 scores of 86%, 90%, and 90%, respectively. Although these simpler models perform reasonably well, our CNN-based approach outperforms them. Furthermore, given the high dimensionality of both inputs and outputs, the CNN model runs and tunes significantly faster, making it more practical for this application. The model input is the three-dimensional composites of the mean flow, which extend +/-60  in the zonal direction and +/-30  in the meridional direction, with resolution of 6  in the zonal direction and 3  in the meridional direction. This relatively low resolution is chosen as the mean flow is associated with large-scale phenomena, and therefore increasing the resolution would not change the results significantly. Therefore, the overall shape of the input is 21 by 21 pixels with 9 channels (three variables times three vertical levels). The model output is the outcome of growth, namely the maximum intensity, the growth time, and the meridional and zonal propagation of the storm from genesis to maximum intensity. The storms' properties in the input and output are based on tracks of storms obtained through a feature tracking algorithm. Tracking all cyclones and anticyclones between the years 1940 and 2023 results in approximately 100,000 cyclones and 50,000 anticyclones, which are divided between train-validation-test in a ratio of 0.6  0.2  0.2. The input images are normalized by the grid-wise average and standard deviation. The input and output vectors are also normalized by the mean and standard deviation. The first part of the model processes the image input to perform feature extraction. The basic building blocks of this part are convolutional blocks, which transform the input using two-dimensional convolution, followed by Relu activation and Max Pooling with a kernel of 2-by-2. Three of the above layers are stacked. The first convolutional layer has 4 filters, and each time a pooling is performed, the number of filters in the next layers is doubled. The second part of the model is built out of dense layers. The output of the CNN is passed to a private dense network for each predicted property (four in total), which is built out of two dense layers with 8 nodes, followed by Relu activation. The outer layer is built out of a dense layer with a single node and linear activation, which allows regression. The performance of the model is assessed using a mean square error loss function. The mean square error is also used to assess the SV. Shortening the averaging window (Fig. S6) significantly affects the model performance only when the averaging time is shorter than the synoptic time scale, as it includes knowledge about the SV. Training the model on part of the storms in the training set shows that about 10,000-20,000 storms are required to achieve convergence (Fig. S5). The CNN's performance is compared to two baseline models: a four-layer Deep Neural Network with 16 nodes per hidden layer, and a Random Forest regressor in Scikit-Learn with 250 estimators, bootstrap enabled, a maximum depth of 30, a maximum feature fraction of 0.05, a minimum of 16 samples per leaf, and a minimum of 20 samples per split. The results of the validation appear in Tab.1 and show that the CNN, DNN, and Random Forest perform very similarly, whereas the K-Means approach while improving with more clusters, remains below the performance of the other methods even with 100 clusters. This research has been and Nimrod Gavriel for the meaningful discussion. This research has been and Nimrod Gavriel for the meaningful discussion. The authors have no conflicts of interest to disclose.
Text S1: Data for the climatological model The resolution of both the input and output images is 6  in the zonal direction, 3  in the meridional direction, and the levels 300, 500, and 850 hPa are used in the vertical. Using this resolution, the inputs and outputs are images of 20 by 60 pixels, with 9 channels for the input (three levels of three variables) and 3 channels for the output (three levels of one variable). The meridional direction of the SH data is flipped to match the NH. The data is sampled in 10-day intervals to ensure significant differences between data points while obtaining a large dataset. Overall, there are about 6,000 data points (36 images per year for each hemisphere for 84 years) divided into train-validation-test set in ratios of 0.5  0.2  0.3. The division between sets is done by splitting the data into chunks of 108 consecutive samples (i.e., about 3 years), with a margin of five samples between chunks that were not included in any data sets. Then, given that the maximum window that was used was 90 days (45 on each side), no data leak was introduced. For training and validation, data augmentation is performed by randomly shifting the central longitude pixel, leveraging the system's symmetry with respect to the central longitude. Each epoch consists of three repetitions of the input data, with each sample shifted by a different randomly chosen integer. The input images are normalized using the mean and standard deviation computed along the zonal direction and over time, ensuring symmetry is preserved in these dimensions. The output is normalized using the global standard deviation (computed across all pixels and time) to prevent feeding extremely large or small values into the model, improving numerical stability. However, the mean is not removed to maintain the positivity of the output. Notably, the loss function itself accounts for per-pixel normalization based on the input (see next section). To capture the variability in each pixel's EKE, which can differ substantially based on the mean flow, we adopt a negative log-likelihood framework that estimates variance on a per-pixel basis. Since EKE is strictly positive, we model its distribution using a Gamma distribution: p(y; , ) = y 1 e y/   () , y > 0,  > 0,  > 0. (3) The authors have no conflicts of interest to disclose. Text S1: Data for the climatological model The resolution of both the input and output images is 6  in the zonal direction, 3  in the meridional direction, and the levels 300, 500, and 850 hPa are used in the vertical. Using this resolution, the inputs and outputs are images of 20 by 60 pixels, with 9 channels for the input (three levels of three variables) and 3 channels for the output (three levels of one variable). The meridional direction of the SH data is flipped to match the NH. The data is sampled in 10-day intervals to ensure significant differences between data points while obtaining a large dataset. Overall, there are about 6,000 data points (36 images per year for each hemisphere for 84 years) divided into train-validation-test set in ratios of 0.5  0.2  0.3. The division between sets is done by splitting the data into chunks of 108 consecutive samples (i.e., about 3 years), with a margin of five samples between chunks that were not included in any data sets. Then, given that the maximum window that was used was 90 days (45 on each side), no data leak was introduced. For training and validation, data augmentation is performed by randomly shifting the central longitude pixel, leveraging the system's symmetry with respect to the central longitude. Each epoch consists of three repetitions of the input data, with each sample shifted by a different randomly chosen integer. The input images are normalized using the mean and standard deviation computed along the zonal direction and over time, ensuring symmetry is preserved in these dimensions. The output is normalized using the global standard deviation (computed across all pixels and time) to prevent feeding extremely large or small values into the model, improving numerical stability. However, the mean is not removed to maintain the positivity of the output. Notably, the loss function itself accounts for per-pixel normalization based on the input (see next section). To capture the variability in each pixel's EKE, which can differ substantially based on the mean flow, we adopt a negative log-likelihood framework that estimates variance on a per-pixel basis. Since EKE is strictly positive, we model its distribution using a Gamma distribution: p(y; , ) = y 1 e y/   () , y > 0,  > 0,  > 0. (3)