Kilometer-scale modeling projects a tripling of Alaskan convective storms in future climate

Convective storms produce heavier downpours and become more intense with climate change. Such changes could be even amplified in high-latitudes since the Arctic is warming faster than any other region in the world and subsequently moistening. However, little attention has been paid to the impact of global warming on intense thunderstorms in high latitude continental regions, where they can produce flash flooding or ignite wildfires. We use a model with kilometer-scale grid spacing to simulate Alaska’s climate under present and end of the century high emission scenario conditions. The current climate simulation is able to capture the frequency and intensity of hourly precipitation compared to rain gauge data. We apply a precipitation tracking algorithm to identify intense, organized convective systems, which are projected to triple in frequency and extend to the northernmost regions of Alaska under future climate conditions. Peak rainfall rates in the core of the storms will intensify by 37% in line with atmospheric moisture increases. These results could have severe impacts on Alaska’s economy and ecology since floods are already the costliest natural disaster in central Alaska and an increasing number of thunderstorms could result in more wildfires ignitions.


Introduction
Alaska is strongly impacted by climate change, mainly because of the polar amplification of warming due to melting sea ice (Ohring and Adler 1978;Markon et al. 2012). Alaska's temperature is projected to increase by 2-3 • C under RCP2.6 to 6-9 • C under RCP8.5 until the end of the century due to anthropogenic emissions of greenhouse gases according to general circulation models Stocker et al. (2013). Recent research found that mean and maximum temperatures in Alaska could increase as much as 8 • C by 2100, and extreme minimum temperatures could increase by up to 22 • C under a high-end emission scenario (Lader et al. 2017;Newman et al. 2020). This has dramatic effects such as thawing permafrost, increased risk of wildfires (Partain et al. 2016), sea ice loss, increased likelihood of rain-onsnow events (Musselman et al. 2018), and more intense precipitation (Bennett and Walsh 2015). Such changes would have major consequences on the local economy, ecosystems, and the life of local populations.
The effects of global warming on extreme precipitation in Alaska have already been well studied. It was shown that a higher fraction of precipitation will come in the form of very heavy precipitation events (Bennett and Walsh 2015), and that the sensitivity of these events to global average temperature is increasing with latitude (Westra et al. 2013), therefore leading to a very strong response in the Arctic region. However, existing modeling studies either use models that have a too coarse resolution to reliably capture extreme precipitating storms (Lader et al. 2017;Bieniek et al. 2014;Bennett and Walsh 2015) or focus on climatological aspects of precipitation changes (Newman et al. 2020). Little is known about the impacts of global warming Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s0038 2-020-05466 -1) contains supplementary material, which is available to authorized users.
1 3 on high latitude intense convective systems, despite the fact that mesoscale convective systems have been observed in high latitude regions such as Finland (Punkka and Bister 2015) or Siberia (Жукова et al. 2018) where they produced high impact weather. These systems are very rare under a current climate, and correspond to extreme events in the Alaskan climate. This differs from lower latitude climates where mesoscale convective systems are common in the warm season.
Considerably more studies have focused on the response of intense convective systems to global warming in mid-latitude regions and showed that climate change can modify the frequency of occurrence of environmental conditions conducive to the development of severe weather (Diffenbaugh et al. 2013;Gensini and Mote 2015;Rasmussen et al. 2017). It was also shown that the rainfall volume from mesoscale convective systems might strongly increase due to temperature increases in the USA (Prein et al. 2017b). Finally, the increases in storm runoff (the excess of hydrological runoff induced by the presence of a storm) across the globe can be much larger than average precipitation increases (Yin et al. 2018), and convective precipitation has a stronger response to climate change than other types of precipitation and therefore increasingly dominates the precipitation extremes (Berg et al. 2013). These results, together with Arctic amplification and the strong warming and moistening of Arctic regions, are the basis for our hypothesis that Alaska's convective systems will experience a strong response to global warming.
In this study, we focus on the impact of global warming on deep convective summer storms in central Alaska. These storms typically occur during the period of highest insolation (from mid-June to July) and produce heavy downpours in central Alaska (Grice and Comiskey 1976;Reap 1991). Moreover, floods were found to be potentially the most costly effect of global warming in Alaska, especially in the central part of the state (Melvin et al. 2017).
To study these storms, we track convective precipitation areas (Prein et al. 2017a) in a present and future climate convection-permitting regional climate simulation that includes Alaska and its surrounding area (Monaghan et al. 2018). The benefits of using convection-permitting modeling for representing convective systems and mesoscale convective systems have been demonstrated in mid-latitude and tropical regions (Prein et al. 2017a;Crook et al. 2019;Chang et al. 2018), and these models have been shown to better reproduce the lifetime, life cycle and diurnal cycle, spatial frequency, speed patterns, and spatial extent distribution of convective systems compared to coarser-resolution models. However, convection-permitting models have the tendency of overestimating peak precipitation rates (e.g., Crook et al. (2019); Prein et al. (2020)) and are computationally very expensive, which typically only allows simulations over decadal time periods.
There are two questions that we answer in this paper. (1) How well can a convection-permitting climate model reproduce summertime hourly precipitation statistics? (2) How do convective systems change in Alaska due to global warming?
After this introduction, Sect. 2 will present the datasets that were used for this study, and provide an assessment of the model with station data to answer the first question. In Sect. 3 we show a convective system climatology of Alaska based on storm tracking. Section 4 will address the impact of global warming on convective systems, and Sect. 5 will provide a summary and a discussion of the main results, as well as insights for the next steps after this study.

