Statistical Projection of the North Atlantic Storm Tracks

A method for empirical–statistical downscaling was adapted to project seasonal cyclone density over the North Atlantic Ocean. To this aim, the seasonal mean cyclone density was derived from instantaneous values of the 6-h mean sea level pressure (SLP) reanalysis ﬁelds. The cyclone density was then combined with seasonal mean reanalysis and global climate model projections of SLP or 500-hPa geopotential height to obtain future projections of the North Atlantic storm tracks. The empirical–statistical approach is computationally efﬁcient because it makes use of seasonally aggregated cyclone statistics and allows the future cyclone density to be estimated from the full ensemble of available CMIP5 models rather than from a smaller subset. However, the projected cyclone density in the future differs considerably depending on the choice of predictor, SLP, or 500-hPa geopotential height. This discrepancy suggests that the relationship between the cyclone density and SLP, 500-hPa geopotential height, or both is nonstationary; that is, that the statistical model depends on the calibration period. A stationarity test based on 6-hourly HadGEM2-ES data indicated that the 500-hPa geopotential height was not a robust predictor of cyclone density.


Introduction
Midlatitude cyclones are connected to heavy precipitation, winds, and flooding (Ralph and Dettinger 2011) and can severely impact infrastructure (Hov et al. 2013;Field et al. 2012).They are often referred to as ''synoptic storms'' or just ''storms,'' and they transport heat from low to high latitudes and represent an important mechanism for distributing energy within the climate system (Caballero and Langen 2005;Shaw et al. 2016).Cyclones also play an important role for the hydrological cycle through redistribution of heat and moisture.In addition, they are associated with cloudiness and hence influence the planetary albedo (Field and Wood 2007).
Storms can be described by a set of features such as the time and location of their tracks; the region of cyclogenesis; their lifetimes, strengths, and spatial extents; and where they dissipate.Furthermore, geographical maps of cyclone densities can be compiled, for example, providing the mean number of storms per season and region.Many studies have been published on midlatitude storms using a vast variety of methods to describe the midlatitude storm tracks based on subdaily reanalysis data, either identifying and tracking individual cyclones (e.g., Ulbrich et al. 2009;Bengtsson et al. 2006;Pinto et al. 2009;Zappa et al. 2013b) or in terms of variance or covariance statistics (e.g., Yin 2005).Neu et al. (2012) noted that different methods sometimes give inconsistent results, which motivated them to assess a range of methods for identifying and tracking midlatitude cyclones.Different tracking methods tend to agree well on the interannual variability and geographical distribution of cyclones, but are less consistent when it comes to the total number of storms, identification of weak storms, cyclogenesis, and cyclolysis.It can be challenging to compare results from studies using different tracking schemes as they may emphasize different aspects of the cyclone life cycle and parts of the storm tracks.Cyclone detection and tracking schemes are also sensitive to the spatial resolution of the input field, and some storm features may not be sufficiently resolved in reanalysis and global climate model data.Intercomparisons of cyclone identification and tracking methods indicate that the uncertainties associated with different identification and tracking algorithms (Neu et al. 2012) are larger than the differences due to different reanalyses datasets (Hodges et al. 2011;Tilinina et al. 2013).
Global climate models (GCMs) are able to simulate midlatitude cyclones, but with some systematic biases.Recent studies of the GCMs participating in phase 5 of the Climate Model Intercomparison Project (CMIP5) have shown that the storms tend to be too weak and that in boreal winter [December-February (DJF)] the storm tracks are too zonal or are displaced southward, with too many in central Europe and too few in the Norwegian Sea (IPCC 2013;Zappa et al. 2013a).However, the position and tilt of the storm tracks have improved in some of the higher-resolution CMIP5 models relative to the previous versions included in the CMIP3 ensemble (Zappa et al. 2013a).
Cyclones exist because of baroclinic instabilities associated with temperature gradients and wind shear (Gill 1982).They tend to form in regions with a large surface temperature gradient, and their trajectories and speed are influenced by the upper-level jet stream (Shaw et al. 2016).In the Northern Hemisphere, cyclones typically occur in the vicinity of the atmospheric upper-level jet surrounding the Arctic.
Present and future changes in the position and intensity of the midlatitude jet stream and storm tracks have been a topic of recent scientific debate.Francis and Vavrus (2015) and Liu et al. (2012) hypothesized that, based on recent historical data, future Arctic amplification would result in an equatorward shift in the jet stream, with more waviness in the jet stream and more blocking events in winter.Shaw et al. (2016) described the situation as a ''tug of war'' between different physical conditions that affect the midlatitude storms and the jet stream, which include the equatorto-pole temperature gradient, the vertical temperature structure, ocean temperatures, and cloud cover.Other studies of the CMIP5 models, assuming a high greenhouse gas emission scenario [representative concentration pathways (RCP) 8.5 (RCP8.5)],(Barnes and Polvani 2015;Cattiaux et al. 2016) have alluded to the effect of greenhouse gas-induced warming in the far future (2071-2100) leading to a stronger warming of the upper troposphere near the equator and greater latitudinal temperature gradient in the upper troposphere [greenhouse gas (GHG) effect], which results in a northward shift of the midlatitude jet stream in all seasons but winter.In winter, the Arctic amplification and greenhouse gas-induced warming have competing effects on the position of the jet stream (Cattiaux et al. 2016;Oudar et al. 2017), resulting in a large spread in the response of global CMIP5 RCP8.5 climate model projections (Barnes and Polvani 2015).
Earlier studies based on CMIP3 found that winter storm tracks over the North Atlantic Ocean were projected to shift poleward (Yin 2005) and extend farther east (Woollings et al. 2012) during the twentyfirst century.Two more recent studies (Cattiaux et al. 2016;Peings et al. 2017) found indications of a narrowing of the midlatitude jet stream in winter in response to climate change, rather than a latitudinal shift.Cattiaux et al. (2016) examined a 24-member CMIP5 RCP8.5 ensemble and Peings et al. (2017) ran the Community Earth System Model Large Ensemble assuming the high emission scenario RCP8.5.Tamarin-Brodsky and Kaspi (2017), analyzing 20 CMIP5 GCMs, also found a poleward shift of the storm tracks in winter, consisting of a simultaneous shift of the genesis regions and an extension of the trajectories before reaching maximum intensity.The trajectory extension affects regions along the northeast ocean border where the storm tracks are strengthened in the longitudinal as well as poleward direction.This finding supported Zappa et al. (2013b), who, analyzing the North Atlantic storm tracks in 19 GCMs, also described an eastward extension of the North Atlantic storm-track region into continental Europe in December-February in the twenty-first century.
In this paper, we describe a new approach to estimate cyclone activity from climate projections.The proposed method is based on empirical-statistical downscaling (ESD) and is adapted to capture features of the North Atlantic storm tracks.The projected cyclone statistics are derived from (seasonally averaged) monthly mean GCM output rather than daily or subdaily GCM output, which is typically required to estimate cyclone activity, for example, in terms of variance and covariance statistics, or identifying and tracking individual cyclones.This makes the method more computationally efficient and therefore suitable for analyzing large ensembles of models (e.g., more than 100 model simulations).Previous studies using subdaily GCM data to describe the storm tracks have been limited to analysis of smaller subsets of GCMs (e.g., Zappa et al. 2013a,b;Tamarin-Brodsky and Kaspi 2017), which can limit the range of possible outcomes and may skew the ensemble mean (Benestad et al. 2017a).The method has the additional advantage of emphasizing the large-scale features of the storm tracks and circumvents the potential problems associated with the low resolution of GCMs that could influence the representation of individual cyclones.

