Polar low variability and future projections for the Nordic and Barents Seas

Polar lows are intense mesoscale cyclones occurring during winter over open sea areas in certain polar sub‐regions. Due to their small size, they are not explicitly represented in present global climate models or Earth system models. In this study 18 members of the CESM Large Ensemble were dynamically downscaled to ∼12 km horizontal mesh width using the quasi‐hydrostatic ALARO model within the HARMONIE script system in climate mode (HCLIM‐ALARO). The domain covers the Nordic and Barents Seas. One historical and two future time‐periods were selected. For validation, the ERA‐Interim reanalysis was also downscaled. A cyclone‐tracking algorithm was used to identify tracks of individual polar lows. Their frequency of occurrence, lifetime, and maximum relative vorticity were estimated. Relative to ERA‐Interim, the historical frequency of occurrence of polar lows was slightly overestimated in the Nordic Seas and underestimated in the Barents Sea, which is likely due to positive biases in sea‐surface temperature and sea‐ice concentration. For future climate projections, the regions of polar low genesis are diagnosed to move northwards in accordance with the sea‐ice retreat. In the Nordic Seas, the number of polar lows decreases at the beginning of the season, while there is an increase in March. In the Barents Sea, a February–April increase in the occurrence of polar lows is seen.


INTRODUCTION
Polar lows are small, intense mesoscale cyclones in the polar areas. Although they have been studied since the late 1960s, we refer to Rasmussen (2011) for a detailed description. Polar lows form in the winter season during spells when the temperature difference between the atmosphere and the sea surface is large. In the Nordic Seas, duty forecasters at MET Norway have diagnosed that there are about twelve polar lows every year on average (Noer et al., 2011). Although considerable improvements in numerical weather prediction (Aspelien et al., 2011;Kristiansen et al., 2011;Jung et al., 2016;Müller et al., 2017) and observational data cover have reduced harmful impacts on society and individuals, polar lows still pose a considerable risk to offshore activities, shipping, fisheries, and coastal communities (West and Hovelsrud, 2010;Gudmestad and Karunakaran, 2012;Jung et al., 2016). Additionally, with increasingly ice-free conditions in the Arctic, economic potential for offshore industries such as shipping, fishing, oil and gas and tourism are expected to increase in the future. It is therefore crucial to have knowledge of future projections of polar lows.
Because of their small size, polar lows are not always detected in reanalyses (Zappa et al., 2014). Neither are they resolved by General Circulation Models (GCMs) or Earth System Models (ESMs), which typically operate on ∼ 100 km horizontal mesh width. However, a downscaling method is an option to produce realistic smaller-scale features required to identify polar lows. Zahn and von Storch (2010) used the regional climate model (RCM) CLM to downscale several global models from the Coupled Model Intercomparison Project Phase 3 (CMIP3) to 50 km horizontal mesh width. They found a future decline in frequency of North Atlantic polar lows as well as a northward shift (∼ 2 • latitude) of the genesis region. Romero and Emanuel (2017) applied a statistical method to the newer CMIP5 ensemble and similarly found a future decline (10-15%) in North Atlantic polar lows, with a strong decrease in the west and south and a small increase in the east.
In this paper we apply an atmospheric RCM with ∼ 12 km horizontal mesh width. This provides approximately four times better resolution than the ∼ 50 km mesh width used by Zahn and von Storch (2010), thus we aim to provide a considerably better representation of polar lows, which can range in scale from less than one hundred to several hundred kilometres (Rasmussen, 2011).
As a coarse-scale proxy to study conditions favorable for polar low formation, Kolstad and Bracegirdle (2008) defined an index of low static stability below 700 hPa to identify cases of Marine Cold Air Outbreak (MCAO). Applied to the future scenarios of the CMIP3 ensemble, they found a projected future decrease in strength of MCAOs over open sea water due to the faster warming in the atmosphere than the ocean, while they also found an increase over the marginal ice zones as a consequence of the ice retreating northwards. However they also noted that many models had large biases in both atmospheric temperature and sea-ice concentration.
While proxies can be a useful way to extract data from a coarser model using a statistical downscaling method, there is no guarantee that the links will remain the same under future climate change. For example, Mallet et al. (2017) found that the statistical link between static stability and weather regimes, which is robust over past decades, weakens in the future due to reduced variability associated with sea-ice retreat. This presents challenges for traditional statistical downscaling methods. We instead apply a more rigorous, albeit computationally demanding, strategy using dynamical downscaling together with a cyclone-tracking algorithm to identify and track polar lows individually.
Any single simulation from a global climate model is intended to represent a combination of responses to an external forcing (e.g. increase of anthropogenic greenhouse gas concentration) and internal, chaotic variability which includes switches between apparent modes of large-scale circulation such as the North Atlantic Oscillation. Since polar lows are relatively rare events, multiple model realizations are needed to achieve a statistically robust result in a situation with non-stationary climate statistics. The Community Earth System Model (CESM) Large Ensemble (CESM-LE; Kay et al., 2015) is one of the most recently produced large ensemble of climate model scenario calculations aimed at separating changes due to external forcing from internal variability. Its horizontal resolution is also high for a large ensemble (1.25 • longitude × 0.9 • latitude mesh width). Created by applying small perturbations to the temperature field initially, CESM-LE provides 40 ensemble members for the years 1920 to 2100 following the RCP 8.5 scenario (Riahi et al., 2011). To our knowledge there have not yet been any attempts to separate external forcing and internal variability in relation to polar lows. Although previous studies (e.g. Zahn and von Storch, 2010;Romero and Emanuel, 2017) have used multiple ensemble members, there were only single or a few members of each model.
We applied dynamical downscaling using a regional climate model in order to identify and track polar lows. By using 18 members of the CESM-LE, we aim at separating responses to external forcing from internal variability, as well as comparing with the ERA-Interim (ERA-I; Dee et al., 2011) reanalysis downscaled using the same model set-up.
In this paper, we intend to answer the following questions: (1) How well are polar lows represented in our specific GCM+RCM model chain compared to downscaled reanalysis? (2) How will polar lows change in the future? (with particular emphasis on frequency, life times and variability).
(3) What is the magnitude of natural variability compared to the forced response?
In Section 2 we describe the methods, including details about the data used, dynamical downscaling set-up and tracking algorithm. The subsequent sections present the results (Section 3), discussion (Section 4) and conclusions (Section 5).