Region of interest
In this study, we focus mainly on continental Alaska (Fig. 1), where convective systems typically occur in the summer season due to solar heating. Thunderstorms are typically observed between May and August (Grice and Comiskey 1976). We define continental Alaska as the five climate regions shown in Fig. 1, which were defined by Bieniek et al. (2012). The North Slope, North of the Brooks range, is dominated by an Arctic climate. Therefore, summer storms are very rare in this region, except over the Brooks range where mountains can trigger storms. The West Coast is exposed to a maritime influence in ice-free conditions (in summer and early autumn) and is also subject to the influence of sub-polar cyclones and rare polar lows hitting the coast and bringing precipitation to the region (Perica et al. 2012). Interior Alaska has a continental climate and reaches the warmest temperatures of the state in summer. This is where the summer convective activity is strongest (Grice and Comiskey 1976;Reap 1991). In all five regions, the storms are typically triggered by solar heating under weak vertical wind shear (Grice and Comiskey 1976) and can be intensified by upper-level troughs, cold fronts, or remnants of tropical systems crossing the state (Perica et al. 2012).
The other climate regions of Alaska that are adjacent to the Pacific ocean are subject to a more maritime climate and are dominated by orographic precipitation and Pacific mid-latitude extratropical cyclones hitting the coastal range (Bieniek et al. 2012).

Model data
This study uses regional climate model data over Alaska covering the period 2003-2015 (Monaghan et al. 2018;Newman et al. 2020). The Advanced WRF model version 3.7.1 (Skamarock et al. 2005) is used to downscale the European Centre for Medium Range Weather Forecasting Interim Reanalysis (ECMWF ERA-Interim) (Dee et al. 2011). The domain covers the entire state of Alaska, as well as a large part of Northwestern Canada including the Yukon territory (Fig. 1). The model horizontal grid spacing is 4 km, enabling convection to be explicitly simulated by the dynamical core of the model and the microphysical scheme (see  for a review). The convective parameterization is turned off, and Thompson microphysics (Thompson et al. 2008) are applied. Of course, there will always be unresolved scales in atmospheric modeling and 4 km grid spacing can only capture the largest modes of convection. The significant improvements when using kilometer-scale models over larger-scale models are well documented in climate  and weather forecasting research (Clark et al. 2016). While there are well known deficiencies of kilometerscale models (e.g., Lebo and Morrison (2015)), the benefits of not using a deep convection scheme at these scales dominate (e.g., Vergara-Temprado et al. (2020); Mahoney (2016); Prein et al. (2020)). Accumulated precipitation is stored at an hourly frequency. A previous study that evaluated this simulation showed that the model can capture the general climate conditions of Alaska Monaghan et al. (2018).
For the future climate, the pseudo-global warming (PGW) method that was first proposed by Schär et al. (1996) is applied. Monthly mean climate change perturbations of temperature, moisture, wind, and geopotential height are added to the 6-hourly boundary conditions provided by  (Amante and Eakins 2009). The five climate divisions that are used in this study (Bieniek et al. 2012) are shown in thick black contours with names in bold. The eight precipitation stations that provide hourly precipitation observations for model evaluation are also shown, as well as dominant geographic features. The top-left inlay shows the simulation domain, that is much larger than the region of interest ERA-Interim. These perturbations are obtained from an ensemble of 19 global climate models from the Coupled Model Intercomparison Project (CMIP5) (Taylor et al. 2012), and are derived from the difference between monthly climate conditions of 2071-2100 compared to 1976-2005 under the RCP8.5 scenario (Riahi et al. 2011). This scenario is a high-end emission scenario; it stands in the upper end of the uncertainty range of emissions under current policies (Capellán-Pérez et al. 2016;Rogelj et al. 2016). A more refined technique has been applied to account for sea surface temperature and sea ice concentration changes, as described in Newman et al. (2020). Sea ice concentrations are derived from the CMIP5 ensemble median. Therefore, the method that is used in this study accounts only for the thermodynamic impacts of climate change, as well as impacts of dynamics at large spatial ( > 2000 km) and temporal ( > 1 month) scales. Circulation changes that happen over time periods shorter than a month are not accounted for, in particular changes in transient eddies that account for most of the moisture transport from the lower latitudes (Naakka et al. 2019). More details about the PGW method can be found in Newman et al. (2020), and more details of the model setup can be found in Monaghan et al. (2018).

Observations
Alaska is a data poor region when it comes to meteorological observations with a very low-density station network and poor satellite coverage. We combined a large range of observational datasets including rain gauges, lightning data, aircraft observations, and radar data to evaluate the performance of the current climate simulation.
Precipitation data The Hourly Precipitation Network (Center and NESDIS, NOAA, 2020) provides sixteen stations for model evaluation. Eight of these stations are within the region of interest (shown in Fig. 1) and provide data for the period 2003-2013. Those stations use tipping buckets to collect precipitation that have a volume of either 0.25 mm or 2.5 mm depending on the stations. Tipping buckets accumulate precipitation until they are full. When emptying, the precipitation amount of the size of the bucket is recorded.
We emulate this behavior with the model precipitation to reduce the effect of the collection method on model evaluation. Precipitation from the nearest grid cell is compared with each station. We also tested distance-weighted precipitation averages of the four nearest grid points to each station, which yielded similar results for all the stations (not shown).
Lightning data The Alaska Lightning Detection Network, managed by the Alaska Fire Service, provides lightning data with the location and time of lightning flashes detected by sensors during the period 2003-2016 (Alaska Interagency Coordination Center -Alaska Fire Service 2020). This data is used as a proxy to evaluate the convective storm climatology, and each flash is counted as one lightning (we do not account for multiplicity). It should be noted that the sensors have been updated in 2007. After this year, they also recorded cloud-tocloud lightning that was not taken into account before. This led to a strong increase in the number of detected lightning and avoids any interpretation of a temporal trend. However, the diurnal and annual cycles, as well as the spatial patterns of lightning strikes, should not be systematically affected.
Radar Data Data from the Doppler detection Radar (RAdio Detection And Ranging) located in Fairbanks was also used in this study. This is part of the NOAA Next Generation Radar (NEXRAD) Level 3 Products (NOAA National Weather Service (NWS) Radar Operations Center 1991). The radar scan time is of approximately 10 min and the detection range is 124 nautical miles (230 km). The dataset provides an estimation of 1-h accumulated precipitation, based on reflectivity to rainfall rate (Z-R) relationships. Since only two weather radars are located in the region of interest (in Fairbanks and Nome) this dataset could not be used for an evaluation of the climatology of storms in Alaska. However, the Fairbanks radar is used for qualitative comparison of the storms between the model and the reality.
Radiosonde data Data from the Fairbanks meteorological station have been retrieved from the Integrated Global Radiosonde Archive (IGRA) available from the National Centers for Environmental Information (NCEI) (Durre et al. 2006). Each day during the simulation, a weather balloon sounding was done at 00:00UTC (16:00 local time) and surface-based Convective Available Potential Energy (CAPE) and Convective Inhibition (CIN) were calculated for this sounding as follows: where T p is the temperature of an air parcel lifted adiabatically from the surface, T e the temperature of its environment, and LFC and EL refer to the altitudes of the level of free convection and the equilibrium level (defined as the two first levels where T e = T p ). This enables comparison with the modeled surface-based CAPE and CIN at the nearest gridpoint in the historical simulation. 22% of the soundings are missing in this dataset. Therefore, we remove dates of missing soundings from the model data when comparison is made to avoid any statistical bias in the comparison.