Method and data
An empirical-statistical downscaling method was adapted with the aim of estimating future storm-track statistics from seasonal mean rather than daily or subdaily GCM output.The analysis has the advantage of being computationally efficient, which makes it possible to apply to large GCM ensembles.The complete analysis can be summarized in the following steps: 1) Identify cyclone trajectories from 6-h reanalysis sea level pressure (SLP) fields and aggregate these to seasonal cyclone densities referred to as Y.

a. Cyclone identification and tracking
The analysis included several steps, described above, starting with the generation of historic cyclone density fields over the North Atlantic.The generation of cyclone density fields involves three steps: first the identification of a cyclone, then the tracking of the cyclone, and finally the aggregation of individual cyclone trajectories into the monthly mean cyclone density fields.The procedures were conducted using two algorithms within the freely available ''esd'' package (Benestad et al. 2015), which is a collection of climate-analysis and empirical-statistical downscaling tools for the R environment that were developed at the Norwegian Meteorological Institute.
The cyclone identification started with applying an updated version of the calculus-based cyclone identification algorithm (CCI; Benestad and Chen 2006) on gridded 6-hourly instantaneous SLP data [ERA-Interim 0.758 3 0.758 grid data (T255 spectral)] for the time period 1979-2016.The CCI algorithm identifies local SLP minima by representing the meridional and zonal pressure profiles as truncated Fourier series and finding the points at which the first derivative was zero and the second derivative was positive in both directions.
The cyclone tracking was conducted by searching for likely trajectories by connecting cyclones in three consecutive time steps and minimizing the total displacement as well as the change in direction and displacement.To avoid spurious SLP minima being identified as cyclones, shallow and short trajectories were excluded that had 1) a lifetime shorter than 8 time steps (2 days), 2) a displacement lower than 500 km from cyclogenesis to cyclolysis, and 3) no times with central SLP lower than 980 hPa.This newly developed tracking scheme is available as the ''track'' function in the esd package.
Based on the historical cyclone trajectories identified as described above, the monthly mean cyclone density was then calculated as the number of cyclones per month within a 500-km distance for each grid point in a field with same spatial resolution as the ERA-Interim data, 0.758 3 0.758.The selected distance was chosen based on the average cyclone radius, which was approximately 500 km for the historical storms (cyclone radius is estimated within the CCI algorithm).

