Variability of daily winter wind speed distribution over Northern Europe during the past millennium in regional and global climate simulations

. We analyse the variability of the probability distribution of daily wind speed in wintertime over Northern and Central Europe in a series of global and regional climate simulations covering the last centuries, and in reanalysis products covering approximately the last 60 years. The focus of the study lies on identifying the link of the variations in the wind speed distribution to the regional near-surface temperature, to the meridional temperature gradient and to the North Atlantic Oscillation. Our main result is that the link between the daily wind distribution and the regional climate drivers is strongly model dependent. The global models tend to behave similarly, although they show some discrepancies. The two regional models also tend to behave similarly to each other, but surprisingly the results derived from each regional model strongly deviates from the results derived from its driving global model. In addition, considering multi-centennial timescales, we ﬁnd


Introduction
Anthropogenic climate change is expected to cause an increase of various types of extreme events, such as heatwaves, but its effects on extreme winds is less clear. Section 3 of the Intergovernmental Panel on Climate Change (IPCC) special report "Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation" states that there is only low confidence in projections of changes in extreme winds (Seneviratne et al., 2012). One way to reduce this uncertainty is to compare the output of paleoclimate simulations over the past centuries with empirical evidence of past wind conditions, for instance derived from historical evidence or natural proxies (Costas, 2013). While there is still a dearth of proxy records reflecting past changes in wind speed, new types of proxy records are being developed (Costas, 2013). A precondition for this comparison is to test whether different climate models provide a consistent picture of past changes in wind speed distribution. In this study we analyse several simulations with global and regional models and investigate to what extent they provide a consistent picture of the relationship between the variations in the wind speed distribution and large-scale drivers. We focus on Northern Europe in wintertime as this region and season are particularly prone to storminess.
Published by Copernicus Publications on behalf of the European Geosciences Union. 318 S. E. Bierstedt et al.: Variability of wind speed distribution during the past millennium The hypotheses put forward to explain changes in storminess are related to the general physical consideration that warmer periods provide more humidity and consequently more (latent) energy for possible storms. However, warmer periods are generally characterized by a weaker meridional temperature gradient due to the stronger warming of the high latitudes with respect to the tropics, and thus a weaker baroclinicity, which should lead to weaker or less storms (Li and Woollings, 2014;Yin, 2005). In addition, the North Atlantic Oscillation (NAO), as the main pattern of troposphere dynamics over the North Atlantic-European sector, is also related to the interannual variability of seasonal mean winds in this region. It is thus very plausible that the NAO is also involved in the variations of the distribution of daily wind speeds. For example, Wang et al. (2011) stated that NAO variations show a relationship with the 99th percentile of wind speed. The NAO itself could also be related to changes in European near-surface winter temperature (Rutgersson et al., 2015). In this regard, scenario simulations indicate a contradicting tendency of the NAO in a warmer future, depending on whether it is defined from the sea-level pressure (SLP) gradient or from the geopotential height at 500 mbar (Z500). Gillett and Fyfe (2013) showed that, in climate simulations, the meridional SLP gradient will tend to become steeper polewards. On the other hand Cattiaux and Cassou (2013) also found a positive trend of the SLP gradient but found a negative trend of the Z500 gradient. Our study focuses on surface winds, which are rather influenced by SLP variations, and therefore we focus on the link between SLP and daily wind distribution. Regarding the climate of the past centuries, the link between the NAO and external climate forcing is, however, unclear in climate simulations (Gómez-Navarro and Zorita, 2013).
Thus, for Northern Europe, from the dynamical point of view it is not clear how the distribution of wind speed would respond to changes in temperature. The analysis of long-term trends in wind extremes and storminess in the observational record has so far yielded inconclusive results, probably due to the difficulty of constructing homogeneous series of wind speed, because of e.g. station relocation or changing measuring techniques. Furthermore, the covered period might be too short to realistically demonstrate trends in the rarely occurring extreme wind events. On the other hand, reanalysis products covering long periods (e.g. the 20th Century Reanalysis 20CR; Compo et al., 2011) may be inhomogeneous due to the assimilation of different types of data through the simulated period. There has been a considerable debate on the storminess trends in long-term reanalysis data sets (Brönnimann et al., 2012;Krueger et al., 2013;Wang et al., 2013). Krueger et al. (2013) stated that before 1950 the time series of 20CR and observational mean sea level pressure are not consistent. They suggested that the increasing density of station data leads to these inconsistencies. Wang et al. (2013), on the other hand, argued that new measurement errors and changing sampling frequencies would produce these inho-mogeneities. This debate has also been discussed by Feser et al. (2015), who concluded that the analysis of long-term reanalysis data, affected by changing station density, should be conducted with caution.
The analysis of the climate of the past centuries can shed light on the question of whether external climate forcing has an effect on the intensity or frequency of wind extremes and whether or not the temperature variability is linked to variability in statistics of wind speeds. Unfortunately, proxybased climate reconstructions in general still do not provide information about extreme wind statistics in the past, except for intense tropical cyclones (e.g. Donnelly and Woodruff, 2007). However, new types of proxy data that may record past wind speed are being retrieved. For instance, coastal dunes at the North Sea coast contain layered structures that can be analysed by ground-penetrating radar. The layered structures contain information about grain size distribution and, indirectly, about intensity or frequency of high winds in the past (Costas, 2013). These types of proxies can potentially be used to test the skill of climate models to simulate the relationships between extreme wind statistics and external forcing or between extreme wind statistics and lowfrequency variability of the large-scale surface climate.
The evolution of temperatures of the past millennium in this region, as reconstructed from proxy and longinstrumental records, exhibits a generally warm period in the early centuries (the Medieval Warm Period) and generally colder centuries around 1700 AD (the Little Ice Age), with the subsequent warming leading to the current warm period PAGES 2k, 2013;Esper et al., 2014). This alternation of warm and cold periods has been likely caused by external climate forcings (Hunt, 2006;Fernandez-Donado et al., 2013), in particular the recent warm period and the Little Ice Age. Therefore this period provides a suitable test of whether the variability of extreme wind statistics may follow a similar alternation.
Climate simulations had been previously used to address the connection between winds and temperatures in the past (e.g. Fischer-Bruns et al., 2005;Schimanke et al., 2012). Fischer-Bruns et al. (2005) analysed two historical climate simulations with the global climate model ECHO-G covering the period from 1550 to 1990 AD. These authors found that storminess and large-scale temperature variations were mostly decoupled in these simulations. However, they reported a connection between storm track variability and temperature over the North Atlantic for one of the two simulations in two periods with extremely low external forcing, namely the Late Maunder Minimum (1675-1710 AD) and the Dalton Minimum (1790-1840 AD). They found no evidence of a linear co-variation between the number of extratropical storms and temperature variations in the simulations analysed.
The spatial resolution of global climate models may not be adequate to realistically represent extreme events, especially over regions with complex coastlines. In this respect, Table 1. Overview of the analysed simulations/reanalysis and their simulation acronyms, underlying atmosphere and ocean models, boundary forcings (only for regional data sets) as well as the spatial resolution of the atmosphere models and time periods, as used for the analysis. regional climate models, driven by the fields simulated by global climate models, should provide a better representation of small-scale processes, topographic influences and of the land-sea contrasts (Gómez-Navarro et al., 2015b;Hall, 2014;Gómez-Navarro et al., 2013), and thus they should be better suited for the simulation of extreme events. Nevertheless, despite the fact that regional models provide an added value (Feser et al., 2011) they are also bound by the circulation biases of the driving global climate model simulations (Hall, 2014). In this study we present an analysis of the variability of daily wind speed statistics over Northern Europe over the past centuries as simulated by different regional and global climate models. We mainly focus on the consistency among the different models in simulating the relationship between large-scale drivers and the statistics of daily wind speed with the goal of identifying robust patterns across models that can be later tested with proxy reconstructions. These results are also compared to a similar analysis of reanalysis data sets. Although these data sets cover a shorter period and, therefore, they cannot properly capture the decadal and multidecadal variability, they at least offer a possibility to ground truth the results obtained from free-running climate simulations.
This paper is structured as follows: Sect. 2 describes the analysed data sets separated into climate simulations of global circulation models, regional circulation models and reanalysis products. Section 3 defines our area of interest and outlines the applied methods and definitions. Section 4 presents the analysis of the relationship of large-scale drivers and wind speed variance, as well as the comparison of the evolution of the wind speed variance in the millennium simulations. A discussion of the results and conclusions closes the manuscript.