Uncertainty quantification
The pseudo-global warming technique adds only multidecadal, multimodel means to the boundary conditions. Therefore, the control and the future simulations both feature similar weather that experiences the same modes of natural climate variability. However, the model is not spectrally nudged, and present and future climate simulations can differ on synoptic scales within the model domain. This introduces differences in the climate change signal beyond pure thermodynamic effects from the pseudo global warming approach. To estimate the range of internal variability and interannual variability, this study uses the two-tailed bootstrap uncertainty test by blocks, presented in Hesterberg (2015). Monthly blocks were used, assuming that summertime precipitation statistics are independent from one month to another, with a sample size of 1000 for each test. For each bootstrap sample, the same shuffling has been used between the historical and PGW simulations which enables to also apply the bootstrapping to the changes between the two simulations. This is made possible by the fact that monthto-month variability is preserved in the PGW method.
The results of this test must be taken with caution since it is only conducted over a 13-year period, whereas Alaska climate is largely impacted by decadal and multidecadal oscillations such as the Pacific Decadal Oscillation (PDO) and the El Niño -Southern Oscillation (ENSO) (Bieniek et al. 2014). The simulation period consists of predominantly negative PDO phases (and negative ENSO phases as well to a lesser extent). However, Bieniek et al. (2014) demonstrated that precipitation is less affected than temperature by these oscillations, especially in summer.

Analysis of the precipitation rate distribution
Much of the previous literature about climate change effects on precipitation intensity often focused on the theoretical scalings of precipitation with temperature/moisture increases. However, this approach does not capture all the subtleties of the changes in precipitation characteristics that can be much more complex. Here we are interested in changes in the full distribution of the precipitation rates.
To manage this highly skewed distribution, the "Analyzing Scales of Precipitation" (ASoP) technique (Klingaman et al. 2017) is used. First, we calculate a histogram of the precipitation rate with exponential bins: the bin width increases exponentially with increasing precipitation rates so that bins have an equal width on a logarithmic scale: b n = b 0 e n . When represented on a logarithmic scale, the usual graphical properties of a histogram are conserved as follows. Let p(r) be the subsequent distribution of the rain rate on a logarithmic scale, so that p(r) is the occurrence frequency of the precipitation rates between r and re . The contribution of the precipitation rates between r and re to the mean precipitation rate P is then : rp(r) . This curve represents the actual contributions of each precipitation rate to the total precipitation. P is the area under this curve. That is, if N bins are used from r 0 to r 0 e N : P = ∑ N n=0 r 0 e n p(r 0 e n ) . This curve can be divided by the mean precipitation rate P, to yield the fractional contributions of each precipitation rate to the mean precipitation. Please see Berthou et al. (2018) for a detailed description of the technique (see their Fig. 2 for a summary).

Tracking of hourly precipitation objects
In this study, a tracking algorithm is used that identifies storms according to the simulated accumulated hourly precipitation. The algorithm is a Python implementation of the method described in Davis et al. (2006); Prein et al. (2017a) and was originally designed to track mesoscale convective systems over the central and Eastern United States. It consists of three main steps : -Smoothing the precipitation field in space and time with a gaussian filter. The smoothing radius is fixed to be 1 h for the time dimension, and an integer of r grid points for the space dimensions (x and y). -Thresholding the smoothed precipitation field by only considering precipitation above a rate p min (in mm/h). This results in a binary file with all grid-cells above the threshold being one and all other cells being zero. Based on the binary file, connected precipitation objects are identified. Within an object cells have to be connected in the longitude, latitude, or time dimension (also diagonal adjacent cells are regarded as connected). This step provides three-dimensional objects. Each contiguous object in the 3D x, y, t space is considered to be a storm. -We only consider the storms that have a volume (number of connected grid cells in the latitude, longitude, and time dimension) exceeding a threshold n min .
This algorithm is applied for all summer months (May-August) during the thirteen simulated summers (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), and all the storms that were formed over the region of interest ( Fig. 1) were selected. The algorithm was applied on a monthly basis. This implies that some storms are cut when they occur between 2 months. However, this effect should be small as the number of hours in a month (about 720) is much larger than the lifetime of storms (about 6 h in average). A total of 527 storms was found in the control simulation.

Diurnal and annual cycles
The model ability to reproduce the general climatology of Alaska, i.e. the seasonal mean and interannual variability surface of temperature and precipitation, the duration of the snow season and the distribution of daily temperature and precipitation, has already been demonstrated in Monaghan et al. (2018). In this section, we, therefore, mainly aim at evaluating the ability of the model to reproduce the summertime precipitation climatology of the region of interest at subdaily time scales.  Figure 2a, c show the simulated and observed annual cycle of precipitation. Precipitation peaks in the summer months in the region of interest, and this highlights the importance of this season for the hydrology of the region. The peak occurs earlier in central Alaska and in the North Slope, in June or July, whereas precipitation peaks in August on the West coast. The model and the stations show good agreement, both in the peak precipitation month (Fig. 2a) and the shape of the annual cycle (Fig. 2c), although the model is wetter than the observations in the region of interest, as already noted in Monaghan et al. (2018) (see their Fig. 4c).
For the summertime (MJJA) diurnal cycle (Fig. 2b,d), the model reproduces the timing of the observed evening peak well. Figure 2d shows that the model has a delay in the peak precipitation of about 1-3 h (depending on the station) and is consideraby wetter. The model also has an exaggerated diurnal cycle amplitude, a problem that has been noticed for other convection-permitting simulations in the Eastern United States (Chang et al. 2018). Agreement between the model and observations is much poorer on the Pacific coast, however this region exhibits a very weak diurnal cycle and is not considered here.