b. Empirical-statistical projection of cyclone characteristics 1) EMPIRICAL-STATISTICAL MODELING
To project the cyclone activity based on monthly mean output from CMIP5 GCMs, we used a similar approach to what has been employed in ESD of temperature and precipitation (Benestad et al. 2016).The basic steps of the ESD method include 1) using multiple linear regression to find a statistical connection between a predictand and a predictor and 2) applying the regression model to GCM output of the predictor to obtain projections for the future.Here, we used as predictand the three leading PCs of the (seasonally averaged) North Atlantic cyclone density (domain: 508W-408E; 508-808N).Using the PCs of cyclone density as predictands had the advantage of making the procedure fast and computationally efficient as well as emphasizing the large-scale patterns of the storm tracks.
The predictors used for calibrating the ESD models were (seasonally or annually averaged) monthly mean SLP or Z 500 from NCEP-NCAR over the North Atlantic Ocean (domain: 808W-808E; 458-858N).We used the NCEP-NCAR reanalysis because it provided a longer data record than most other reanalyses, which is beneficial to the model calibration in ESD.To combine and identify similar spatial patterns in the reanalysis and GCM output, common EOFs were used to represent the predictors.To avoid overfitting only the 10 leading EOFs of the predictors were used in the downscaling exercise, and the number of predictor EOFs were further reduced through a stepwise screening based on the Akaike information criterion (AIC; Akaike 1974).The linear regression model can be described as where ŷs j are the estimated values of the predictand PC y s j for PCs j 2 [1, 2, 3] and season s, N is the number of PCs of the predictor (the value for N was less than 10 and depended on the stepwise screening based on AIC), xi s j are the PCs of the predictor, and c0 s j and ci s j are the regression coefficients.Here, both predictands and predictors involved PCs, the former being the PCs representing gridded cyclone densities and the latter being PCs representing SLP or Z 500 from the combined data matrix with reanalysis and GCM results.Separate regression models were fit for each season and for the annual mean.Both the predictand and predictor PCs (only the reanalysis part) were detrended by subtracting the linear trend before linear regression to avoid that existing trends influence the calibration.

2) VALIDATION OF THE DOWNSCALING MODELS
The skill of the ESD models was evaluated by a fivefold cross validation (Gutiérrez et al. 2013), meaning that a part (one fifth) of the predictand and predictor datasets was withheld from the multiple linear regression, so that an independent prediction could be made by using the regression model with the excluded part of the predictor.The procedure was repeated five times so that an independent estimate of the predictand PCs was obtained for the full length of the data.As a skill score, we used the Pearson's correlation coefficient of the PCs obtained through cross validation and the corresponding original cyclone density PCs.

c. Global climate simulations
Ensembles of model projections of SLP and Z 500 from CMIP5 (Taylor et al. 2012) were used as predictors in the ESD models to project future changes in the midlatitude cyclone densities.We considered two possible climate futures (RCP): the intermediate scenario RCP4.5 in which emissions peak around 2040 and then decline, and the high emission scenario RCP8.5, which represents business as usual.The ensemble sizes were 108 and 81, respectively, although several model simulations were carried out with the same GCM.

