Trends of atmospheric circulation during singular hot days in Europe

The influence of climate change on mid-latitudes atmospheric circulation is still very uncertain. The large internal variability makes it difficult to extract any statistically significant signal regarding the evolution of the circulation. Here we propose a methodology to calculate dynamical trends tailored to the circulation of specific days by computing the evolution of the distances between the circulation of the day of interest and the other days of the time series. We compute these dynamical trends for two case studies of the hottest days recorded in two different European regions (corresponding to the heat-waves of summer 2003 and 2010). We use the NCEP reanalysis dataset, an ensemble of CMIP5 models, and a large ensemble of a single model (CESM), in order to account for different sources of uncertainty. While we find a positive trend for most models for 2003, we cannot conclude for 2010 since the models disagree on the trend estimates.


Introduction
Extreme event attribution (EEA)  aims at evaluating how the properties of a specific extreme climate event have been affected by anthropogenic forcings. Climate change may play a role on either-or both-the dynamics and the thermodynamics explaining the event. The influence of climate change on the thermodynamics of European heatwaves has been largely studied and proven for both specific events (e.g. Stott et al 2004, Christidis et al 2015, Russo et al 2015 and types of events (e.g. heatwaves in Russo et al 2014). The evolution of the dynamics related to heat-waves is still a debated subject.
The atmospheric dynamics in the Northern Hemisphere mid-latitudes are driven by the vertical static stability (e.g. Lim and Simmonds, Walland and Simmonds 1999) and by the latitudinal temperature gradient. This gradient could be modified by climate change through two processes: surface arctic amplification (AA) and upper-tropospheric tropical warming (Peings et al 2017). The evolution of those two factors is still very uncertain, with a wide range of responses across climate models (Zappa and Shepherd 2017), and even across different members of a single model ensemble due to internal variability (Deser et al 2017, Peings et al 2017. Over Europe, the link between long-lasting anticyclonic circulation, called blockings (e.g. Ruti et al (2014)), and high summer temperatures has been established (e.g. Jézéquel et al 2017a, Pfahl andWernli 2012, Sousa et al 2018). Francis and Vavrus (2012) detected the emergence of a significant increase in the persistence of blockings over the recent years using a reanalysis dataset. They explain this emergence by a mechanism based on the AA. Coumou et al (2015) found similar results focusing on summer and using satellite data. However, both Barnes (2013) and Screen and Simmonds (2013) argue that the results of Francis and Vavrus (2012) depend on the methodology they used and could be subject to ambiguous interpretations. Cattiaux et al (2016) used global climate models to extend the search of trends to the twentyfirst century. They found no evidence of an increase of persistence of blockings. Those studies evaluate the evolution of the circulation on large scales, on either the whole Northern Hemisphere or the North Atlantic region. In contrast, we are interested in capturing trends related to specific heatwave events, and we hence focus on a much smaller scale. Ruti et al (2014) calculated summer trends of the blocking index defined by Tibaldi and Molteni (1990) over the Euro-Russian region using a reanalysis dataset and an atmospheric-only model for the 20th century. They found a statistically significant increase in the duration of blocking episodes for the second part of the century, which they attribute to climate change, using different forcings as inputs of their model. However, the 20th century might not be long enough to evaluate trends on blockings. Indeed, using a large ensemble from a single model representing internal variability, Peings et al (2017) found a decrease in the blocking index over the 1920-2100 period for the North Atlantic region, which includes Ruti et al's Euro-Russian region. Those differences could be related to an inconsistency between different models or to different evaluations of the internal variability. This led us to use a set of different models and a large ensemble to account for both.
In the context of EEA, Trenberth et al (2015) argued that due to the large internal variability of dynamical processes, it is best to focus only on thermodynamical processes for a fixed dynamical state in order to extract the signal related to climate change. A few attribution studies that condition the signal to the circulation follow this approach to extract thermodynamical signals hidden in a large internal variability (e.g. Cattiaux et al 2010, Meredith et al 2015. However, this does not allow to calculate the complete influence of climate change on the events of interest . Shepherd (2016) highlighted that it is possible to study the dynamic and thermodynamic contributions separately. Few papers have studied the influence of climate change on the dynamics applied to a singular event (Vautard et al 2016, Yiou et al 2017. Both of those articles calculate the dynamical difference between two worlds (with and without climate change). Here we focus on detecting whether there is an evolution between 1950 and 2100 in the occurrence of circulations related to a given day.
Jézéquel et al (2017b) proposed to calculate a trend on the number of close days to the observed flow of December 2015 in Western Europe using a single model ensemble. In the present article, we refine this approach to single day atmospheric circulation patterns. We detail the proper statical methodology to calculate dynamical trends with a focus on the calculation of the statistical confidence interval, of multi-model uncertainties, and of internal variability. We seek to detect changes in the occurrence of circulation patterns related to specific hot days. We leave the attribution of those changes to further studies.
We first present the methodology to estimate trends of the circulation for a given daily event. We then apply this methodology to two case studies: the 2003 heatwave in Western Europe and the 2010 heatwave in Russia. These two heat-waves have been ranked first and second in Russo et al (2015) list of top ten European heat-waves since 1950. We finally discuss those findings and potential larger applications of our methodology to other types of events.

Datasets
In this study, we assume that the geopotential height at 500 hPa (Z500) is a proxy for the extra-tropical atmospheric circulation. We focus on the summer season (June-July-August: JJA). We use daily averages of Z500 from three datasets over two , called Western Europe (WE) hereafter and [10 E-68 E; 45 N-70 N], called Russia (RU) hereafter.
The first dataset is the National Centers for Environmental Prediction/National Center for Atmospheric Research, NCEP/NCAR, reanalysis I dataset (Kalnay et al 1996) between 1950 and 2016. Its horizontal resolution is 2.5 by 2.5 degrees. This dataset, called A 1 hereafter, allows us to assess whether dynamical trends are detectable in a short dataset, which is as close as possible to the observations. The second dataset is an ensemble of 18 models from the fifth Coupled Model Inter-comparison Project (CMIP5) ( Taylor  We use three types of data in order to compare reanalysis data with a single model ensemble (CESM-LENS) that reflects the internal variability of a climate model and a multi-model ensemble (CMIP5) that reflects the uncertainty due to the model formulation. This allows to estimate different components of the uncertainty (section 2.3).
Historical runs over 1950-2005 are merged with RCP8.5 runs over the 2006-2016 period to allow the comparison with reanalysis data over the whole 1950-2016 period. The choice of RCP8.5 is (1) coherent with In this article, we focus on very hot days, which are related to anticyclonic blocking situations. We are therefore interested in finding close Z500 patterns to those types of circulation. The Z500 is however related to lower-tropospheric temperatures, so that a global surface warming implies a generalized Z500 increase. In order to focus on the dynamical signal and ensure that our method would not interpret a uniform Z500 rise as a change in circulation, we choose to remove this background thermal effect (contrarily to Horton et al 2015). This way, we aim at dealing with dynamical changes unrelated to thermodynamical trends. This is done by subtracting a spatially uniform Z500 trend, calculated on the mean seasonal (JJA) spatial average on the region of interest, using a cubic smoothing spline in time (similarly to Jézéquel et al 2017a). By subtracting a uniform field, we do not alter the horizontal gradients of Z500 that depict the circulation. An alternative to using Z500 would have been to use SLP but in summer, the SLP field is affected by a heat low effect that blurs the dynamical signal (Jézéquel et al 2017a).

Dynamical trend estimation
Our goal is to determine whether a given circulation pattern has become more or less frequent during a given period. We consider a Z500 reference pattern Z belonging to the dataset A 1 that occurs on a day d. For all the days d ′ in the dataset A , we compute the set of Euclidean distances between ′ ∈ and the reference ∈ 1 , defined as the root mean square of the differences between each grid point within the region of interest. For the reanalysis dataset, we exclude the days within the same year as the event of interest. We determine the xth quantile q of those distances for each separate dataset A (the value of q can hence differ depending on the dataset). The value of x can be chosen heuristically, e.g. the 5th quantile. From Z and q we define the class of days or patterns D(Z ) in the ensemble A that are similar to Z : The class D(Z ) is shown for August 13th 2003 over the WE region for one model of A 2 (MPI-ESM-MR) in figure 1(a) (blue dots). figures S3-S8 in the supplementary material show that even if the exact anomaly of Z500 is not captured by the days in D(Z ), they all display blocking patterns within the regions of interest. This means that the 5th percentile chosen to define D(Z ) is relevant to study the evolution of blocking patterns in those regions.
For each year y in A , we count the number N ) in order to study potential trends in . This requires to properly model the evolution of this variable. The first step is to find a suitable distribution to describe it. The variable is discrete and bounded. can only take integer values between 0 and N tot = 92 (the number of days in JJA). We display the evolution of with time in figure 1(b) for one model. As Var( ) is 2.0-15.2 times larger than the expected value E( ), we conclude that the distribution of is systematically overdispersed with respect to a Poisson or to a binomial distribution (with parameter p), for which the variances would be respectively equal to E( ) and (1 − p)E( ).
Once there is one day in D( ) in a given summer, there is a high chance that the following days will also be in D( ), because of the persistence of atmospheric circulation. Hence the odds of having another day in D( ) within a given year increase with the number of days already in D( ) within the summer.This explains why is overdispersed. We chose to model the distribution as a beta-binomial distribution, which fits well bounded discrete distributions that are overdispersed, so that: where B is the beta function (Whittaker and Watson 1996), and and parameters which allow to account for possible overdispersion. We tested the goodness of fit of the beta-binomial distribution for each dataset using a Pearson 2 test. The p-values are all greater than the 0.05 significance level, meaning that we cannot reject the hypothesis that follows a beta-binomial distribution.
The second step is to find a statistical model to describe the evolution of with time. We used a generalized linear model (glm, see (equation 4)) to determine the temporal trend of (Nelder and Wedderburn 1972). The glm is a generalization of the linear regression through the use of a link function g allowing the transformed mean to vary as a function of predictors. We transform the mean as g(E( /N tot )) where with u ∈ [0, 1] and E(.) is the expected value. g is called the logit link function. We used the R package VGAM (Yee 2010), which includes the function vglm that fits a glm to betabinomial distributions (Prentice 1986).
For a year y in A , we assume that where and are the regression coefficients. The interpretation of regression coefficients is not straightforward because the glm uses the logit link function, which produces a non-linear regression. We therefore present the results using fitted values of E( ). We used the inverse link function ( ) = × −1 ( + ) and the regression coefficients to obtain the fitted values of E( ) for year y, which gives the solid red line in figure 1(b). We then calculated the difference between the fitted values of E( ) between the end and the beginning of the time series, in order to analyze the evolution of E( ).
This regression is a way to determine whether the days similar to get more (or less) likely with time. However, it does not discriminate whether any change detected is related to the fact that days close to happen more regularly every summer, or if they are more numerous within a given event. Decomposing those two parts of the signal is beyond the scope of the present article.

Uncertainties
In order to derive a confidence interval on the estimated trend, we first calculated a confidence interval for -this is done assuming that̂follows a Gaussian distribution. This confidence interval on can then be translated into a confidence interval on the average number of days belonging to D( ), by calculating the fitted values of E( ) corresponding to the upper (resp.lower) bound of . We consider that the change is significant if the confidence interval on does not include 0.
Besides the statistical uncertainty, the two ensemble datasets allow to evaluate the uncertainty due to internal variability in the case of CESM-LENS A 3 and the multi-model uncertainty in the case of the CMIP5 ensemble A 2 .
The comparison of those three sources of uncertainties allows us to detect whether the circulation undergoes a significant evolution. It also weighs the sources of uncertainties and assesses the confidence in the methodology. We cannot attribute any detected evolution to climate change with this methodology, as we do not compare our results to those which could be obtained in a world without climate change.

Two case studies
We chose two epitomes of heat-waves of the 21st century, largely studied in the literature to apply our method: summer 2003 (e.g. Beniston 2004, Fischer et al 2007, Stéfanon et al 2012 in the WE region and summer 2010 (e.g. Dole et al 2011, Rahmstorf and Coumou 2011, Trenberth and Fasullo 2012, Otto et al 2012, Hauser et al 2016 in the RU region. The thermodynamical component of climate change has been identified by those authors, but the dynamical contribution has not been as emphasized. We used those two cases as examples to apply our methodology to detect circulation trends. The hottest day of the NCEP reanalyses in the WE region was recorded on August 13th 2003, and the hottest day in the RU region was recorded on August   In the case of August 13th 2003, CanESM2 and three runs of CESM-LENS have significant positive trends, and one run of CESM-LENS has a significant negative trend from 1950-2016. The other models and runs display no significant trend. The bigger uncertainty comes from the internal variability assessed with CESM-LENS. This means that we cannot judge the quality of a model with respect to the simulation of dynamical trends by comparing it to the NCEP reanalysis, which is just one realization of what could have happened for the same background state of the climate. In the case of August 7th 2010, no model detects either a positive or a negative significant trend on the historical period. The statistical uncertainty is larger than for 2003. The multi-model uncertainty equals the internal variability. Using only reanalyses or historical runs of 67 years is not sufficient to detect any significant signal. This is coherent with the findings of Deser et al (2017) who have shown that SLP trends over the North Atlantic region have different signs for different runs of CESM-LENS even over 50 years, although the focus of their study was the winter season. We get past the internal variability using 151 years (from 1950-2100) and RCP scenarios.
For the longer periods, the results differ between 2003 and 2010. For the former and RCP 4.5, seven models detect a significant positive signal. For RCP 8.5, ten models detect a significant positive signal. Out of the 30 runs of CESM-LENS, 29 detect significant positive difference between 1950 and 2100. With the exception of MIROC models, the models which detect a significant positive trend reproduce best the observed anomaly ( figure S5). Although the response differs from one model to another, there seems to be an agreement on a positive difference of approximately 5 days in 151 years. With the choice of the 5th percentile to define D( ), the mean number of days in D( ) for each summer is approximately 4 days. Therefore a difference of 5 days is not negligible. The models do not agree for 2010. For RCP4.5, we find two models with a significant positive trend. For RCP8.5, we find four models with significant positive trends, and three models with significant negative trends. Out of the 30 runs of CESM-LENS, 27 yield a significantly positive trend. The models hence disagree, which questions the robustness of trends found in studies where only one model is used. The models that find a significant positive trend (including CESM-LENS) are less able to reproduce the intensity of the observed Z500 anomaly (supplementary figures S4 and S6).

Discussion
Our methodology gives different results for the 2003 and 2010 events. While we find a positive signal with most models for 2003, the models do not show coherence for 2010. This is not surprising, as both events happened in two different regions, in which there is no reason for the dynamics to evolve in a similar way. We can see in figures 2(c) and (d) that the atmospheric pattern in Western Europe in 2010 is almost the inverse of the 2003 pattern. However, we are more confident in the ability of the models to reproduce the 2003 pattern than the 2010 one because of the difference in intensity and extent of both blockings and the larger spread in figure 3(b) compared to 3a. The trends are more pronounced in RCP8.5 than RCP4.5, which is an argument to attribute the significant changes in the weather pattern of one central day of the 2003 heatwave to climate change. If the models with significant positive trends are to be believed, and for the 2003 case these models are the ones simulating the most realistic patterns (see supplementary figures S5 and S7), this could mean longer and more frequent heatwaves similar to 2003 in Western Europe, without even taking into account the thermodynamical effect of climate change on temperature. This thermodynamical effect has been largely proven in the literature (e.g. Meehl and Tebaldi (2004), Bador et al (2017)), and is stronger than the dynamical effect. We however stress that the Z500 anomaly is not a sufficient condition for a heat-wave to develop (Boschat et al (2016), Quesada et al (2012)). Peings et al (2017) find a decrease in the one-dimensional blocking index as defined by Tibaldi and Molteni (1990), which would indicate a lesser importance of the dynamics in the years to come, using the CESM-LENS dataset. There is no reason to expect the same results from both studies, since we focus on a specific dynamical event through the use of a two-dimensional Z500 field over a rather small region, while Peings et al (2017) looked at circulations leading to heat-waves in general over a much larger region.
All the Z500 fields were detrended to remove from Z500 the thermodynamical influence of climate change. However, the shape of the modeled Z500 distribution can differ from the observed one. We tested four types of normalization: no normalization, a simple normalization (division by the standard deviation) on every grid-point, a simple normalization on the mean of the Z500 field and a quantile-mapping (e.g. Panofsky and Brier 1958, Déqué 2007, and Gudmundsson et al 2012. We normalized using the 1950-2005 period which is common between historical runs and NCEP. Although the normalization changes results for a few individual models, it does not change the collective results of the ensemble of CMIP5 and CESM-LENS models (not shown here). Since the normalization does not fundamentally change our results, we use non normalized Z500 anomaly fields.
We also tested how the results change when we choose a different percentile to define D( ). We tested four percentiles: the 2nd, the 5th, the 10th and the 25th percentiles. The differences detected between the 1950 and 2100 values of monotonically increase with the percentile. The results get more significant (further from 0 and in some cases become significant) for higher percentiles.
There are a few limitations to this methodology. We only considered daily events, which are not the heat events with the largest impacts. In the supplementary material, we calculated the dynamical trends for each day of both events (figures S1 and S2). In terms of dynamical trends, we find that August 13th 2003 and August 8th 2010 are typical of the whole heatwaves. We also observe that for both cases RCP4.5 has less statistically significant models than RCP8.5 which could mean that the dynamical signal is enhanced with a stronger climate change. Another caveat is related to the internal variability of the dynamics. Given that 70 years are not enough for any signal to exceed the range of observed natural variability, we have to rely heavily on models that might not accurately reproduce some aspects of the dynamics of the atmosphere.
The biggest advantage of this methodology is that it is easy to implement and very cheap in computation time. It would be possible to do those calculations in a few minutes time each day for a region of interest, and hence give an idea of whether climate change might make dynamically driven events more or less likely in the future for very specific types of circulation. It could serve for other types of events than hot days, e.g. for atmospheric patterns leading to daily extreme precipitations. In further studies, we intend to use it more systematically to see if it helps us to identify types of circulation whose probabilities evolve according to an ensemble of models.