Hourly precipitation rates distribution
The distribution of hourly precipitation rates for four representative stations and a comparison with the model results is shown on Fig. 3. The lowest precipitation rates contribute the most to total precipitation. The smallest precipitation rate possible, i.e. the size of the station bucket, was not included as it has a very high contribution (much larger than all the other bins) to total precipitation, because that is actually the integrated value of the contribution of all the precipitation rates below the bucket size to total precipitation. Therefore, this figure only focuses on a part of the precipitation distribution.
The model reproduces the shape of this distribution well. At low precipitation rates, the model is much wetter than the observations, both for annual precipitation and summeronly precipitation. This reflects two biases: first, the frequent undercatch of gauges, especially in Alaska and Arctic regions where the frequent snowy and windy conditions can lead to very strong undercatch (Kane and Stuefer 2015) in addition to re-evaporation of precipitation from the bucket. Second, the model is likely overestimating drizzle, which is a well-known problem in atmospheric models. This problem is solved or even oversolved in some convection-permitting simulations (Berthou et al. 2018) whereas it persists in other kilometer-scale models Dai et al. (2017).
However, the agreement between the model and the observations for higher precipitation rates (over 1mm/h) is good, especially during the summer season. In particular, the tails of the distributions are well reproduced by the model. The better agreement in summer might be partly due to the smaller undercatch from precipitation stations rather than better performance of the model. Therefore, the results presented in Figs. 2 and 3 show that the model is able to reproduce the main features of summer precipitation in continental Alaska at hourly timescales reasonably well. The model is overestimating precipitation amounts, as noticed in Monaghan et al. (2018), however Fig. 3 shows that this is due to an overestimation of drizzle in the model (precipitation below 1mm/h occurs more often in the model than in observations). These results also suggest an undercatch of snow in the stations record since the difference between model and observations are smaller in the warm season compared to the whole year, especially in the station of Nome. The effect of these differences on this study should be minor, since we focus on intense precipitation from summer convective systems.

CAPE and CIN distributions
Supplementary figure S1 shows the distribution of surfacebased CAPE and CIN in the model and as measured by weather balloons in Fairbanks. Modeled and observed CAPE distributions are very close, and the uncertainty ranges overlap considerably. This shows that the model is able to represent observed atmospheric instability with very good confidence, and in particular represents well the frequency of occurrence of moderate ( ≈ 100 J/kg) and large ( ≈ 400 J/ kg) instability. This is crucial for representing well the initiation and structure of convective systems. The model is also able to reproduce the distribution of CIN well. There is less overlap of the uncertainty intervals of the two distributions, which suggests that strong values of CIN are probably too infrequent in the model compared to observations. However, this should not have a strong effect on convective initiation, as very small values of CIN that do not represent a strong barrier to convection are typical in this region as shown by the distribution. (By comparison, CIN is typically larger by an order of magnitude in the US Great Plains (Rasmussen et al. 2017)).

Evaluation of the tracked storm systems
After evaluating the models ability in simulating fundamental characteristics of the rainfall climatology in Alaska we will focus our analysis on the evaluation and future changes of organized convective storms in the rest of the paper. Figure 4 depicts the evolution of an organized thunderstorm that was observed in the vicinity of Fairbanks on July 16th, 2005. Isolated, intense convective cells developed East of Fairbanks around noon, and then grew in size and merged while moving northwestward during the afternoon, eventually producing a convective line of 200 km in length at 4:30 PM. This storm then vanished within a few hours. Another example is provided in Figure S2 that shows a system that formed in the evening of July 4th, 2015 West of Fairbanks and took a roughly circular shape with a diameter of around 70 km while moving eastward. The system experienced two active phases at 9 PM and 1 AM making it long-lived despite its moderate precipitation rates. These examples show the morphology of convective systems in central Alaska. Those systems can live for several hours with spatial extensions of dozens or even hundreds of kilometers. Although a direct comparison between the observed and modeled storms is not possible due to the chaotic nature of convection we can compare the simulated and observed storm characteristics qualitatively.

Tracking organized convective systems
An example of the result of the storm tracking is shown on Precipitation produced by the model (Fig. 5a) reproduces this situation well, with patchy and intense precipitation, that corresponds likely to convection or deep convection. The algorithm identified four storms on this day: two storms south of the Brooks Range, one storm south of the West coast, and one in the plains of Central Interior Alaska. The location of these storms corresponds well to the location of the deep convective clouds in the satellite data. This indicates that the simulation has skill in capturing observed deep convective systems despite the fact that no spectral nudging was used. Only the largest and the most intense convective cells or clusters are selected by the tracking algorithm, as the model also produces other cells of limited extent on the southern slopes of the Brooks range and in interior Alaska. We decided to focus our analysis on the largest and most intense storms due to their potential impacts on society and the ecosystems. Figure 5c shows the position of the storms at different times. The storms move very slowly. The tracker identifies multi-cell, organized convective systems. They have a large spatial extent (a few hundred of kilometers) and live for multiple hours and therefore can be classified as mesoscale convective systems (Houze 2004). Figure 5d-g shows the evolution of one storm: it consists of three different cells that merge and intensify, reaching its maximum intensity around 18:00-20:00 and then slowly decays. Qualitatively, this storm is similar to the ones observed in radar observations in the surroundings of Fairbanks, in terms of storm size, shape, precipitation intensity and lifetime (see Figs. 4 and S2). The life cycle is also similar as both observed and modeled storms are formed by the merging of individual convective cells. This indicates that the model is able to simulate the morphology and dynamics of organized convective systems in Alaska.