a. Cyclone density
The cyclone density of the North Atlantic storm tracks varied throughout the year, and the highest number of cyclones occurred in the winter season (Fig. 1).The winter cyclone density map (Fig. 1a) showed a high frequency of low pressure systems over the ocean east of Canada and south of Greenland and Iceland.The storm tracks were relatively strong and narrow over the ocean (Fig. 1a), but approaching Europe the cyclone density diminished-indicating fewer, not necessarily weaker storms-and the storm tracks widened slightly, covering all of Scandinavia and parts of central Europe with the highest cyclone density in the Barents Sea.In autumn and spring, there were fewer strong storms (Figs.1b,d), and the highest concentration occurred south west of Iceland.There were very few cyclones identified in the summer season (Fig. 1c).
Figure 2 shows the rate of change in the cyclone density over 1979-2016 based on a linear regression against time.The winter trend map (Fig. 2a) indicated a strengthening and southward shift of the storm tracks, while the changes autumn could be interpreted as a decrease in storm activity and a northward shift (Fig. 2d), and in spring (Fig. 2d) an increase in cyclone activity along the center of the storm tracks.However, a test of statistical significance (Student's t test) showed that the changes in cyclone density could not be distinguished from natural variability.The trend analysis indicated small regions with trends significant on the 5% confidence level (p , 0.05; empty circles in Figs.2a,d,e), a few grid points with trends significant at the 1% level (p , 0.01; filled circles), and no trends significant at the 0.1% level (p , 0.001).Given the size of the domain (6363 grid points), these results are within the expected range of trends under an unchanging climate (Wilks 2006).
The seasonal mean cyclone densities were validated against independent data through a correlation analysis with rain gauge data [see online supplemental material (SM) Fig. S1].The main results of this validation suggested a skillful reproduction of past storm tracks, having high correlations with precipitation along the coastal regions of Norway and mainly during winter when the midlatitude storms are more prominent.The evaluation also indicated lower correlations in the interior and during summer corresponding to regions where precipitation is more from convective processes (Benestad et al. 2017b).The rain gauge data used in the evaluation of the seasonal storm tracks were taken from the daily European Climate Assessment precipitation data (Klein Tank et al. 2002) and included stations in both Norway and Sweden over 1979-2012.

b. EOF analysis of the cyclone density
EOF analysis of the cyclone density was performed and the three leading PCs used as predictands in the empirical-statistical downscaling procedure.The EOF analysis was also used to investigate the spatiotemporal structure of the cyclone density (results for the annual mean results are displayed in Fig. 3 and seasonal analyses are included in the SM and Figs.S2-S5).
The first EOF pattern of the annual mean cyclone density (Fig. 3c) resembled a dipole with its positive center of action along the center of the storm tracks in the North Atlantic Ocean (Fig. 1e), and a negative center over Great Britain.Similar leading patterns were identified in all seasons (see the SM).A positive correlation was found between the first PC (PC1) of the cyclone density and the NAO index (Barnston and Livezey 1987;van den Dool et al. 2000;Chen and van den Dool 2003) in the annual mean, winter, spring, and autumn (correlation coefficient r 5 0.69, 0.68, 0.60, and 0.44, respectively; all statistically significant with p , 10 24 ).Like the NAO index, the positive phase of PC1 is associated with increased cyclone activity along the North Atlantic storm tracks.
EOF patterns of higher order than 1 were not as easily interpreted and varied depending on season, but in essence they can be described as east-west-or northsouth-oriented dipole or tripole patterns that influence the extension and tilt of the storm tracks.
Based on the eigenvalues (Fig. 3b), the three leading PCs accounted for 59%-70% of the variance in cyclone density, depending on the season (Table 1).An additional analysis was included in the SM (and Figs.S6-S11), comparing the original cyclone density and the corresponding field reproduced from a varying number of EOFs and PCs.This comparison looked at the ratio of the standard deviation of the original and reproduced field, and the correlation between the two.The results indicated that the first three patterns of the cyclone density reproduced the cyclone activity in the North Atlantic well, but including 5 or 10 EOFs improved the representation and widened the region that was well represented.This was true for all seasons except summer, in which the correlation was low in the whole domain, even when the standard deviation ratio was close to one.These results indicated that the summer storm tracks may not be readily described by large-scale patterns and hence may not be suitable for analyses such as EOF or ESD.

c. Validation of the empirical-statistical downscaling
As described in section 2, empirical-statistical models of cyclone density were calibrated using linear regression between the PCs of the predictand (the seasonal mean cyclone density) and predictor (seasonal mean SLP or Z 500 ).A fivefold cross validation was then applied to evaluate the skill of the ESD models in terms of predicting independent data.
The cross validation showed that the empiricalstatistical models had considerable skill in predicting the first three PCs of the cyclone density in the northern North Atlantic in all seasons except summer (Table 2).The results of the cross validation were similar for the regression models using SLP and Z 500 as predictor, with SLP scoring only slightly higher correlation values.The cross-validation correlation coefficients were around 0.8-0.9 for the first PC, and between 0.4 and 0.7 for the second and third PCs, with slight differences depending on the season and predictor.The ESD models for the winter season had the highest predictive skill.However, in the summer season when cyclones are less frequent and not well described by large-scale pattern (as discussed in section 3b), the empirical-statistical models had no predictive skill and the correlation scores were small or negative.
To further compare the predicted and original cyclone density, we calculated the cyclone density based on the first three cross-validation PCs and compared with the original cyclone density field.The cross validation showed very similar results for both predictors (correlation maps for the analysis using Z 500 as predictor are shown in Fig. 4 while the results using SLP as predictor are included in the SM: Fig. S12).On the one hand, the correlation between the predicted and observed cyclone density was high along the center of the North Atlantic storm tracks for all seasons except summer.On the other hand, the empirical-statistical models had difficulties predicting variations in areas where storms are less frequent, like over southern Scandinavia, central Europe, and Greenland.For the annual mean and the winter season, correlation values were as high as 0.7-0.9along large areas of the storm tracks.The values were slightly lower in spring and autumn than in winter, but, the correlation map in summer indicated that the regression models have no predictive skill.
Additional cross validation showed that including higher-order PCs of the predictand (cyclone density) in the ESD models did not contribute to a more skilled or complete prediction than ESD models using the three leading predictand PCs (see cross-validation correlation maps in the SM: Figs.S12-S15).The reasons are twofold: 1) higher-order PCs carry less information than the leading ones (Table 1 and Fig. 3) and 2) higher-order PCs are poorly represented by the ESD models because they represent finer-scale variability rather than large-scale climate patterns.