Data
Our study focuses on the statistical relationship between spatial and temporal mean temperature/pressure and daily wind statistics. We use monthly mean 2 m temperature (T2M) val-ues, monthly mean values of mean sea level pressure (MSLP) and daily mean 10 m wind speed (WS) for our analysis. These values are taken from a set of five simulations performed with five different models, with different spatial and temporal resolution. These models include global as well as regional models. Additionally, we analyse one global reanalysis and one regional reanalysis. Table 1 summarizes the information about spatial and temporal resolution and time periods of the data sets used. Figure 1 shows the land-sea-masks for all data sets. For the regional data sets the whole available domain is visible and the study area is marked by a red rectangle. In the following we introduce the different data sets in more detail and describe their main differences.

Global climate model simulations
The coupled GCM ECHAM4-HOPE-G, also denoted in previous literature as ECHO-G (Legutke and Voss, 1999) consists of the atmospheric component (AGCM) ECHAM4 (Roeckner et al., 1996) and the ocean-ice component (OGCM) HOPE-G (Wolff et al., 1997). Both sub-models have been developed at the Max Planck Institute for Meteorology (MPI-M) in Hamburg. The ECHO-G millennium simulation (1001-1990 AD) is part of an ensemble of simulations conducted by the Helmholtz-Zentrum Geesthacht (HZG). The simulation is forced with changes in total solar irradiance, the dimming effect of volcanic eruptions on solar irradiance, and changes in greenhouse gases (CO 2 , NO 2 , CH 4 ). The simulation was started with a cold ocean initial condition taken from a previous simulation corresponding to a situation representative of the Little Ice Age (Hünicke et al., 2011;Zorita et al., 2004).
The coupled GCM ECHAM5/MPI-OM consists of the atmospheric component ECHAM5 (successor of ECHAM4, Roeckner et al., 2003) and the ocean and sea-ice component MPI-OM (Marsland et al., 2003). The analysed simulation ECHAM5/MPI-OM (850-2005 AD) is part of an ensemble of simulations conducted by the Max Planck Institute for Meteorology in Hamburg and will hereafter simply be named ECHAM5. It is driven by changes in solar irradiance (Bard et al., 2000), volcanic eruptions (Crowley et al.,  We also include a climate simulation with the model MPI-ESM with the P configuration (MPI-ESM-P, Giorgetta et al., 2013), which consist of the successor of ECHAM5 and the newest version of ECHAM, named ECHAM6 (850-2005 AD). ECHAM6/MPI-OM was chosen due to its higher spatial resolution similar to the resolution of the RCM simulations analysed, but also for being the next generation of the previous versions of ECHAM4 and ECHAM5. The main differences between ECHAM6 and ECHAM5 are the higher vertical resolution (i.e., 47 instead of 31 vertical levels), increased horizontal resolution, the incorporation of new aerosol and surface albedo climatology and the use of a new shortwave radiation scheme in ECHAM6 . The MPI-ESM-P simulation is part of the Climate Model Intercomparsion Project version 5 (CMIP5, Taylor et al., 2012). The boundary layer and turbulence parametrization in ECHAM6 is based on the eddy diffusivity/viscosity approach . The model was driven by changes in greenhouse gases and spectrally resolved solar irradiance, volcanic activity (Crowley et al., 2008), changes in Earth's orbital parameters and land use changes (Pongratz et al., 2008).
We outlined above the key-properties of the analysed global simulations. However, the differences in the external forcings used to drive the simulations play a crucial role in our analysis. Hence, we additionally provide a comparison of these differences. An overview can be found in Table 2. A more comprehensive comparison between ECHO-G and ECHAM5 was given by (Fernandez-Donado et al., 2013), which we here summarize and extend for ECHAM6.
While all three simulations include total solar irradiance (TSI), greenhouse gas (GHG) and volcanic forcing as external forcings, only ECHAM5 and ECHAM6 incorporate also anthropogenic aerosols and land use changes. There are various estimations of past TSI, which can be broadly divided into a strong (S; > 0.2 % TSI change since the Late Maunder Minimum (LMM)) and a weak (s; < 0.1 % TSI  Schmidt et al. (2011) change since LMM) amplitude of variations (Fernandez-Donado et al., 2013). Strong (S) solar forcing is applied for ECHO-G and ECHAM5, in which ECHO-G uses higher values than ECHAM5. Weak (s) solar forcing is applied for ECHAM6. Estimations of changes in the main well-mixed GHG concentrations (CO 2 , CH 4 , N 2 O) are obtained from air bubbles enclosed in Antarctic ice cores. The CO 2 concentrations were prescribed in ECHO-G after Etheridge et al. (1996), in ECHAM5 and ECHAM6 they follow Marland et al. (2003). The incorporation of the volcanic forcing into ECHO-G was done by including the net effect of stratospheric volcanic aerosols in an effective solar constant in terms of a reduction in incoming shortwave radiation. For ECHAM5, latitudinally resolved changes in optical depth in the stratosphere were used. This might have an impact on climatic changes due to volcanic eruptions e.g. on the atmospheric circulation, especially over the extratropics during winter (see Fernandez-Donado et al. (2013) and included references). The orbital changes and land-use changes (only included in ECHAM5 and ECHAM6) are the same in both simulations. Orbital variations follow Bretagnon and Francou (1988) and vegetation changes Pongratz et al. (2008). Note that for the RCM simulations the same external forcings were applied as for their driving global simulations.

Regional climate model simulations
The RCM MM5 model consists of a slight modification of the non-hydrostatic Fifth-generation Pennsylvania-State University-National Center for Atmospheric Research Mesoscale Model. Such modification allows this meteorological model to perform long climate simulations. This setup has been used to conduct a long high-resolution climate simulation of the European climate during the last millennium, driven at the domain boundaries by the coupled GCM ECHO-G (Gómez-Navarro et al., 2013, 2015a. For the planetary boundary layer formation parametrization, this simulation used the medium-range forecast (MRF, Hong and Pan, 1996) scheme. The RCM was driven by the same set of external forcing as the driving GCM ECHO-G (Sect. 2.1) to avoid physical inconsistencies. The model output is available on a daily scale and covers the period 1001-1990 AD. The analysed millennium simulation MM5-ECHO-G will hereafter be named MM5. A second regional simulation was carried out with the nonhydrostatic operational weather prediction model COSMO in CLimate Mode (CCLM) (Rockel and Hense, 2008). The CCLM model was driven by initial and boundary conditions of the global ECHAM5 simulation. This regional model was free to produce its own small-scale spatial variability. The COSMO model uses a boundary layer approximation by implying horizontal homogeneity of variables and fluxes, which ignores all horizontal turbulent fluxes (Doms et al., 2011). The CCLM simulation over Europe had roughly the same spatial resolution as the MM5 simulation, but it covers only the period from 1655 to 1999 AD. The simulation CCLM-ECHAM5 will hereafter be abbreviated as CCLM.

Reanalysis data
The NCEP/NCAR reanalysis covers the period from 1948present and is available at 6-hourly intervals (Kalnay et al., 1996;Kistler et al., 2001). In addition we analyse the highresolution regional product coastDat2 (Geyer, 2014) resulting from a simulation with the regional model CCLM and driven by the global NCEP/NCAR reanalysis using a spectral nudging technique (after von Storch et al., 2000). The regional reanalysis coastDat2 covers Europe and the North Atlantic for the period from 1948 to present.
All wind speed data were daily averaged to proceed with the analysis.

Methods and definitions
Our analysis concentrates on the distribution of daily wind speed in wintertime (December, January, February -DJF) over Central and Northern Europe. The area of investigation has approximately the same extension from 45 to 65 • N and 0 to 30 • E for all data sets analysed.
The statistics of daily wind speed were evaluated over gliding time windows for the different simulation periods. These wind speed statistics include the standard deviation (SD) of the distribution, its 50th, 95th and 99th percentiles (P50, P95, P99) and the differences P95 minus P50 (diffM) and P99 minus P95 (diffE) as a measure of the width of the distribution in the high wind ranges. The analysis of several percentiles and their differences allows the determination of basic changes in the characteristics of wind speed distributions, hence it is possible to investigate if it shifts with time with unchanged shape and/or whether its width changes. Thus, we can discriminate between a change in the mean and/or in the extreme of the wind speed distribution. In our case "extreme" means the tail of the distribution, which includes values above the 95th percentile. For instance, increasing diffM and diffE values would show a broadening of the distribution which means higher extreme wind speeds. The three climate parameters analysed regarding their influence on wind speed are (1) mean seasonal near-surface air temperature (mTemp), (2) mean seasonal meridional temperature gradient (tGrad) and (3) the North Atlantic Oscillation index (NAO). Table 3 presents a summary of these statistical relationships derived from the different model simulations analysed. The presented time correlation coefficients are obtained by calculating the parameters of the daily wind probability distribution at gridcell scale, followed by averaging over the whole spatial domain.
The temperature gradient is calculated as the absolute value of the difference between the northern (N) and the southern (S) half of the investigation area tGrad = abs(N-S) for each model simulation. Due to the different resolutions the exact geographical domains of N and S differ. Therefore the border between both varies from around 53 to 57 • N.
The North Atlantic Oscillation (NAO) index is defined as the leading pattern resulting from principal component analysis (PCA) of the winter mean sea-level pressure (MSLP) field. This dominant pattern of variability is characterized by a low pressure system over Iceland and a high pressure system over the Azores (exemplarily shown for ECHAM6 in Fig. 2). As MSLP field we used the GCM gridded pressure information of the North Atlantic and European area. This domain approximately spans from 70 • W to 30 • E longitude (from ≈ Boston to Istanbul) and from 70 to 20 • N latitude (from ≈ Tromsø to the southern part of Morocco). For computations conducted for the RCM simulations we used the NAO patterns derived from the driving GCM fields as well.
Because we are interested in the relationship between the slowly changing mean climate and the variability of the distribution of daily wind speed, the wind statistics are calculated considering gliding time windows over the respective time series for each model simulation before they are correlated with the running mean values of the atmospheric drivers. The climate parameters analysed are considered as means over the respective time windows and for Figs. 3-8 these values are spatially averaged before the correlation computation with the wind statistics. For the long climate model simulations (ECHO-G, ECHAM5, ECHAM6, MM5, CCLM) we use 30-year running windows and for the shorter reanalysis data (coastDat2, NCEP) 5-year running time windows.

Test for significance: random-phase bootstrap
A random-phase bootstrap method (Schreiber and Schmitz, 1996;Ebisuzaki, 1997) is applied to determine the significance of the correlation coefficients shown in Table 3 with a significance level of p = 0.05. This method allows us to take into account the autocorrelation structure of the series. For this method a Fourier transformation of the time series is conducted. The phases of the Fourier-transformed series are then replaced by random phases, and the result is transformed back to the time domain to obtain new surrogate time series. The surrogate time series has the same spectrum and auto- Table 3. Time correlation coefficients between the following parameters of the probability distribution of daily mean wind speed: standard deviation of wind speed (SD), the 50th, 95th and 99th percentile (P50, P95, P99) and the differences between P95-P50 (diffM) and P99-P95 (diffE) and some large-scale drivers: spatially averaged December-February air temperature (mTemp), the spatial air temperature gradient (tGrad) and the North Atlantic Oscillation index (NAO). The parameters of the probability distributions have been computed in 30-year sliding windows for the simulations and in 5-year sliding windows for the reanalysis products. The time series of the drivers have been smoothed with a running mean filter. Significant (p < 0.05) coefficients (tested with a random phased bootstrap method) are written in bold. correlation as the original time series, but has a random time evolution. By generating a large number of surrogate time series, an empirical distribution of the correlation coefficient under the null-hypothesis (that the correlation is zero) can be constructed and used to determine the statistical significance of the correlation coefficient.

Results
In the following, we first present the general findings for each of the climate drivers analysed (mTemp, tGrad, NAO) by comparing all considered model simulations and reanalysis products. This is followed by a more detailed presentation of the results with a focus on (a) the regional model simulations and their corresponding driving global models (b) the simulation with the global model ECHAM6/MPI-OM and (c) the reanalysis products. In addition, (d) we compare the results for the overlapping time periods without reanalysis (1655-1990 AD) and with reanalysis  AD) data and (e) for some of the long simulations a comparison of time slices is presented.

Relationship between mTemp and tGrad
A common characteristic shared by all analysed simulations is the negative correlation between the mean winter temperature (mTemp) and the mean meridional temperature gradient (tGrad). Hence, in warmer decades the northern regions warm more strongly than the southern regions, and in colder decades the northern regions also cool more strongly than the southern regions. This "high-latitude amplification" is also found in climate simulations for future scenarios. In those simulations, it is caused by several positive feedbacks that operate more strongly at high latitudes, such as ice-snow-albedo feedback (Hall and Qu, 2006). In addition, the (negative) black-body radiation feedback is weaker at high latitudes so that the associated temperature response is stronger (Pithan and Mauritsen, 2014). European temperature reconstructions over the past centuries also indicate that in colder periods such as the Late Maunder Minimum (around 1700 AD) temperatures in higher latitudes cooled down more strongly than further south (Luterbacher et al., 2002).

Relationship between mTemp and wind speeds
A positive link between these two variables would support the idea that in a warmer atmosphere, holding more humid-ity and being more energetic in general, stronger winds are more probable. In fact, the relationship between the mTemp and the median winds (P50) is positive in all analysed simulations and reanalysis products, but with the exception of the regional simulation with MM5. However, the correlations, taken individually, are not always statistically significant at the 5 % level. Warmer air temperatures are also strongly linked to larger values of the high percentiles of the distribution of daily wind, P95 and P99, for most of the simulations. Again, the exceptions relate to the regional model simulations MM5 and CCLM. MM5 presents a negative correlation and CCLM a weak positive correlation.
Variations in the width of the daily wind distribution are described by the differences between the high percentiles, P95 or P99, and the median wind P50. The correlations between mean temperature and these measures of the distribution widths tend to be small for all simulations with the exceptions of the regional models MM5 and CCLM. For these two regional models the correlations are strongly negative, and more strongly so for the MM5 model, indicating that in periods with warmer air temperatures the wind distribution gets narrower at the same time that it shifts to lower values of wind speed, as indicated by the negative correlation with P50.
Briefly summarized, the relationship between winds and mean regional temperature in the global models does support the idea that warmer temperatures are associated with a shift of the distribution of daily winds as a whole. However, this link is not very strong and is contradicted by the regional models.

Relationship between tGrad and wind speeds
The temperature gradient should modulate the atmospheric baroclinicity, which should be reflected in the distribution of wind speed. The correlation coefficients between the distribution of wind speeds and tGrad are summarized in the third block of Table 3. In general, the correlations tend to be weak, with some exceptions. In the MM5 simulations they are stronger and positive, whereas in the ECHAM6 simulation they are somewhat weaker but negative. In the NCEP reanalysis the correlations between tGrad and the median wind P50 or the higher percentile winds P95 and P99 are negative and statistically significant. Therefore, this analysis does not support the idea of a relationship between stronger temperature gradient should cause stronger mean winds or more frequent extremes.

Relationship between NAO and wind speeds
The NAO is a large-scale winter circulation pattern that describes the mean strength of the seasonal mean westerly winds in the North Atlantic-European sector and therefore it is plausible that it is also related to the distribution of the daily wind speed in Northern Europe. The correlations between the NAO index across the different simulations yield, however, an incoherent picture. Most simulations do display a positive and relatively strong correlation between the NAO index and the spatially averaged P50, but again with the exception of the two regional models, MM5 and CCLM.
Thus, the regional models behave differently to their respective driving GCMs. In the case of MM5 the correlation between the NAO index and P50 is strikingly negative whereas in the case of CCLM the correlation is weakly positive. A positive phase of the NAO is linked to stronger westerly winds over Northern Europe and hence a negative or weakly negative correlation of the NAO with P50 is surprising. We show later that the negative sign of this correlation in the regional simulations can be explained by the behaviour of the regional models over land areas, whereas the sign of the correlation between NAO and wind over the ocean is the expected one.
The correlation between the NAO index and the width of the distribution (SD) of wind speed averaged over the study region tends to be also positive for most simulations, indicating that stronger mean westerlies tend to concur with a wider distribution of daily wind speed. However, there are exceptions. Again, the regional model simulation MM5 displays a strong negative correlation and the regional model simulation CCLM shows a positive but weak correlation. These negative (MM5) or positive but weak (CCLM) correlations also contrast with the link between the NAO index and the width of the wind speed distribution in their parent global models, ECHO-G and ECHAM5, respectively, both of which display positive and statistically significant correlations. Similarly to the global models, in both reanalysis products the NAO index is strongly and positively correlated with the width of the wind speed distribution.
Briefly stated, the global models display the expected link between the NAO and the median winds over the ocean, but this link is distorted over land. The regional models, with a domain mostly located over land, now show the expected relationship between NAO and wind speed.

Relationship between NAO and mTemp
It is well known that the winter NAO index is positively correlated with air temperatures in Northern Europe. The link between the parameters of the wind speed distribution on one side, and the NAO or the mean air temperature on the other side may thus be just a reflection of the same physical relationship. This is also supported by paying attention to how the correlations with the NAO and with the mean temperature vary across simulations (third line in Table 3). It seems clear that this line in the table displays a similar, though not identical, pattern of correlations across the simulations analysed. However, the spatially aggregated analysis does not allow to disentangle which one of both factors, NAO or mTemp, is the physical driving factor for the variations in the distribution of wind speed.
The results of this section can be summarized in two main points. All models reproduce the link between mean temperature and temperature gradient, and therefore the regional analysis seems to be physically consistent regarding the thermal variables. The link between these thermal variables and the distribution of wind speed is much weaker and very much model dependent, with regional models deviating from their respective global models.

Simulations of the regional models CCLM and MM5
versus their driving global models ECHAM5 and ECHO-G In the following subsections we investigate in more detail the links between the large-scale atmospheric drivers and the distributions of daily wind at a grid-cell level, which allows us to better understand the spatially aggregated correlations included in Table 3. The following figures (Figs. 3-6) display the correlation patterns between the different large-scale climate indices and the parameters of the wind speed distribution at grid-cell scale. The upper two panels in the figures show the results derived from the regional models for direct comparison with their driving GCMs. The lower two panels include the results derived from the GCM ECHAM6 (not used to drive any regional model in this study) and the two reanalysis products.
In the CCLM simulation (1655-1999 AD) the relationship between P50 wind speed and mTemp (Fig. 3) in general shows a negative correlation or correlations close to zero, with the exception of a land area east of the Baltic where the correlation is positive. Hence, in this regional simulation, colder periods can be related to generally stronger winds. In contrast, the correlation between the mean temperature and P50 in the ECHAM5 simulation is positive in the northern part of the domain, with some regions showing negative correlations in the south, in-between correlations are around zero. The MM5 and the ECHO-G simulations also display a qualitatively similar, but even more clear, contrast. In the MM5 simulation, median wind speeds over land areas are clearly negatively correlated with mean temperature, whereas over oceanic areas P50 is positively correlated with mTemp. The ECHO-G simulation displays clearly positive correlations over the whole area, similar to the ECHAM6 simulation and the regional reanalysis products.
Concerning the link between mTemp and the width of the distribution (SD) the regional simulations and ECHAM5 the SD correlates negatively with the mean temperature (Fig. 4a, b and c), indicating that colder periods tend to be related with a wider distribution of wind speed over the whole area. This is also supported by the negative correlation between mTemp and the difference P99-P50 (not shown, Table 3). Since P50 in the regional models was also negatively correlated with mTemp, the regional models tend to be associated with broader wind speed distributions and with stronger mean winds.
The ECHO-G simulation (used to drive MM5) displays, in contrast, positive correlations in a region along the North and South Baltic Sea, straddled by regions of zero or negative correlations in Scandinavia and central Europe (Fig. 4d).
We assume that the results showing colder periods correlating with stronger winds may be induced by a stronger meridional temperature gradient (tGrad; see Sect. 3). This assumption is strengthened by the negative correlation between mTemp and tGrad, which shows a spatially averaged value of −0.56 in the CCLM simulation and −0.47 in the MM5 simulation. This negative correlation indicates that lower mean temperatures tend to occur with stronger meridional temperature gradients in the regional simulations, and thus, tGrad could be the primary driver for the changes in the wind speed distribution. Therefore, we also analysed the relationship between the parameters of the wind speed distribution and tGrad (Table 3). This analysis reveals contrasting results between the two regional simulations: the correlations in the CCLM simulation are only minor and not statistically significant, whereas in the case of MM5 they are relatively strong for all parameters of the wind speed distribution, except for the difference P99-P95.
We additionally investigate the relationship between the mean NAO index and the distribution of wind speeds. The correlation coefficient between the NAO index (Sect. 3) and the median wind P50 displays clear similarities of the regional simulations and the driving global simulations, respectively (Fig. 5). However, the differences between the results provided by the regional models and those provided by the global models are profound. The correlation patterns derived from the regional models display positive correlations over oceanic areas, but in general negative correlations over land areas. This is reflected by a west-east correlation dipole (in the CCLM simulations the zero of this dipole does not coincide with the coast line). In contrast, the global simulations display positive correlations over the whole region (Fig. 5b, d, e). The idea that a positive NAO index should be associated with stronger winds in general is only confirmed in the global simulations. The regional models indicate that a stronger NAO tends to be linked to stronger winds only over the ocean and coastal areas, but not over Central and Eastern Europe. This difference may hint to an influence of the boundary layer parametrization in regional models over continental areas.
The spatially averaged correlation between the NAO and SD is very low in the CCLM simulation (0.12), but much stronger and negative in the MM5 simulation (−0.42). This difference can now be explained by the different correlation patterns shown in Fig. 6. The spatially resolved map showing the correlation between the NAO index and the grid-cell SD in the regional models CCLM and MM5 and in the global model ECHO-G is remarkably similar to the correlation patterns between the mean temperature and the SD (compare Figs. 6 and 4). Again, the CCLM simulation shows positive correlations over North-Western Europe whereas the correlations in MM5 are negative over the whole region. The global model ECHO-G tends to show higher positive correlations over the North Sea and the Southern Baltic Sea, straddled by negative correlations over Central and Eastern Europe. An exception is the result for ECHAM5, although the spatially averaged correlation between NAO and mTemp shows a high positive value of 0.57 (Table 3), the patterns for mTemp-SD and NAO-SD show completely different signs. A comparison with a different ECHAM5 simulation (not shown) with a weaker solar forcing (Krivova and Solanki, 2008) showed rather comparable results between mTemp-SD and NAO-SD. Which leads to the conclusion that the external forcing of each simulation plays a crucial role for the relationship of mTemp and the wind speed distribution.
As already known by the scientific community NAO and mTemp are correlated, hence it is not surprising that in most simulations both show comparable relations to the wind speed distribution. Nevertheless, due to internal model variability and different resolutions the global and regional models show different spatial fingerprints for the correlation maps.

The global model ECHAM6/MPIOM
In general the correlation patterns between the large-scale drivers and the parameters of the distribution of wind speed resemble those obtained with ECHAM5 and ECHO-G (ECHAM4/HOPE-G), but some clear differences exist. The higher spatial resolution of ECHAM6 does not, however, lead to correlation patterns that resemble those derived from the regional model simulations CCLM and MM5, pointing towards changes in the physical parametrization (i.e. PBL -Planetary Boundary Layer -scheme) as the main factor explaining the differences in the simulations. The correlation patterns between the median wind speed P50 and the mean temperature or the NAO index in ECHAM6 are indeed similar to the ones derived from ECHAM5 and ECHO-G, displaying generally positive, albeit weak, correlations between the median winds and temperature (Fig. 3). The correlations between the median wind and the NAO index are positive and strong (Fig. 5). However, the correlations of these two driving factors with the width of the wind speed distribution, represented by SD, differ in ECHAM6 from the other two versions of ECHAM ( Fig. 4 and Fig. 6). ECHAM6 displays correlation patterns that are positive and spatially more homogeneous, whereas the ECHAM5 and ECHO-G simulations show higher correlations in the Southern Baltic surrounded by negative correlations in Scandinavia and Central Europe. Again, an exception is ECHAM5 correlation between mTemp and SD, which shows an over all negative relationship.
Physically, the relationship between mTemp and median wind is positive (but statistically not significant) in the ECHAM6 simulation. This indicates that warmer temperatures are accompanied by a shift towards higher wind speeds and by a slight tendency to a broader wind speed distribution (see also Fig. 4e). tGrad shows negative correlation with the median wind, consistent with the link between a decreased temperature gradient in warmer periods.

The reanalysis data coastDat2 and NCEP
We present the link between the large-scale drivers and the wind speed distribution for the two reanalysis products NCEP and coastDat2. It can be argued that these two reanalysis data sets should be closer to the real climate, because they both incorporate information based on meteorological observations. On the other hand, the reanalysis models are integrated over a relatively short period of time of about 60 years. Therefore the decadal-scale links between the large-scale climate drivers and the probability distribution of wind speed derived from these data sets is most likely afflicted with a higher degree of uncertainty. The correlation patterns derived from NCEP and coastDat2 are based on gliding 5-year windows, instead of 30-year windows as for the longer simulations described before. In addition, this period has witnessed external climate forcings that are quite different from the natural forcings of the previous centuries. Nevertheless, the links between temperature and wind speed distribution should be independent of the radiative forcings that drive the surface temperatures.
The correlation between mTemp and the parameters of the wind speed distribution for coastDat2 and NCEP shows generally significant positive values with P50, SD, P95 and P99, but no significant correlation for diffM and diffE. The spatially resolved correlation between mTemp and P50 (Fig. 3f, g) and the SD field (Fig. 4f, g) is also dominated by positive values, with highest coefficients over southern Norway, and weaker or slightly negative correlations in the southern regions of the domain. Therefore, the correlation patterns in the reanalysis data represent a shift of the wind speed distribution from low to high wind speeds during warmer periods, with a tendency to have a wider wind speed distribution in the northern regions and a small influence of temperature in the southern regions.
The relationship between tGrad and mTemp is negative (−0.25 for coastDat2, −0.52 for NCEP), again showing that colder periods are related to stronger meridional temperature gradients in agreement with all other models analysed here. Thus, the results obtained from the reanalysis products resemble more closely the ones derived from the global climate model simulations (ECHAM5, ECHAM6, ECHO-G) than from the regional simulations (MM5 and CCLM).
The correlation between tGrad and the distribution of wind speed is found to be predominantly weak in the coastDat2 data, with the only mentionable value (0.34) for the correlation tGrad-diffM. This result means that higher temperature differences between North and South are slightly correlated with a broader wind speed distribution. In contrast, for the NCEP reanalysis, the link is strong but opposite: weaker meridional gradients are linked to stronger median winds and wider wind speed distributions. Regarding the link to the NAO, both reanalysis data sets display a consistent picture, with a positive NAO closely linked to stronger median winds and wider distributions (SD) in most of the domain. This link is stronger over the northern regions and becomes smaller and even negative over the southern fringes of the domain. Again, this spatial structure resembles more closely the structure provided by the global models and differs from the pattern provided by the regional models.
4.5 Results in the overlapping time periods 1655-1990 and 1948-1990 This section is dedicated to the comparison of the above explained results with results for the overlapping time periods without (1655-1990 AD, 30-year running mean) and with reanalysis data (1948( -1990. Regarding the overlapping period 1655 to 1990 the main conclusions remain the same for both, Table 3 and Figs. 3-6. In Table 3 some correlations become higher (mTemp), whereas some stay at the same level (NAO, tGrad -beside ECHAM6 which now shows values around zero) and some are lower (mTemp-tGrad). The conclusions concerning the spatial correlation maps are also very comparable: due to almost the same time period all four results for CCLM in the overlap show almost identical patterns. The GCM results are also very similar between both periods with slightly higher values for the overlapping period. Only MM5 shows a difference above Western and Central Europe where positive values oc-cur for the overlapping period (1655-1990 AD) and negative or values around zero for the whole period (1001( -1990. Regarding the overlapping period 1948 to 1990 the values and patterns change. Nevertheless, each model simulation still shows different results, and they do not become more similar to the reanalysis data (Figs. 7 and 8). Figure 7 shows the correlation pattern of mTemp and P50 for the period 1948-1990 which shows more negative (positive) areas for MM5, ECHO-G, ECHAM6, coastDat2 (CCLM, ECHAM5) compared to the results of the whole available time periods.
The results for NCEP change only marginally. For the correlation between mTemp and SD (Fig. 8) we also see more negative (positive) areas for MM5, ECHO-G, coastDat2, and NCEP (ECHAM5). CCLM shows a shift from a positive region over the Benelux area to Scandinavia. ECHAM6 stays comparable in both periods. Note that for all data sets these results are less robust due to the fewer analysed values.

Centennial-scale evolution of the wind speed variance over the past millennium
In the previous sections we analysed the links between largescale atmospheric drivers and the distribution of wind speed at decadal and multidecadal timescales. The time series of the width of the wind speed distribution over the past millennium indicate, however, that the slowly changing soil boundary conditions may also have a strong influence on the longterm evolution of the variability of wind speed in Northern Europe. Figure 9a shows the time series of the spatially averaged standard deviation of the wind speed distribution at each model grid-cell for the simulation conducted with the model ECHAM6. The most remarkable feature of the averaged SD is its almost continuous increase during the simulated period. This monotonous increase is also seen in the corresponding time series calculated with the output of the ECHAM5 simulation (not shown), but not in the data of the ECHO-G (based on ECHAM4) simulation (Fig. 9b). A sug-gestion about the origin of the increase in the standard deviation of the wind speed distribution can be obtained by comparing the spatially resolved SD in the last decades versus the initial decades in the simulation. Figure 10a shows the ratio between the spatially resolved standard deviations for the periods 1871-1990 AD (P1) and 1001-1091 AD (P2), respectively. The values of this ratio are higher than unity (SD larger at the end of the simulation) over the land areas of central Europe, with a maximum at about 25 degrees east. The standard deviation over oceanic grid-cells and over Scandinavia does not change significantly between these two periods in the simulations. Again, this spatial pattern of increase in the width of the wind speed distribution is also simulated by the ECHAM5 simulation (not shown), but not in ECHO-G where the ratio is scattered around 1 (not shown).
The spatial pattern of changes in SD between the beginning and end of the ECHAM6 and ECHAM5 simulation suggests that the increase in the width of the wind speed distribution may be related to surface-boundary processes. This suggestion is supported by the changes in forest cover in the course of the last millennium as reconstructed by Pongratz et al. (2008). This reconstruction was used to drive the models ECHAM6 and ECHAM5 (see Sect. 2). The difference in tree fraction in each model grid-cell between the periods 1871-1990 AD (P1) and 1001-1091 AD (P2) is shown in Fig. 10b. The spatial agreement between the reduction in tree fraction and the widening of the wind speed distribution between the beginning and the end of the millennium is remarkable and strongly supports the hypothesis that the distribution of wind speed is mainly affected by land-use changes and related changes in surface roughness length. This is supported by the analysis of a third time period from 1571 to 1690 AD (P3) which also shows a strong agreement for the SD ratio (P3 / P2; Fig. 10c) and the tree fraction difference (P3-P2; Fig. 10d), albeit with less intense values presumably due to the less intense deforestation in period P3. Hence, we conclude that a less extensive forest cover causes a widening of the wind speed distribution, and vice versa. This is also visible for the time series in Fig. 11a, which shows the temporal evolution of the tree fraction (black line) and the 30-year running mean SD for ECHAM6 (green line in Fig. 9a multiplied by −1), both lines show a remarkable agreement in the long-term evolution. The simulation with ECHAM5 shows a comparable evolution for the SD (not shown). The simulation with the model ECHO-G, which was not driven by changes in land use, does not show a longterm increase or change in the width of the distribution of wind speeds (see green line in Fig. 9b), supporting the strong influence of land cover changes on the distribution of wind speeds.
Therefore, at multi-centennial timescales the correlation between the wind speed distribution and temperature that was explored in the previous sections, for ECHAM5 and ECHAM6, could have been indirectly caused by land-use changes. At these timescales, anthropogenic deforestation and mean temperature exhibit a positive trend. Thus the expansion of the wind speed distribution and the increase of temperature in these decades might be induced by physically different factors, leading to positive correlations in our analysis. Figure 11a shows the temporal evolution of the tree fraction used to drive ECHAM5 and ECHAM6 (black line). This proves that on shorter timescales (e.g. 100 or 200 years) the effect of this surface-boundary process is negligible, especially before the 17th century when the trend is very weak. Figure 11b, c exhibit the correlation between mTemp and the SD for the periods P1 and P2 calculated with a 5-year running mean. Hence, these figures show the relationship between mean temperature and the wind speed distribution independent of the deforestation effect. Therefore, this statistical effect can be disentangled by separating the analysis of these simulations into two time periods (TP1: years 850-1500 AD, TP2: years 1500-2005 AD).
The correlations between mean temperature and the width of the wind distribution does show a difference in the correlation. For TP1 of ECHAM6 mTemp-diffM is around 0 and mTemp-diffE −0.25, for P2 0.63 and 0.57, respectively. For TP1 of ECHAM5 mTemp-diffM is −0.26 and mTemp-diffE is −0.38, for TP2 0.26 and 0.12, respectively. Both time periods in the GCMs show different signs for the re-lation between mean temperature and the shape of the wind speed distribution. These results suggest that specifically for ECHAM6 the correlation between the mean temperature and the width of the wind speed distribution is a statistical artefact mediated by deforestation.

Discussion and conclusions
This study investigates and compares different simulation data sets and reanalysis products, on timescales covering the last decades to the past millennium, regarding the probability distribution of the daily wind speed in winter time over Northern Europe. Our investigation is aimed at identifying the large-scale factors that drive changes in this probability distribution. The study is based on correlations between different parameters of the wind speed distribution and different climatic indices related to mean temperature, meridional temperature gradient, and the North Atlantic Oscillation. The overlaying question is whether and how the wind speed distribution may change during varying climate conditions and hence whether these conditions may provoke more and/or stronger wind speed extremes.
One prominent result is that the link between the thermal indices and the North Atlantic Oscillation appears physically consistent in all data sets, and thus all models are consis- Figure 11. (a) tree fraction after Pongratz et al. (2008) (black) averaged over the investigation area. ECHAM6 30-year running mean SD of wind speed (multiplied by −1) (blue). (b) Correlation between field mean temperature and SD for 1871-1990 with a 5-year running mean computation (c) correlation between field mean temperature and SD for 1001-1091 with a 5-year running mean computation. tent in this regard. The relationship between the NAO and mean temperature over Europe is a well-known effect (Rutgersson et al., 2015), in wintertime this is also positive in all models. This effect is caused by the advection of maritime air masses by stronger westerly winds in a more positive NAO state (also discussed in Gómez-Navarro et al., 2015a). Also, the correlation between mean Temperature (mTemp) and the temperature gradient (tGrad) shows negative values for all models, indicating that warmer periods are linked to a weakened meridional temperature gradient. This might be explained by the fact that northern regions warm (cool) more strongly when the overall temperature is higher (lower). In climate change projections this is referred to as polar amplification, although it has been also identified in paleoclimate simulations and reconstructions over the past millennium (Luterbacher et al., 2002).
A second important result is that the correlation between the large-scale indices and the parameters of the wind speed distribution exhibit markedly different results among the data sets analysed and it is difficult to derive general conclusions on the effect of these large-scale drivers on the distribution of daily wind. Comparable difficulties are reported by Fischer-Bruns et al. (2005). All in all the expected link between the mean temperature and the wind speed distribution and between the mean temperature gradient and the wind speed distribution Li and Woollings (2014) is not confirmed in each simulation. In our study the three global simulations present similarities that warrant to place them in one group. This can be expected to a certain degree because the three atmospheric models included in the GCM belong to the same ECHAM family. Further analyses with completely different GCM families would be necessary to reinforce our findings. Likewise, the correlation patterns obtained from the two regional models are also generally similar, although in this case both models, MM5 and CCLM, are structurally different.
The striking result is that the regional models do not seem to inherit the dynamical properties of their respective global models, but produce instead different correlation patterns between the large-scale drivers and the wind speed distributions. It is plausible that the higher resolution and the different parametrization schemes of the boundary layer shape the link between the large-scale dynamics and turbulent processes that modulate the width of the daily wind distribution. As Hall (2014) and Gómez-Navarro et al. (2013) already stated, the RCMs should provide a better representation of small-scale processes, topographic influences and of the land-sea contrasts, and thus should be better suited for the simulation of extreme events.
Another indicator for the influence of the spatial resolution on our results might be the fact that only the regional simulations MM5 and CCLM show strong negative correlations between mean temperature and the width of the probability distribution as measured by diffM/diffE. These correlations suggest that colder periods are connected with stronger wind speed extremes. In contrast, the GCM data show no clear correlations between these parameters. This study does not allow us to provide a comprehensive dynamical explanation for the different behaviour of wind speeds in changing temperature or pressure conditions: the models show different results although each model seems to be dynamical consistent in itself. Therefore, a detailed analysis of each of the simulations, and maybe of the physical parameterization and computer codes, becomes necessary to understand how the different correlation patterns arise.
On centennial timescales, we identified land-use changes as a very important factor modulating near-surface wind in the simulations. Note that anthropogenic changes in land-use are prescribed only in the ECHAM5 and ECHAM6 simulations, whereas for ECHO-G land-use is kept constant during the whole simulation. The analysis of the ECHAM5 and ECHAM6 millennium simulations reveals a strong increase of the standard deviation of wind speed for the last decades since the industrialization, and in areas that coincide with larger deforestation along the last centuries. The impact of land-use changes on wind conditions was also shown by Pessacg and Solman (2013) in simulations with the regional model MM5 over South America for different idealized landuse scenarios.
The main conclusion that can be drawn from this study is that the link between large-scale climate drivers and the distribution of daily wind speeds in wintertime in this region is complex and not fully constrained by currently available simulations. All models analysed here have been individually profusely used in climate simulations and the data sets have been used in a number of other previous studies, and no gross deficiencies have been pointed out so far. We conclude that, although climate models may be dynamically sound in the large-scale contest, the impact of climate change on variables like near-surface wind speed distribution possibly depends more strongly on the details of the physical parametrization and changes in surface forcing, like deforestation, than on the large-scale dynamical drivers, such as large-scale temperature or sea-level-pressure changes.