Dependence of results on the tracking parameters
The tracking algorithm can be modified by three parameters: a precipitation threshold p min (in mm/h), a smoothing radius of r grid points, and a minimum storm volume of n min grid points. In the original version of the algorithm, these parameters were set respectively to 5 mm/h, 8 grid points, and 2000 grid points. These parameters need to be adapted, since Alaska's summer convective systems are typically smaller and less intense than their mid-latitude and tropical counterparts. Four different configurations are tested that produce a reasonable number of tracked storms (i.e., several hundreds over the simulation period). Then, the tracking algorithm is run with these different parameters for the month of July 2015 (a month with strong convective activity across the state) in order to be compared. The parameters for these simulations are shown in the top left insert of Fig. 6.
The different configurations exhibit similar tracks and identify storms at the same locations (not shown). In particular, all configurations of the parameters lead to a spurious identification of storms on the Pacific coastal range and the Alaska range, that are actually areas of long-lived, intense orographic precipitation. Removal of these areas susceptible to intense orographic precipitation led us to our region of interest defined in Fig. 1. In particular, the Pacific Coast and the windward slopes of the adjacent mountains have been removed, whereas mountain ranges protected from the Pacific Ocean have been conserved in the Southeast Interior region. Within the region of interest, 72% of the rainfall produced by the storms that are identified by the tracker are classified as convective according to the separation   Figure 6 shows the distribution of the lifetime, hourly rainfall volume, and maximum hourly precipitation rate for the different sets of parameters tested for tracking precipitation objects. The different configurations produce very similar distributions of maximum precipitation rates. The distributions of lifetimes are also similar, except for configuration A that identifies many long-living storms of 10 hours or more. This seems unrealistic and is probably due to a misclassification of orographic rainfall as convective systems, and might be due to the lower precipitation threshold used in this configuration. Finally, the distribution of precipitation areas varies much more between configurations. This appears to be an effect of the smoothing radius r and the precipitation threshold p min : the larger the radius and the smaller the threshold, the larger the identified storms.
Configuration B was selected for this study, mainly because it selects storms with slightly higher precipitation rates, and avoids misclassifying many orographic storms as convective systems. Finally, it has a smoothing radius large enough to avoid splitting one connected storm system into several individual storms.
Moreover, Fig. S3 shows that the storms identified by the tracking algorithm with the different configurations have very similar tracks and locations over the region of interest. Some configurations produce longer tracks or split one track into several tracks, however, it seems that the physical systems identified as storms are similar using the described settings, and do not depend a lot on the algorithm parameters. Therefore, the tested parameters of the tracking algorithm do not have a substantial influence on the number and characteristics of identified storms.