d. Projected changes
Estimates of cyclone density were calculated for the future by applying the empirical-statistical regression  2. Cross-validation correlation scores of the empiricalstatistical regression models for cyclone density.The correlation scores r PC1 , r PC2 , and r PC3 are the Pearson's correlation coefficient between the first three PCs of the cyclone density fields and the corresponding PCs calculated with the empirical-statistical regression models.The scores are based on a cross-validation method that repeatedly excludes one part of the cyclone density data, thus allowing independent comparison with data that have not been used to tune the regression models.The table shows results for all seasons and for regression models using sea level pressure or 500-hPa geopotential height as predictors.

SLP
Z The results from the multimodel ensemble are summarized here in terms of the 5th, 25th, 50th (median), 75th, and 95th percentile of change from the present day  to the far future (2070-99) within the full ensemble (Figs. 5 and 6 show results assuming emission scenario RCP8.5).
Using SLP as predictor (Fig. 5), a majority of projections indicated an increase in cyclone density in the northern part of the North Atlantic and in the Barents region.This change may not be statistically significant for all ensemble members.The results were similar in all seasons, but stronger changes were seen in the winter season compared to spring and autumn.In winter the strongest increase was seen farther southwest in the storm tracks.
Projected annual means of the cyclone density based on the geopotential height at 500 hPa (Z 500 as predictor) showed changes somewhat similar to the projections using SLP as predictor, with increasing cyclone density over the northern part of the domain (Figs.6a-e).In autumn and spring (Figs.6k-t), the projections based on Z 500 indicated an increase in cyclone density in the north part of the storm tracks and decreasing in the southeastern part of the domain, with a larger spread within the ensemble and more pronounced changes compared to the annual mean.In winter (Figs.6f-j), a majority of ensemble members showed a southeastward shift of the storm tracks, with a decreasing cyclone density north of Iceland and in the Barents region and an increase over the southern part of the domain, which differs considerably from the projections based on SLP as predictor.