Data
We retrieved CESM-LE data from the Climate Data Gateway at NCAR (https://www.earthsystemgrid.org/dataset/ ucar.cgd.ccsm4.CESM_CAM5_ BGC_LE.html; accessed 31 July 2019). The following variables were used: 6-hourly data on all model levels -temperature ( ), zonal wind ( ), meridional wind ( ), specific humidity ( ) -and 6-hourly data at the surface -surface pressure (PS) and surface skin temperature (TS). For sea surface temperature (SST) and sea-ice Atmospheric physics ALARO (Gerard et al., 2009) Land surface model SURFEX (Masson et al., 2013)   concentration (SIC), 6-hourly data were not available, so daily fields were used instead. A total of 18 ensemble members were used. We started out with members 1-15, and members 17, 19 and 25 were added based on the decision to include members with both low and high sea-ice concentrations ( Figure S1).

Dynamical downscaling
We used the HARMONIE script system in climate mode (HCLIM; Lindstedt et al., 2015). The model set-up is detailed in Table 1. The HARMONIE system allows for a choice of different atmospheric physics packages available in the shared ALADIN-HIRLAM (Aire Limitée Adaptation dynamique Développement InterNational-High-Resolution Limited-Area Model) numerical weather prediction (NWP) system (ALADIN, ALARO and HARMONIE-AROME). Our choice, ALARO, is a model designed for grey-zone dynamics, providing consistent results at horizontal resolutions with mesh width ranging from a few tens of kilometres down to less than 4 km (Gerard et al., 2009).
The HARMONIE system uses the surface modelling platform SURFEX to parametrize surface processes. Within SURFEX sea-ice is modelled by the Simple Ice scheme (SICE; . In its standard set-up, SICE keeps sea-ice thickness set to a constant value, which may be acceptable in numerical weather prediction mode (with typical forecast ranges from a few days to a week), but it leads to large biases in the energy fluxes over longer time-scales. Therefore an ice mass balance (prognostic sea-ice thickness) parametrization has been implemented (Batrak, personal communication, 2018). Sea-ice concentration was updated from the driving, global data every 6 hr. The sea-ice thickness initial state was set to 1 m. Two reference datasets were used for comparison: Global Ice-Ocean Modeling and Assimilation System (GIOMAS; Zhang and Rothrock, 2003) and TOPAZ4 (Sakov et al., 2012). As shown in Figure 1, our simulated sea-ice thickness is between the two most of the time, and an improvement compared to using constant thickness. For snow-on-ice, our simulation has a positive bias compared to TOPAZ4, but there is likely also a problem with TOPAZ4 being saturated at ∼ 0.4 m snow depth. The model domain ( Figure 1) covers the Norwegian, Barents and Kara Seas with 192 × 360 grid points on ∼ 12 km horizontal mesh width and with 65 levels in the vertical. The downscaled regional output was stored at 3-hourly time resolution, in order to improve cyclone tracking ability over the 6-hourly global data. The choice of time periods was limited by the availability of the 6-hourly data from CESM-LE; 15 years in the historical period and two 10-year periods in the future scenario periods. For the run driven by ERA-Interim, a 20-year period was used. For all simulations, soil variables were initialized from the climatological September mean calculated from ERA-Interim 1990-2005. While this may not be the best approach for the future periods, we considered it sufficient since our main interest is ocean areas, where SST, SIC as well as skin temperature over ice (weighted mean of SST and ice surface temperature) are all initialized from CESM. Furthermore, we initialize in September, which is the climatological sea-ice minimum.

Biases
An overview of the biases of CESM-LE in this region is available in Landgren et al. (2019). In particular there is a positive SIC bias, strongest in the Barents Sea, as well as a positive SST bias east of Iceland. In our simulations, HCLIM interpolates SIC and SST directly from CESM-LE, so our simulations will inherit any biases in these fields. For polar lows, a positive bias in SIC means that the genesis region may be slightly displaced geographically. In particular we may expect somewhat fewer polar lows near the eastern coast of Greenland and in the eastern Barents Sea compared to reanalysis. It is not obvious if or how this affects polar low frequency, intensity or lifetimes. Differences after downscaling are shown in Figure 2. As expected, the largest changes are in mountainous and coastal regions. Along the Norwegian west coast, the downscaled data give a higher and more realistic precipitation than the coarser ERA-Interim. Over open ocean the changes are small, mainly because SST and SIC are passed through the downscaling model unchanged.
When the simulations were performed, the fifth-generation ECMWF reanalysis, ERA5, was not yet available for the full simulation period, which is why we used ERA-Interim. It has subsequently become available, offering ∼ 2.5 times higher horizontal resolution, which will be useful for future work.
Comparing the 18 downscaled members of the CESM-LE to downscaled ERA-Interim, the largest differences ( Figure 3) are related to differences in the SIC. The 500 bias means that CESM-LE, like many other CMIP5 models, has North Atlantic/European storm tracks which are too zonal, as discussed in (e.g.) Day et al. (2018). Biases for individual months October-May are shown in Figure S2.

Tracking
We applied a cyclone detection and tracking method used operationally for polar low forecasting at MET Norway. The detection algorithm is based on three criteria: vorticity at 925 hPa, SST-500 difference, as well as a threshold on wind speed at 10 m above sea level. In their comparison of the effectiveness of different criteria for detecting polar lows, Stoll et al. (2018) ranked the SST-500 difference third best among the 18 MCAO criteria, only marginally worse than the first two. In their study, the threshold including 90% of all polar lows was 41.4 K. Zahn and von Storch (2008) used 43 K, while MET Norway uses 40 K as threshold in operational polar low forecasting. This makes sense because a duty forecaster deals with only a handful of cases simultaneously, and a false positive can easily be manually discarded, while a false negative could lead to a missed forecast.
For the wind speed criterion we have used 15 m/s as threshold (as in Heinemann and Claud, 1997), while Zahn and von Storch (2008) used 13.9 m/s (near-gale, as in Rasmussen, 2011). For this we use a box of 132 × 132 km (11 × 11 grid cells), similar to the 2 • latitude radius used by Stoll et al. (2018). It should be noted that the wind speed threshold in the definition is based on observations, while modelled 10 m wind speed is a diagnostic field which is often subject to tuning. Modelled values of extreme winds also increase with higher resolution (e.g. McInnes et al., 2011;Moore et al., 2015) which may also support using the higher value, though a thorough investigation of the scaling of the threshold is beyond the scope of this study. (Sensitivity to the choice of wind threshold for the downscaled ERA-Interim simulation is shown in Table S1.) Regardless, Stoll et al. (2018) found that the 10 m wind speed was not very effective at identifying polar lows. As a measure of intensity we use the maximum relative vorticity (at 925 hPa) for each polar low track, but no threshold is applied to this quantity. For the spatial filter we select wavelengths between a lower and an upper threshold.
To study the sensitivity of the results to this parameter, three different intervals were chosen: (a) 150-600 km, as used in the operational forecasting at MET Norway, (b) 200-600 km, as used by Zahn and von Storch (2008;, as well as (c) 250-850 km to include slightly larger systems, noting that Rasmussen et al. (1993) and Rasmussen (2011) have the upper bound at 800 and 1,000 km, respectively.
We applied the tracking algorithm to data from our HCLIM-ALARO model set-up at 12 km, while in the operational setting the method uses data from the 2.5 km HARMONIE-AROME model. Due to the difference in horizontal resolution, some parameters were scaled accordingly, but the method is otherwise unchanged. A summary is given here, while additional internal documentation and Fortran code are available upon request.
The relative vorticity is calculated from the and fields at 925 hPa, and then transformed using a discrete cosine transform (DCT) method based on Denis et al. (2002). A forward DCT is performed twice on the vorticity field, once in each of the and dimensions. A filter matrix S is set up such that only the relevant wavelengths are kept. After multiplying the transformed vorticity with S, a backwards DCT is performed (also in the and dimensions).
Track points are defined as follows. First, a potential point is found as a local maximum in the relative vorticity field. If the difference between the value at this point and the area average is larger than a threshold, it is accepted as a track point. The area is defined as a circle of radius seven grid cells, and the threshold is set to 5× 10 −5 s −1 . After defining the track points, the temperature difference SST-500 is calculated and track points below a threshold are discarded. Track points are also discarded if no grid cells within an 11 × 11 grid cell box have 10 m wind speed above 15 m/s. The track points are grouped into tracks by considering the coordinates of individual track points in three steps which are repeated for each model output time step (in our case every 3 hr): 1. If a track has two or more points, its position in the next time step is extrapolated from the last two and a spiral search is performed around the extrapolated position. If a local vorticity maximum is found within a radius of eleven points, it is included in the track. Otherwise the track is marked as ended. 2. Next, tracks containing only one point are considered, and the nearest local vorticity maximum (up to a radius of eleven grid points) is selected as the next track point. 3. New tracks are formed by local vorticity maxima that are at least eight grid points away from other local maxima.
After all simulation time steps have been considered, only tracks containing at least four points (i.e. 9 hr from beginning to end) are kept.
Analysis was performed in two regions: (1) Nordic Seas: 20 • W-20 • E, 65-78 • N; and (2) Barents Sea: 20-54 • E, 70-78 • N. For a track to belong to a region, it was required to have at least three track points within the region. Tracks were allowed to be assigned to both regions simultaneously.
We also applied an optional criterion to check whether the track is south of the polar front or not. This was done by taking the 500 hPa wind speed in a meridional band of width 2 • longitude centred on each track point (±1 • longitude). If within this band the southernmost location above 24.8 m/s (threshold taken from Stoll et al., 2018) was north of the track point, the whole track was discarded. However, because of the size of our regional domain, the polar front may not always be well defined, and we chose to present the results (Section 3.2) without this criterion. It is included in the discussion in Section 4. Regarding the thresholds for temperature difference, both 40 and 43 K were used. The resulting trends were similar (for details see the discussion regarding Figure 7), however the lower value naturally yielded a larger number of tracks. In order to reduce the number of false positives, we decided to focus on the more strict threshold in the results.
In order to visualize the track density, filled circles of radius 100 km were drawn around each track point. By summing over all tracks and dividing by the simulation length, the unit becomes number of tracks per month. In order to avoid double-counting when circles of slow-moving polar lows overlap, each track could only expose each grid point once. Of course, multiple exposures from multiple tracks were allowed. While drawing using a constant radius does not take into account any variations in the actual horizontal extent of individual polar lows, we argue it is nevertheless useful in determining areas of high and low exposure to polar lows.

Change in extreme winds
Before applying the tracking algorithm, we wanted to see if there were any general trends in extreme winds. Figure 4 shows the change in 98%ile 10 m wind speed for the 18 downscaled members. Strong winds decrease in the Nordic Seas, while they increase in the Kara Sea. Most of the Barents Sea does not have significant changes. The decrease is in agreement with the general static stability increase, while the Kara Sea increase is likely due to the decreasing SIC.

Polar low frequency and lifetime
We calculated the sum of all identified polar low tracks to produce mean track densities for the three time periods ( Figure 5). A future increase in track density is found along the retreating sea-ice edge. There also seems to be a weakening along the Norwegian coast. Due to the rare occurrence of polar lows, no conclusions should be drawn from finer spatial structures in the maps. Due to a positive SST bias in the CESM-LE in an area east of Iceland, the tracking algorithm produced many false positives in this region (because of the SST-500 criterion incorrectly classifying other mesoscale cyclones as polar lows). Tracks originating in this area were therefore filtered out before the statistical analysis could take place. The frequency, lifetime and vorticity of tracked polar low events are presented in Figure 6. In the Nordic Seas, the polar low frequency from the downscaled CESM-LE members is slightly higher than the downscaled ERA-Interim. However, in the Barents Sea, due to the positive sea-ice bias in CESM-LE, too few polar lows are formed in September-February, while the results for March-May well resemble the downscaled ERA-Interim. The largest biases are in November-January, and likely caused by the sea ice freezing too early in the CESM-LE, as shown in Landgren et al. (2019). This can be partly confirmed by comparing the ERA-Interim results with the near-future period (2026)(2027)(2028)(2029)(2030)(2031)(2032)(2033)(2034)(2035), for which the differences in both polar low frequency and ice concentration are smaller. In the later future period (2071)(2072)(2073)(2074)(2075)(2076)(2077)(2078)(2079)(2080), the effect of the temperature change driving a static stability increase is likely the dominant mechanism behind the polar low frequency decline. Comparing the historical period with the two future scenario periods, there is a statistically significant reduction of polar low frequency in the Nordic Seas in October-January, with ensemble mean change of -36, -29, -21 and -23% and -values 3.4 × 10 −3 , 1.2×10 −4 , 5.4×10 −3 and 4.3×10 −4 , respectively. In the Barents Sea the melting ice may open up more areas for polar low development in the late season. In terms of polar low lifetimes, there is a statistically significant reduction in the Nordic Seas in November-January (ensemble mean -12, -10 and -6.6% with -values 4.2×10 −3 , 4.0×10 −3 and 4.8×10 −3 , respectively), and in December in the Barents Sea (-13.4%, =0.029). The intensity, as measured by the maximum relative vorticity along the track, sees a statistically significant decrease for the ensemble mean for three out of four months in the peak season December-March in both regions.
A sensitivity study was carried out by performing the detection and tracking with different combinations of track selection criteria. Figure 7 shows statistically significant changes for the polar low frequency of occurrence, their lifetimes and maximum vorticity, separately for the months October-May. Each row shows a different combination of track selection criteria, and the colours indicate when there is a statistically significant change from the reference period to the future periods.
Changes in ensemble median statistics from the historical to the future periods are due to changes in external forcing, while the ensemble spread represents the internal variability of the model system. Figures 6 and 7 show the changes for the ensemble as a whole (spread shown is between members), but we also calculated the change in interannual variability on a per-member basis. This was done by taking the standard deviation of the monthly values, for example, using nine October values for member 1 for years 2026-2034. Because the historical period is longer than the future period, only the first nine years were used when calculating the standard deviation, so that the sample size would be the same for each period. The standard deviations were then compared between the historical period and the future period. A -test was applied to assess the significance of the change, comparing the 18 values (one for each ensemble member) of the historical period to the future period. The results are summarized in Figure 8 to highlight the sensitivity to track selection criteria, in the same way as Figure 7. There are not as many months with significant changes as in Figure 7. However, there are significant changes in March, with increased interannual variability of frequency in both regions, and decreased interannual variability of vorticity in the Nordic Seas for the larger systems.

DISCUSSION
In order to assess the sensitivity of the detected trends to the choice of track selection criteria, we applied twelve different combinations and compared when statistically significant changes were detected. As expected, the number of detected cyclones is sensitive to the combination of track selection criteria used (Table S1). However, as can be seen in Figure 7, the projected future changes in annual cycle is not particularly sensitive. The general tendencies are towards an early season decrease in frequency in the Nordic Seas and late season increase in the Barents Sea, with an increase in March for the domain as a whole. The largest differences are due to the size criterion. There are more months with significant changes when the spatial selection includes smaller cyclones, in particular the frequency decrease in the Nordic Seas and the increase in the Barents Sea both occur in more months.
Using the polar front criterion, the number of polar lows is lower (Table S1) but there is no obvious improvement in terms of biases relative to ERA-Interim (e.g. comparing Figures 6  and S3). In terms of trends there are some differences particularly in December, but most months are unchanged.
The probability distribution of all tracked polar lows during September-May is included in Figure S4. While there is a clear reduction in the vorticity, the changes in polar low frequency and lifetimes are very small for the season as a

F I G U R E 8
As Figure 7, but for interannual variability. Calculated as the standard deviation of values for each member and period. Then the ensemble distribution of the 18 values of standard deviation (one for each member) were compared between the historical and future periods using a two-sided -test to assess the statistical significance. For the Barents Sea, a few combinations had ensemble members with no tracks in May, and are masked out in black whole, indicating that the decrease in early winter is partly compensated by an increase in spring.
In our previous study, Landgren et al. (2019), we applied an MCAO index similar to Kolstad and Bracegirdle (2008) as a coarse-scale proxy for polar lows in the CESM-LE and found a wintertime decrease of strong MCAO in the Nordic Seas. The reductions were strongest at the beginning of the winter season, which we can now confirm by downscaling and tracking. While we found that there were more months with statistically significant changes for smaller size criteria, we also note that Day et al. (2018) studied Arctic synoptic-scale cyclones in CESM-LE and saw a reduction in winter cyclogenesis and intensity. Terpstra et al. (2016) and Michel et al. (2018) cautioned that the SST-500 criterion may favour selection of so-called reverse-shear polar lows (over forward-shear), but they also noted that the most intense polar lows were of this type. Forward-shear cyclones are cyclones developing on regular baroclinic flows where the thermal wind has the same direction as the mean flow. Such systems are difficult to separate from regular but vertically shallow cyclones, and are not the main interest of the present paper. Michel et al. (2018) detected ∼ 61% of the cases from Noer et al. (2011), of which only about 12% had a genesis SST-500 below 40 K. It is nevertheless a possible limitation in our methodology that we may exclude some weaker forward-shear polar lows. This could be more of importance in the Barents Sea, where a larger fraction of the polar lows are forward-shear than in the Nordic Seas (Terpstra et al., 2016). As with any statistical analysis based on thresholds obtained in historical climate, it is not certain that the SST-500 threshold used for polar low detection may hold under future climate change. For example, it could be that a future elevated tropopause invalidates this relation. It could also be that a reduction in intense reverse-shear systems is balanced by an increase in weaker forward-shear systems (Michel et al., 2018).
As high SIC prevents surface heat fluxes from the ocean to the atmosphere, we initially speculated that ensemble members with more ice would form fewer polar lows. However comparing the polar low frequencies of the four members with highest bias and the four with lowest bias does not show any statistically significant differences (not shown). This suggests that for polar low frequency the influence of variability is larger than that of mean SIC bias. Still, additional sea ice in the Barents Sea shifts the genesis region towards the west and may cross the border between the two areas, which could also partly explain the positive frequency bias in the Nordic Seas. We performed a brief analysis of individual polar lows that go from one domain to the other ( Figure S5) and found that the future periods have relatively fewer Barents Sea polar lows originating from the Nordic Seas; instead, more form locally, in agreement with more open water in the Barents Sea. Rae et al. (2017) used 30 years of reanalysis data as well as 100 years of climate model simulations and cautioned against drawing conclusions on potential links between Arctic sea ice and Arctic (summer) cyclone statistics from only one tracking algorithm, one model resolution and limited periods (typically 30 years). In our approach, we believe we have partially taken care of this, using over 600 years of simulations in total, as well as showing the sensitivity of the trends to some of the thresholds in the tracking algorithm. However it is worth noting that the years cannot be considered as completely independent samples because they may be influenced by longer-term (e.g. decadal) variability. However, with a relatively large number of members and the fact that the ensemble members were perturbed in the year 1920, we consider it unlikely that enough of them would be synchronous in decadal or multi-decadal oscillations to contribute a significant difference. While the fact that the model members are from the same ensembles is to our advantage when separating variability and external forcing, it is also a disadvantage that they may share the same deficiencies related to model processes, e.g. perhaps as manifested in the sea-ice bias discussed earlier.
While we have only used one scenario (RCP8.5), results for the near future (2026-2035) are expected to be similar regardless of choice of scenario. The frequency decrease in the later period (2071-2080) could be less pronounced in RCP4.5 or RCP2.6. Zahn and von Storch (2010) used three future scenarios from CMIP3 (A1B, A2 and B1) and found that the differences between them were small compared to change from the historical period, although we should note that in particular those three scenarios have much smaller differences than the RCP scenarios.
Because these results are based on one ESM and one RCM, each of which may have its own biases due to model process uncertainties, further studies could include other model combinations to improve the estimate of range of likely outcomes. It is also clear that many GCMs/ESMs still have large sea-ice biases. While it may be possible to use bias correction of SIC (e.g. together with an atmosphere-only GCM as a step between ESM and RCM), any adjustments break the physical consistency provided by the ESM between variables, and is very challenging in the marginal ice zone. It may also be worthwhile to consider producing downscaled results from a limited multi-model subset of CMIP5 with smaller biases to see if any trends in polar lows in the Barents Sea can be detected. Nevertheless, studies using a large number of perturbed simulations from a single model complement multi-model ensembles by allowing separation of internal variability and external forcing.

CONCLUSIONS
Eighteen members of CESM-LE downscaled with HCLIM-ALARO over the Nordic and Barents Seas showed a future decline in Nordic Seas polar lows in the RCP 8.5 scenario. Statistically significant changes were found both in frequency, lifetimes and intensity. Analysed over the October-May season, the strongest decreases were found at the beginning of the season (November-January), for which they were statistically significant for most combinations of track selection criteria. Similar to what was found for CMIP3 models in Zahn and von Storch (2010) and Kolstad and Bracegirdle (2008) and CMIP5 models in Romero and Emanuel (2017), the future retreat of the sea ice in the CESM-LE leads to a northward shift in the genesis region. However, with the large number of members from the same model, in this study we were additionally able to assess the statistical significance of the results in relation to internal variability. Although the variability itself was largely unchanged, it could now be better separated from the external forcing. By considering variability as noise on top of the externally forced signal, we were able to separately identify months with statistically significant contribution of external forcing.
The decline in frequency in the beginning of the winter season and increase in March leads to a delayed polar low season in the future. This is in agreement with the atmospheric static stability increase being distributed differently throughout the year, as discussed by Landgren et al. (2019). For lifetimes of polar lows, there was a statistically significant reduction in November-January in the Nordic Seas.
The results for the Barents Sea were less conclusive. Although the simulations show a statistically significant increase in polar low frequency in February-April, the results are difficult to interpret in relation to the real world because the historical simulations have too few polar lows (relative to the downscaled reanalysis) due to excess sea ice in the CESM-LE in this region. Future sea-ice reduction opens up new areas for polar low development, balancing the reduction expected from the atmospheric static stability increase, making the frequency change uncertain in the Barents Sea.