Climatology of central Alaska summer storms
The annual and diurnal cycle of storms are shown on Fig. 7 in thick solid lines, as well as the cycles for lightning data in dotted lines. We see a strong annual cycle with few storms in winter and spring, and a sharp peak in June and July corresponding to the season of maximum insolation. This result is compatible with lightning observations, airplane observations conducted by Grice and Comiskey (1976), and the lightning data studied in Reap (1991), which show a very short season for convective systems with storms mainly in June and July, and a few storms in May and August. A qualitative comparison with lightning data and the previous studies show that the model seems to slightly overestimate the number of storms from August to October. This might be due to intense extratropical cyclones that are advected from the sea and are identified as storms. Indeed, these late summer and autumn storms are mainly found in the West Coast and Southeast interior regions (Fig. S4) whereas the Central Interior, Northeast Interior, and North Slope regions exhibit a much narrower stormy season. Additionally, the observations and model periods are fairly short and do not overlap, which introduces differences due to climate internal variability and climate change. This sharp peak of convective storm activity in the mid-summer is also found in Siberia (Жукова et al. 2018) whereas the season for MCSs is extended from April to November in Finland (Punkka and Bister 2015). We also find a pronounced diurnal cycle of storms frequencies that is peaking in the evening at local 9PM. The cycle is tarnished by a strong relative uncertainty, making it statistically non-significant. However this uncertainty is largely reduced in the PGW simulation that has a more pronounced diurnal cycle with a similar shape and peak timing (see Fig. S4). The diurnal cycle of storm systems is compatible with the rain gauge observations and the diurnal cycle of precipitation in general as shown in Fig. 2. However, lightning data shows a diurnal cycle that peaks earlier during the late afternoon, which is more in agreement with previous studies that have focused on the diurnal cycle of thunderstorms in Alaska : Grice and Comiskey (1976) noticed a peak of the deep convective activity in the early afternoon, and Reap (1991) in the early evening. The strong delay between the model and the observations is unclear and might be due to several factors. First, the model has a tendency to produce a diurnal cycle peaking slightly later than the observations in the West Coast and Central Interior (see Fig. 2). Second, here we focus on the largest and most intense storms, mainly organized convective systems, that generally peak later than the more frequent, smaller deep convective cells that were considered as well in the studies mentioned above and might dominate the signal of lightning data (we recall that the events we focus on in this study are extreme events and do not dominate warm-season hydrology of Alaska). Finally, it has been shown that lightning is strongly correlated with the convective area of the storms and therefore peaks at the beginning of the storms. Mattos and Machado (2011) have found that 80% of lightning in mesoscale convective systems in Southern Brazil occurs before the storms mature, which might explain why lightning in our analysis occurs when the model produces storm genesis. As visible on Fig. 7b, storm genesis has a much sharper, statistically significant diurnal cycle than the number of storms. This cycle peaks earlier at local 7 PM and matches better the lightning observations, although a delay of about 2 h persists, which is consistent with the observed delay in the diurnal cycle of precipitation (Fig. 2). After genesis, most of the precipitation in the system is likely stratiform. It can still be heavy, but does not produce any lightning (Houze 2004).
We also underline that this high latitude region experiences a weaker diurnal cycle that mid-latitude regions. For example, in Fairbanks in early July (period of maximal convective activity), the sunset happens around 00:30 AM local time and at the zenith (around 2PM), the sun is about 45 • above the horizon. Therefore solar forcing could help sustain or even trigger storms until late into the evening. Figure 8 shows the distributions of lifetime, maximum precipitation rate, area, and speed of the tracked storms in the control simulation. The speed is determined using the hourly change in the location of the center of mass of the storms. The storm lifetime is on average 5-h but some storms exist for more than 10 h. The maximum precipitation rates are very high compared to usual precipitation in Alaska, reaching values of up to 40 mm/h. Finally, the distribution of storm area is highly skewed, with most storms having areas of 1000 km 2 to 4000 km 2 (corresponding to a diameter of 35-70 km assuming storms are circular). Some storms can reach a very large extent, up to 12000 km 2 (i.e. a diameter of 120 km). This shows that we are tracking large, intense, and long-lived storms that are typically organized convective systems, or even mesoscale convective systems for the largest ones. Finally, these systems generally move very slowly at an average of 16 km/h and very few storms have a speed exceeding 30 km/h. The slow speed, in addition to their large size and high rainfall rates, can potentially lead to very large amounts of precipitation in mesoscale catchments. In comparison, convective systems in the central United States move much faster (approximately 35 km h −1 ) Prein et al. (2017a), have higher precipitation rates (up to 100 mm h −1 ), and can reach an extension of several hundreds of kilometers. Figure 9a shows the storm track density, calculated from the control simulation. Figure 9b shows thunderstorm density obtained from airplane observations during the 1970s (Grice and Comiskey 1976) for a qualitative comparison of spatial patterns. Note that there are likely major differences between the observed and tracked thunderstorms (e.g., in terms of their size, likelihood of detection, underlying climate conditions) which only allows a qualitative comparison of the spatial patterns. Some spatial patterns agree well, for example, a thunderstorm maximum over the Yukon-Tanana uplands (this maximum is much stronger in the observations than in the model) (refer to Fig. 1 for the names of the natural features of Alaska). Also, both the model coupled with There are however some discrepancies as well : a local maximum is found over the Brooks range in the simulation, that is not present in the observations (however this is likely impacted by the Brooks range being out of the usual visual limit of the flight track that was used in Grice and Comiskey (1976)). Finally, it is obvious that the algorithm does identify spurious storms on the Alaska range and Wrangell mountains. This is likely orographic precipitation over the highest tops. This kind of problem has already been reported in Poujol et al. (2019), where intense orographic precipitation was often misclassified as convective precipitation over the highest tops of Southern Norway. Figure S5 compares the storm track density with the lightning density over the region of interest. Although lightning density exhibits much less small-scale variability than storm density from this study and Grice and Comiskey (1976), lightning is also prevalent in Interior Alaska South of the Yukon River, which comforts our results.
Finally, Figure S6 shows that track density is correlated with topography both in the historical and PGW simulations. The Pearson's correlation coefficient is about 0.35 and is statistically significant (p value<0.001). This is consistent with previous studies (Grice and Comiskey 1976;Reap 1991) that observed that more convective systems occur at higher altitudes and are likely triggered by orography. However, this dependence appears to be very strong in the model compared to Grice and Comiskey (1976) and storms are probably too tied to orography, as it is also suggested by comparison with lightning data in Fig. S5. This issue has been reported in other convection-permitting models (e.g., Purr et al. (2019)) 4 Impact of global warming on Alaska's storms Figure 9b shows the storm density in the pseudo-global warming (PGW) simulation for the summer months (MJJA). The general pattern of storm frequency is preserved, with a well-pronounced maximum in Yukon-Tanana uplands. Most importantly, the average frequency of storms is three times higher in the future climate simulation (527 in the control simulation and 1558 in the PGW simulation), leading to significant increases in the storm track density. Figure 9d shows the relative change in the tracked storm frequency with an increase of a few dozens to a few hundred percent South of the Yukon River, in the regions regularly hit by storms in the control simulation. Some regions where no storms occurred in the current climate simulation experience a regime shift with frequent storms in the PGW simulation. This is the case for some portions of the West Coast like the coast of the Kotzebue Sound and the Yukon delta, as well as the North Slope region. In particular, the region north of the Brooks range might experience a shift during Arctic summers, from current nearly precipitation free to more continental midlatitude-like summers with intense convective precipitation occasionally occurring by the end of the century. Figure 9e and f shows the annual cycle of storms in the control and PGW simulations. The dramatic increase of the number of storms is mostly concentrated in the summer season. The convective system season is extending in some regions. Firstly, the number of storms increases substantially in May in nearly all the subregions, so that future Mays will have almost the same amount of storms as the current August. Second, one can see an extension of the season towards September and October on the West Coast and, to a lesser extent, in Central Interior Alaska. (The storms identified in autumn and winter in the Southeast interior likely do not correspond to convection but to orographic precipitation as already noted). Individual examination of those autumn storms revealed that they are mostly systems that form over the Bering Sea and are advected onshore. They are likely extratropical cyclones with large regions of intense embedded convection.

Storm frequency
Our results are in agreement with high latitude observational studies. Some studies focused on the influence of global warming on continental high-latitude deep convection and noted a strong increase in the frequency of cumulonimbus clouds in most of Siberia, with significant increases of about 10% per decade in summer and up to 20% per decade in spring and autumn (Chernokulsky et al. 2011;Алексеев et al. 2014). The larger increases in the transition season result in a widening of the deep convective season. This differs from this study for Alaska, where the relative increases are the same during all seasons and the extension of the transition season is therefore less pronounced. Moreover, a large number of studies noted an increase in convection in general over the high latitudes (e.g., Ye et al. (2017), among many others). Fig. 9 a Density of storms in the present climate simulation, defined as the number of hours per year (MJJA only) that a grid-point is covered by a storm. b Same, for the PGW simulation. c Thunderstorm density (in days per year) from airplane observations during a 6-year campaign in the 1970s. The figure is reproduced from Grice and Comiskey (1976). d Track density difference between the PGW and the historical simulations is shown at locations where the average track density is above 0.25 h/year. Hatchings show regions with storm frequencies of less than 0.25 h/year in the historical simulation and with more than 0.25 h/year in the PGW simulation. e Annual cycle of storm frequency in the region of interest. Shadings show the results of a two-tailed 98% bootstrap uncertainty test. f Same as e, but separated by climate region. Uncertainties can be seen on Fig.  S4c, d. In a-d Note the logarithmic scales of the color bars. Brown contours show smoothed topography with contours at 250, 500, 750,1000,1500, 2000 and 3000 m. The outline of climate regions and rivers have been added for easier comparison between the figures. Note that a gaussian smoothing with a standard deviation of 5 was applied to storm track density to remove small scale variability ◂