e. Analysis of the CMIP5 ensembles of SLP and Z 500
To address the difference between the predictors, we investigated the CMIP5 ensembles of SLP and Z 500 .The CMIP5 projections indicated that Z 500 increases during the twenty-first century in the North Atlantic region for all GCMs, seasons and assuming both RCP8.5 and RCP4.5 (not shown), although stronger changes are associated with the high emission scenario.The zonal mean trends for the winter season assuming RCP8.5 are displayed in Fig. 7b.In the winter season (December-February), the increase in Z 500 is stronger in the polar region, north of about 708N, which indicates a decrease in the meridional (north-south) gradient, as the mean Z 500 tends to be lower at high latitudes.
For SLP, the CMIP5 models indicate no consistent trend in SLP in the North Atlantic domain.Some GCMs indicate a negative trend in the polar region, in the winter (Fig. 7a) as well as in spring and autumn (not shown), which constitutes a strengthening of the mean meridional SLP gradient, while others show no notable trend during the twenty-first century for either emission scenario.

f. Test of stationarity
An important assumption for the empirical-statistical method presented here is that the relationship between predictand (PCs of the cyclone density) and predictor (PCs of SLP or Z 500 ) is stationary (i.e., that it does not change over time).To test this assumption, the cyclone density was derived from 6-hourly SLP data for one of the CMIP5 simulations (HadGEM2-ES with emission scenario RCP8.5) for the period 1949-2100.The empirical-statistical projection method was then applied to this cyclone density record using the seasonal mean SLP or Z 500 from the same GCM as predictor with three different calibration periods (1979-2010, 1949-2010, and 2069-2100).
The test indicated that Z 500 was not a robust predictor.The empirical-statistical model changed depending on the calibration period, revealing that the relationship between the cyclone density patterns and the Z 500 patterns was not stationary (see the analysis in the SM and Figs.S16-S18).Depending on the calibration period, the projected changes in cyclone density also changed drastically and the projected trends differed from the calculated trends.Using a longer calibration period  did not result in a more skillful prediction.FIG. 6.As in Fig. 5, but based on statistical models using Z 500 as predictor.
SLP appeared to be a more robust predictor in the sense that the models trained on SLP showed less sensitivity to the choice of calibration period than the models trained on Z 500 .The model trained on HadGEM2-ES's seasonal mean SLP was skillful in predicting historical to far future changes in HadGEM2-ES's cyclone density in spring and autumn.However, for the winter season, the projected changes displayed a different pattern from the calculated changes.

Discussion
Comparison with independent data (cross validation) was an important step to assure that the ESD models of cyclone density had a basis in reality.For example, based on the poor cross-validation correlation, the summer season could be excluded from further analysis because of the lack of predictive power.However, the cross validation did not give any indication that one choice of predictor, SLP or Z 500 , would be a better choice.Nevertheless, when applied to the CMIP5 ensemble the two predictors gave different results in projected cyclone density for the high latitudes during the winter season (DJF; Figs.5f-j and 6f-j).
The most plausible explanation for the diverging projections for the cyclone density in winter is that the spatiotemporal patterns of SLP and Z 500 evolve differently with a global warming (as seen in the analysis of CMIP5 SLP and Z 500 ) and that the predictorpredictand relationship (i.e., the empirical-statistical model) between one or both of them and the cyclone density change depending on the calibration period, which would violate the assumption of stationarity that the ESD method is based upon (Hewitson et al. 2014).This suspicion was confirmed by downscaling the cyclone density calculated from 6-hourly HadGEM2-ES SLP data.A stationarity test indicated that the empiricalstatistical models of cyclone density were more sensitive to the calibration period when using Z 500 as predictor than SLP.In winter, both predictors gave spurious trends in the projected cyclone density.The results highlight the importance of testing the assumption of stationarity in empirical-statistical methods by utilizing the long simulations of GCMs and RCMs.
It has been suggested that when choosing an ESD approach, the choice of predictor should include the direct radiative influence of greenhouse gases and all sources of climate change (Hewitson et al. 2014).While SLP represents the dynamic response to external climate change forcings (i.e., changes in the large-scale atmospheric circulation), Z 500 is sensitive to the changes in temperature of the atmosphere between the surface and the 500-hPa level and thus more closely related to the thermodynamic (direct radiative effect on the base temperature) climate change signal (Christidis and Stott 2015).By this logic, it would make sense to choose Z 500 over SLP as a predictor, assuming that there was a plausible physical reason why the projected change in Z 500 should translate to a change in the cyclone density and that the predictand-predictor relationship is stationary.However, as illustrated by the stationarity test, Z 500 is not a robust predictor of future cyclone density in a changing climate.
The analysis of the CMIP5 projections of Z 500 and SLP provided some clues about the reason for the discrepancy.The mean Z 500 increased in the North Atlantic region for all CMIP5 members and all seasons during the twenty-first century assuming RCP4.5 and RCP8.5 (Fig. 7b shows results for winter assuming emission scenario RCP8.5).However, given the large spread in the response in the projected cyclone density (Fig. 6), with different results in different seasons and a large spread between CMIP5 members, it is probable that it is the change in the Z 500 spatial pattern rather than background mean that drives the changes in the statistical projections of cyclone density.In winter, the projected weakening and southward shift of the storm tracks assuming RCP8.5 (Figs. 6f-j) may be related to the projected decrease in meridional gradient in Z 500 seen in many of the CMIP5 members (Fig. 7b).On the other hand, using SLP as predictor, some of the CMIP5 members indicated a strengthening of the projected winter storm tracks (see the 75th and 95th percentile of change in Figs.5i,j), and this could be linked to the strengthening of the SLP gradient over the twenty-first century seen in some of the GCMs (Fig. 7a).
The changes in CMIP5 Z 500 and SLP assuming RCP8.5 described here are consistent with previous studies, such as Cattiaux et al. (2013) who found that Z 500 indicated a negative trend of the winter northern annual mode (NAM) under anthropogenic climate change, while the SLP-based NAM did not change consistently but varied between GCMs.The negative phase of NAM is associated with suppressed subpolar westerlies and storm activity at high-and midlatitudes (Thompson and Wallace 2001).The projected changes in the winter cyclone density using Z 500 as predictor-which indicated a weakened and southeastward shifting storm track-are thus consistent with a negative trend in NAM and weakening of the meridional gradient of Z 500 .The nonstationary link between cyclone density and Z 500 suggests that the projected negative trend in NAM based on Z 500 may not actually translate to a weakening of the storm tracks.
A comparison between the projected cyclone densities for the autumn, winter, and spring seasons can give an indication about the quality of the downscaled results when there are no known reason that they should evolve differently.In this case, one reason for why the projections for winter could differ from other seasons is that the Arctic seas are covered by sea ice over a larger area.Furthermore, the downscaling could be repeated with different predictor domains to explore how robust the results are and whether the choice of domain size was too large (Benestad 2001).
One possible problem with the method described in this paper is that only the first three PCs of the cyclone density were downscaled, which may not be enough to capture the changes in cyclone density in all regions.This could be solved by representing the storm-track statistics in a different way, for example, calculating the cyclone density for a coarse grid and downscale each grid cell rather than the PCs.This would be more computationally expensive but may give more robust projections.
The main reason for developing the adapted ESD method for cyclone density was the computational efficiency of the method, which allowed analysis of large ensembles of (seasonally averaged) monthly mean GCM model output.Previous downscaling analysis suggest that model differences are often less pronounced than internal variability on regional scales (Benestad et al. 2016), and it has previously been demonstrated that one single GCM is able to produce widely different regional and local projections (Deser et al. 2012).This implies that even if you would have a perfect GCM, you would need a large ensemble of GCM runs to represent the internal variability of the climate system.Using a subset of GCMs can, if the subensemble is too small and is not selected carefully, result in a skewed or too-narrow view of the consequences of climate change.Here, we used so-called ensembles of opportunity, that is, all available GCMs from the CMIP5 experiment.This approach provided a large sample size (many GCM projections), but a more carefully constructed sampling protocol might have better represented the full range of possible outcomes.
The results of the statistical projections may potentially be sensitive to the cyclone identification and tracking method and how the cyclone density was calculated.For example, the exclusion of shallow (SLP .980 hPa at all times) cyclones may leave out many of the summer storms, and the exclusion of short trajectories (displacement , 500 km) could be problematic for detection of the Iceland lows, which are sometimes stationary.Here, the cyclone density was calculated as the monthly frequency within 500 km of each grid cell (which subsequently were aggregated to seasonal and annual means), and the distance to the cyclones taken into consideration influenced the smoothness of the cyclone density field.The evaluation of the seasonal cyclone densities through correlation with independent rain gauge data suggested a high degree of realism, however, because the pattern of correlation matched well the expectations of precipitation predominately associated with midlatitudes over the coastal areas and less over the interior parts where local convection is more predominant.The regression analysis involved in the ESD also suggested high skill through a good match between predictand and predictors.
The spurious trends in winter cyclone density found in the GCM over the Arctic could be related to the choice of predictor domain.When the predictor covers an extensive area characterized by a dipole pattern, anomalies of opposite polarity may carry different weight and a change in their relative strengths can introduce such spurious trends (Benestad 2002).Other predictands and predictors could be used with the statistical projection method described here.Various features of the storm tracks, such as the strength and the radius of cyclones, may be associated with other large-scale climate features.However, it is important that there is a valid physical reason for assuming a statistical connection, that the statistical relationship is independently validated, that the statistical relationship does not change in the future, and that the predictor represents all relevant sources of climate change forcings.The choices used here, using both SLP and Z 500 as predictors, were based on the expected physical dependency between these fields and the storm activity, since storms are associated with low pressure systems and a result of baroclinic instability (Gill 1982).