Precipitation rate distribution
According to the Clausius-Clapeyron law Clausius (1850), the atmospheric moisture content is increasing nearly exponentially with temperature, by approximately 7% per degree Celsius of warming. This means that saturated moisture of air will increase strongly with climate change, which is expected to lead to a similarly large increase in extreme precipitation rates Trenberth et al. (2003). High-latitude regions might be especially vulnerable to this increase of extreme precipitation since they experience much larger warming than the global average, and therefore a larger relative increase in saturated specific humidity. This motivates us to investigate the impact of warming on the precipitation rates from Alaska convective systems. Figure 10a shows the changes in the precipitation distribution for the summer season for different subregions in Alaska. There is a general intensification of precipitation by the end of the century with high precipitation rates over a certain threshold p contributing more to total precipitation, at the expense of the lower, more moderate precipitation rates. This intensification is well marked in all the subregions, but the threshold p varies between the subregions and is largest in interior Alaska. This transition seems to be a robust feature of climate change and has been documented in Europe (Kendon et al. 2014;Chan et al. 2018), Africa (Kendon et al. 2019), and the contiguous United States of America (CONUS) (Prein et al. 2017c;Rasmussen et al. 2017). Figure 10b focuses on precipitation changes in convective systems. The tracked storms mostly produce precipitation between 2 and 30 mm/h, with a peak contribution around 7 mm/h. This peak contribution stays similar in the future climates simulation, however, the future distribution has a higher skewness and its tail extends towards much higher, more extreme, precipitation rates. Therefore, precipitation rates above 10 mm/h contribute more to the storm precipitation in the future simulation, at the expense of the rates below this threshold, as shown by the red curve. It is important to mention that Fig. 10b shows normalized curves to highlight the structure of precipitation changes in future storms. The actual contribution of convective precipitation rates increases dramatically across the entire distribution driven by the large increase in convective storm frequency that was discussed in the previous section.
Other storm characteristics (lifetime, area and speed) do not exhibit any significant change between the control and PGW simulations, unlike the results of Prein et al. (2017b) for the contiguous United States (CONUS) that showed an increased rainfall area for future storms as well. Finally, the storm peak precipitation exhibits an intensification very similar to the precipitation coming from storms in general, Fig. 10 a Change in the fractional contribution of precipitation rates to total precipitation, for summer (MJJA) and different subregions as indicated in the legend. b Same as a, but only for the precipitation coming from storms. The shown statistics are based on the Analyzing Scales of Precipitation (ASoP) technique of (Klingaman et al. 2017) with exponential bins b n = 0.02e 0.16n , n ∈ [[10, 50]] (in mm/h). Shadings show the results of a bootstrap two-tailed 98% uncertainty test i.e. Fig. 10b curves are very similar but shifted to the right when done with storm peak precipitation only (not shown).
This increase in storm precipitation intensity can also be visualized with a storm composite, as shown on Fig. 11a. This figure shows a composite of average precipitation in the storms, in radial bands centered on the storm's maximum precipitation at each hour. Storms have a typical core radius of 10-20 km with precipitation rate decreasing sharply with distance within the core, and decreasing slower beyond. The intensity of precipitation is increasing in the core region of the storms ( ∼15 km radius) for precipitation rates larger ∼ 4 mm h −1 but this core region does not become larger. This result is not surprising since a robust increase in the intensity of extreme precipitating convective systems is found in many studies all over the world, whereas cell size changes are more uncertain. Indeed, Prein et al. (2017b) found an increase in the area of mesoscale convective systems in the US, whereas Peleg et al. (2018) found a decrease in the size of convective cells with temperature in Israel, pointing out that the response of storm area is dependent on the local climate, the morphology of the storms (e.g., single-cell storms v.s. organized convection), and the applied methodology (e.g., tracking algorithm and precipitation threshold). Figure 11b shows the proportion of precipitation coming from storms identified in this study among total precipitation exceeding a certain precipitation rate R min . In other words, with p(r) defined as in section 2.3.2 and in the limit → 0 it shows :

Contribution of storm precipitation to total precipitation
where p storm (r) is the precipitation rate distribution for storm-produced precipitation only. Therefore, the left-most values (at 0.1 mm h −1 ) show the proportion of total precipitation coming from storms. This proportion is about 1.2% in the historical simulation, and 5% in the PGW simulation. This increase is not surprising, given the fact that storms are both more numerous and more intense in the PGW simulation.
Storms have a relatively large contribution to extreme precipitation. For example, they account for about 15% of precipitation exceeding 5 mm/h and one-third of precipitation exceeding 10 mm/h. This proportion is significantly increasing in the PGW simulation. This means that precipitation from storms is becoming a much more important . 11 a Radially averaged hourly storm precipitation composite of precipitation, in the control and PGW simulations. The center point for the composite was chosen to be the location of maximum precipitation in each storm and each hour. b Percentage of precipitation coming from storms for precipitation exceeding a certain pre-cipitation rate only, in the region of interest and during the summer months. (Example : in the historical simulation, among all the model precipitation coming from rates above 7mm/h, 20% comes from storms.) Shadings show the result of a two-tailed 98% bootstrap uncertainty test component of the future hydrology of Alaska across a wide range of precipitation rates. Figure 11b also shows that extreme precipitation from storms increases significantly faster with global warming than extreme precipitation coming from other mechanisms, such as orographic enhancement or isolated, small convective cells. Figure 12 shows the mean and maximum precipitation of storms median, upper and lower quartiles binned by dewpoint temperature, calculated as : with standard notations, q being the water vapor mixing ratio, T 0 an arbitrary temperature and e 0 the saturated water vapor pressure at this temperature. Dewpoint temperature is a measure of moisture present in the air, and represents the temperature at which the air reaches saturation when it is cooled at constant pressure. Previous work has shown that it is beneficial to use dewpoint temperature instead of air temperature in precipitation scaling analyses (Lenderink and Attema 2015).