Conclusions
Our main finding was that, although the cyclone density can be predicted from large-scale predictors such as SLP and Z 500 , the choice of predictor affects the future projections because SLP is less sensitive to increased temperatures than Z 500 is.Additional analysis indicated that the relationship between predictand and predictor was not stationary when using Z 500 as predictor.For future studies, we suggest using pseudoreality datasets such as GCM or RCM simulations to test stationarity in the predictor-predictand link for empirical-statistical downscaling analyses.There is a large spread within the cyclone density projections based on CMIP5 SLP simulations; however, the multimodel ensemble median indicates an increase in the storm-track activity in winter, and a poleward shift during autumn and spring.

FIG. 1 .
FIG. 1. Seasonal average cyclone density in the North Atlantic for the (a) winter (DJF), (b) spring [March-May (MAM)], (c) summer [June-August (JJA)], (d) autumn [September-November (SON)], and (e) the whole year.The cyclone density in each grid box was calculated as the number of cyclones per month within a 500-km radius, based on the cyclone trajectories identified and tracked from ERA-Interim reanalysis SLP data for the period 1979-2016 as described in section 2.

FIG. 3 .
FIG. 3. EOF analysis of the annual mean cyclone density, calculated from ERA-Interim SLP data for the period 1979-2016 as described in section 2: (a) the three leading PCs that represent the temporal variation in (c)-(e) the corresponding spatial EOF patterns; (b) the variance associated with each pattern.The first three EOFs explain 67% of the total variance of the annual mean cyclone density.

FIG. 4 .
FIG. 4. Cross validation of the statistical models of cyclone density, with Z 500 as predictor: the correlation between the predicted and independent original values for (a) winter (DJF), (b) spring (MAM), (c) summer (JJA), (d) autumn (SON), and (e) the annual mean.

FIG. 5 .
FIG. 5. Statistical projections of the change in cyclone density from the present (1981-2010) to the far future (2070-99) assuming the RCP8.5 scenario.The projections are based on a statistical model applied to the CMIP5 ensemble, using SLP as predictor.The maps show the (a) 5th percentile, (b) 25th percentile, (c) median, (d) 75th percentile, and (e) 95th percentile of change in the ensemble.The projected absolute change is shown in the same units as in the maps of present-day cyclone density (Fig. 1), i.e., the number of cyclones per month within a 500-km radius.

FIG. 7 .
FIG. 7. Zonal mean of the winter (DJF) trends of (a) SLP and (b) Z 500 for the CMIP5 model projections in the North Atlantic domain (608W-508E; 458-858N) during the period 2000-99 assuming the high emission scenario RCP 8.5.Note that the spread toward the high latitudes may be an artifact resulting from lower degrees of freedom (Benestad 2005).

TABLE 1 .
Accumulated variance (%) accounted for by the leading PCs of the annual and seasonal mean cyclone density.Only the three leading PCs were used for the predictand in the downscaling and projection here, accounting for 59%-70% of the total variance in the cyclone density.