Precipitation scaling with temperature
In the historical simulations, the mean precipitation and maximum precipitation of the storms increase with dewpoint temperature, at a sub-CC (Clausius-Clapeyron) and CC scalings respectively. However, it has been shown that the temperature-precipitation relation does not necessarily follow a CC scaling as it can depend on the link between dewpoint temperature and large-scale atmospheric conditions (Lenderink et al. 2017), and it is not representative of the climate change signal (Prein et al. 2017c;Zhang et al. 2017). One should rather focus on the relative shift of this curve from the historical to the PGW simulation, as it is done in Prein et al. (2017c); Zhang et al. (2017); Drobinski et al. (2018). It appears that mean precipitation is scaling at sub-CC rates, while maximum precipitation scales close to CC expectations : the red curve is almost a copy of the blue curve shifted along the dashed lines in Fig. 12b. This shows that extreme precipitation from convective systems in Alaska does follow the increase in local atmospheric moisture, subsequent to the local increase in temperature. This CC scaling is in agreement with most kilometer-scale modeling studies in mid-latitude and tropical regions (Fowler et al. 2020).

Conclusion
In this study, we analyze the skill of convection-permitting climate simulations over Alaska in representing precipitation characteristics under present climate conditions. Additionally, we assess climate change impacts on convective systems by comparing storms under current conditions with storms at the end of the century under a business as usual emission scenario using the pseudo-global-warming technique. A tracking algorithm was applied in order to identify convective systems across Northern, Western, and Interior Alaska during the summer season (May-August). The model accurately reproduces hourly summertime precipitation statistics from observations, and the annual cycle of the identified convective systems agrees with the few available observations in Alaska, and other high latitude regions such as Finland or Siberia. However, the diurnal cycle of precipitation and of convective storms has a delay of 1-3 h relative to observations. We track large, organized convective systems that exist for several hours and reach a size of several dozens of square kilometers up to scales of mesoscale convective systems. These storms have maximum precipitation rates on the order of 30 mm/h and move slowly (on average 10 km h −1 ), which might produce flash flooding in mesoscale catchments. These storms occur mainly South of the Yukon River, with an average recurrence period of about once a year at each location. Storms in the Northern and Western parts of Alaska are much rarer.
Under future climate conditions the frequency of storms will triple. South of the Yukon River the storm frequency changes from once a year to almost once a month during summer. Additionally, convective systems start to appear in regions where they were not simulated in the historical run, especially in the North Slope and on the West Coast. Besides the dramatic increase in the frequency of these systems they also intensify in terms of their hourly rainfall rates by 5.5%. Highest increases are found in the core of the storms (their maximum precipitation increases by 37%) whereas their size does not change. This intensification in addition to the frequency increase results in an increased contribution of extreme precipitation to the total precipitation. The extreme hourly precipitation increase is in line with increasing atmospheric saturation vapor pressure following roughly a Clausius-Clapeyron scaling (i.e., 7% per degree Celsius). Due to these changes, convective systems will have a larger contribution to moderate and extreme precipitation in the future climate, making them the dominant process for extreme precipitation in Alaska. This implies a large increase in the potential magnitude and frequency of summertime flash floods in Alaska. Also, although we tracked systems producing heavy precipitation, our results raise the question of high latitude dry thunderstorms, that might strongly respond to climate change as well, potentially causing additional wildfire ignitions (Partain et al. 2016). This might have potentially large impacts on ecosystems and deserves further work.
The presented results mainly reflect the thermodynamic response of convective systems to climate change due to the applied PGW approach. We do not account for potential changes in synoptic circulation, which could alter these results. This study only accounts for changes in the monthly mean circulation, for example strengthening of the polar cell. In particular, it does not account for potential changes in the frequency of intrusion of transient eddies, moisture plumes or atmospheric rivers from the low latitudes, that could have a significant influence on convective storm dynamics. Also, changes in dynamical structures such as fronts or low level pressure systems that could favor large-scale ascent, and therefore help the onset of convection, are not taken into account if these features are not formed within the domain. Finally, changes in wind shear are likely not well represented in this approach, (and indeed no statistically significant change in the wind shear occurs in these simulations). This is still a large unknown of climate change, especially given the fact that most of the observed historical changes in the frequency of convective systems in other regions of the world has been attributed to circulation changes (e.g., changes in the monsoonal circulation in Bangladesh (Habib et al. 2019) or changes in dominant synoptic winds in the Sahel (Fitzpatrick et al. 2020).) Future kilometer-scale simulations are in progress that will focus on including changes in synoptic conditions and covering longer time periods that allow investigating effects of climate internal variability. This will give an opportunity to compare the pure thermodynamic results of this study to results accounting for dynamical changes as well, and to have an assessment of climate variability of convective storms and how this variability might change in the future.
We also highlight that these results have been obtained for a high emission scenario, that is unlikely to happen if current climate policies are strengthened in the coming decades (Rogelj et al. 2016). Climate change impacts by the end of the century might be much smaller if emissions are mitigated in the coming decades.
Our results differ from those reported in similar studies conducted in the mid-latitudes, that show only a small increase (Prein et al. 2017b) or even a decrease (Dai et al. 2017) in the number of convective systems, for regional convection-permitting simulations (Dai et al. 2017) or idealized large-eddy simulations (Lochbihler et al. 2019). In the mid-latitudes, more intense storms consuming more moisture, combined with limited moisture availability, are expected to result in fewer but more intense storms Dai et al. (2017);Fitzpatrick et al. (2020). This study shows that the response of deep convection might be very different for high latitude regions.
Future analyses will focus on the physical reasons for the increase in the frequency and the intensity of convective systems in Alaska to better understand the contribution of thermodynamic and dynamic changes, as well as the role of specific local mechanisms such as sea ice loss, or uneven warming across the state.