Near-surface mean and gust wind speeds in ERA5 across Sweden: towards an improved gust parametrization

The ERA5 reanalysis product has been compared with hourly near-surface wind speed and gust observations across Sweden for 2013–2017. ERA5 shows closer agreement than the previous ERA-Interim reanalysis with regard to both mean wind speed and gust measurements, although significant discrepancies are still found for inland and mountainous regions. Therefore, attempts have been made to improve formulations of the gust parametrization used in ERA5 by adding an elevation-dependency and by adjusting the convective gust contribution. Major improvements, especially over mountain regions, are achieved when the elevation differences among the stations are considered. Closer agreement between the observed and parametrized gusts is reached when the convective gust contribution is also tuned. The newly designed gust parametrization was also tested for Norway, which is characterized by more complex topography. Wind gusts from the selected Norwegian stations are more realistically simulated when both the elevation-dependency and the tuned convective contribution are implemented, although the parametrized gusts are still negatively biased. Such biases are not explained by the different in gust duration in recorded wind gusts between Sweden and Norway.


Introduction
Globally, storms with heavy rainfall and strong winds have been identified as by far the costliest among various types of climate-related and geophysical disasters, with twice the reported losses for either flooding or earthquakes (Wallemacq et al. 2018). Regionally, windstorms and extreme winds contribute to more than half of the economic losses associated with natural disasters in Europe (Ulbrich et al. 2013). Across Scandinavia, the occurrence of extreme winds can affect aviation security, as well as damage buildings and forests, representing severe hazards to people, properties and transport (Achberger et al. 2006;Suomi et al. 2014). For example, in Sweden, where forests cover 56% of the land, with 95% of those forests being used for the timber industry, wind damage could seriously affect the national economy (Hannon Bradshaw 2017).
Driven by the twenty-century surface air temperature rise, extreme wind events could change in their frequency and magnitude of occurrence (Pryor et al. 2010;Azorin-Molina et al. 2016;Wu et al. 2018;Safaei Pirooz et al. 2019). To understand past's extreme wind conditions in a warming climate, multi-decadal wind observations are needed. To define the occurrence of extreme wind events, just looking at the near-surface (~ 10 m height) mean wind speed (hereafter, WS), i.e. the mean wind over the last 10 min (WMO 1987), is not sufficient. To capture the abrupt increase in wind and its turbulent signature, the World Meteorological Organization (WMO 1987) suggests also recording the so-called near-surface (~ 10 m height) peak or wind gust (hereafter, WG), defined as the maximum 3 s wind speed. By definition, 1 3 WG can capture the turbulence fluctuations due to the short averaging time of the WS calculation and can provide complementary information to WS climatology, particularly for determining the occurrence of severe wind events.
Unfortunately, the quality of extreme wind measurements in Sweden is poor due to the lack of homogeneous and continuous gust observations. Only with the introduction of automatic weather stations using Thies ultrasonic anemometers in 1996 (Minola et al. 2016;Lekander 2019, personal communication), have consistent and continuous wind gust measurements become available across the country. Therefore, alternative datasets should be explored to properly investigate long-term WG changes across Sweden. Reanalyses, by combining past observations with models, are able to generate consistent series of different climate variables and nowadays are among the most-used datasets in geophysical research (Dee et al. 2011;Hersbach et al. 2018). Before global reanalysis products are used for assessing changes and impact of extreme winds, however, their ability in representing observed near-surface wind statistics must be demonstrated. Different reanalysis products have already been explored for modelling wind power across Sweden (Liléo and Petrik 2011;Olauson 2018), but no studies have tested if reanalyses like ERA-Interim and the new ERA5 are able to realistically simulate near-surface wind gusts over the country.
In models used in current reanalyses, eddies of all scales responsible for gusts cannot be resolved explicitly and WG is thus parametrized. In fact, numerical models are not able to resolve explicitly those processes with scales of motion smaller than those resolvable by the models (i.e., the gridscale). For this reason, it is necessary to relate the cumulative effects of these subgrid-scale processes to the resolvable scales of motion using the given scale variables: this approach is named parametrization (Anthes 1985). In the reanalysis products of the European Center for Medium-Range Weather Forecasts (ECMWF), wind gusts are parametrized and calculated by post-processing (ECMWF 2016). In fact, the parameter gust is computed at every time step and its maximum since the previous post-processing time is written out and archived as WG. The parameter gust is calculated at each time step of the model as the sum of three terms: (1) the instantaneous 10 m wind speed; (2) the turbulent gustiness (i.e., the standard deviation of the horizontal wind; hereafter ∆WG turbulent ); and (3) the convective contribution (hereafter ΔWG convective ). It was Panofsky et al. (1977) who proposed, on the basis of the similarity theory, an empirical relationship to express ΔWG turbulent in terms of the shear-stress (expressed in term of the friction velocity) and the stability of the boundary layer (expressed in term of the Monin-Obukhov length). This hypothesis has been tested over uniform surfaces. The convective term in the sum was added later by Bechtold and Bidlot (2009) to include the gustiness in deep convective situations with strong wind shear, as frontal systems and long-lived organized mesoscale convective systems.
This study aims to compare the new ERA5 reanalysis product for hourly observations of WS and WG across Sweden. In particular, ERA5 wind simulations are compared against the previous ERA-Interim reanalysis dataset. Furthermore, this study examines the gust parametrization adopted in ERA5 to understand if the estimated gusts with the current formulation match the measured gust. An improved formulation for the parametrization of gusts is established across Sweden and is also tested for Norway to check if such improved parametrization could be applied successfully to other regions.

Wind observations across Sweden
In this study hourly wind observations (WS and WG) across Sweden were provided by the Swedish Meteorological and Hydrological Institute (SMHI). All wind series are available online at the SMHI open data page (https ://www.smhi.se/ data/utfor skare n-oppna -data/; last accessed 8 May 2020). Observed wind series have been selected for the 5-year period 2013-2017 because during this time almost all measurements are labelled as "controlled and approved" by SMHI (as specified in the metadata information of the downloaded series). Before 2013, measurements are all marked as "suspected", i.e. their quality has only been roughly controlled, or even not controlled at all, by SMHI. The few observations marked as "suspected" during 2013-2017 have been discarded for this study. Therefore, no homogenization is needed as the quality of the recorded WS and WG values was ensured by the SMHI. Wind series from 90 weather stations have been chosen as those stations: (1) are fairly evenly distributed across the whole of Sweden (Fig. 1); (2) have both hourly WS and WG measurements for 2013-2017, when the high-quality of hourly measurements is guaranteed by SMHI; and (3) have no more than 6% (around 100 days or 3.5 months) of missing or "suspected" hourly observations. Following the WMO's guidelines (WMO 2014), WS were measured at 10 m height as the average wind over a 10 min interval for each hour (Lekander 2019, personal communication); while WG refers to the maximum 2 s gust recorded in the last hour at a 10 m above the ground. This deviates from the standard 3 s averaging time proposed by WMO (1987), and is shorter by 1 s. Automatic weather stations across Sweden measure wind using Thies ultrasonic anemometer 2D (Wern and Bärring 2009; sensor specifications can be found at https ://www.thies clima .com/en/Produ cts/Wind-Ultra sonic -Anemo meter /?art=809; last accessed 8 May 2020).

Wind observations across Norway
To test the newly formulated gust parametrization for a region with more complex topography than Sweden, hourly WS and WG observations across Norway have been selected. Wind series were provided by the Norwegian Meteorological Institute (MET Norway) through the eKlima web portal (https ://shark i.oslo.dnmi.no/porta l/page?_pagei d=73,39035 ,73_39049 &_dad=porta l&_schem a=PORTA L; last accessed 8 May 2020). Wind series across Norway have been chosen to cover the 3-year period 2015-2017 because the quality of the recorded wind measurements is ensured to be "OK" by MET Norway (as reported in the metadata info of the downloaded series), with only a few observations labelled as "uncertain/erroneous" (which are excluded here). For this study, 15 weather stations have been selected as they: (1) span the 3-year period 2015-2017; (2) have measured both hourly WS and WG; (3) cover both coastal and inland regions (Fig. 1); and (4) have no more than 5% (around 55 days or 2 months) of missing or "uncertain/erroneous" hourly observations. Similar to the Swedish weather stations, WS across Norway is calculated as the average wind over the last 10 min of each hour; but differently, WG is the maximum 3 s wind speed in the last hour, following the output-averaged time recommended by WMO (1987). MET Norway uses both Gill Wind Observer II and Young Wind Monitor-MA anemometers to measure near-surface wind; a few stations are still equipped with Vaisala WAV151 (Husebye 2019, personal communication).

ERA5 and ERAINT reanalyses
Reanalysis data provide spatially complete and physically coherent simulated data of various global climate variables by using a forecast model in which information from observations of various types and from multiple sources are assimilated (Dee et al. 2011). Physical coherence is ensured as estimated variables are made to be consistent with the law of physics of the forecast model as well as with the observations; spatially completeness is achieved by using the model equations to extrapolate information from observed variables to unobserved parameters at nearby locations. In this study, two reanalysis datasets from ECMWF are tested for wind observations: ERA5 and ERA-Interim (hereafter, ERAINT). ERA5 is the latest climate reanalysis produced by ECMWF, which has replaced the highly successful ERAINT (started in 2006). ERA5 is based on the Integrated Forecasting System (IFS) Cycle 41r2, in operation in 2016, while ERAINT Fig. 1 Elevation map of the study area with the location of the 15 weather stations across Norway (red circles) and the 90 stations across Sweden classified as coast (blue-green circles), inland (yellow circles), and mountain (violet circles) stations uses the 2006 forecasting model IFS Cycle 31r2. Therefore, ERA5 benefits from 10-year of developments in model physics and data assimilation (Hersbach et al. 2018). The new ERA5 includes hourly outputs (versus 3-hourly ERAINT outputs), with an increased horizontal resolution of 31 km, compared to the 80 km of ERAINT. Besides a considerable increase in horizontal and vertical resolution and a decade of improvements in representation of model processes and data assimilation, ERA5 has improved in the ingested number of observations, with newly available observed datasets assimilated. In fact, ERAINT has become progressively outdated in representing physical processes and the number of observations which it can take has declined. In ERA products like ERAINT and ERA5 all near-surface wind observations over land are excluded as they cannot be usefully interpreted by the data assimilation system (Dee et al. 2011). Instead, across land, vertical wind speed profiles from multiple sources (radio-and aircraft-sondes, satellites, etc.) are used as inputs (https ://confl uence .ecmwf .int/displ ay/ CKB/ERA5+data+docum entat ion#ERA5d atado cumen tatio n-Obser vatio ns; last accessed 8 May 2020).
For both ERA5 and ERAINT 10-m height zonal (u) and meridional (v) wind components are outputted and WS at each grid point is calculated as the square root of the sum of the square of the wind components: WS = √ u 2 + v 2 . The parameter gust is computed at every time step using the parametrization of Panofsky et al. (1977) and Bechtold and Bidlot (2009) (see Sects. 1 and 4.2) and its maximum during the last post-processing period (3-h for ERAINT; 1 h for ERA5) calculated as WG (ECMWF 2016). ERA5 also provides the instantaneous gust (hereafter, ERA5 instantaneous or ERA5_inst), i.e. the parameter gust calculated for a given time (each hour).
For this work, to investigate similarities and discrepancies between wind observations and outputs from the selected ERA reanalyses, the observed wind series from a given measuring station is compared with the wind series from the closest reanalysis grid point, assuming that this closest series matches the observed one better than any other series from more distant grid points.

Spatial scale definition
The ability of gridded datasets of 80 (ERAINT) and 31 km (ERA5) in representing the wind spatial variability of the study region is investigated by defining the WS and the WG spatial scale. Following Chen et al. (2016), the Pearson's correlation coefficients between the WS (WG) series of each station and all the others are calculated and plotted against their distance. The e-folding distance of the correlation decay is thus defined as the WS (WG) spatial scale and it is quantified by fitting an exponential decay function (hereafter, ζ for WS and ψ for WG) over the scatter between correlation and spatial distance. In fitting the exponential decay functions ζ and ψ, we force them to be 1 for distance equal to 0 (highest possible value for autocorrelation) and to never be negative (by assuming that it is not possible to have any inverse correlation: the temporal correlation of two series which strongly differ in their temporal variability tends to be equal to 0). Therefore, the fitted exponential function is calculated by using only those points in the scatterplot which give a decay function positive for any spatial distance.
The value assumed by ζ (ψ) for the horizontal spatial resolution of a given gridded datasets tell us the temporal correlation that observed WS (WG) series have on average at the dataset spatial scale, which can or cannot justify the comparison of station data with the investigated gridded dataset.

Model output statistics
The gust parametrization currently used by ECMWF is improved in this study by using the so-called Model Output Statistics (MOS) technique. Defined by Glahn and Lowry (1972), MOS refers to the statistical post-processing which is used to improve the numerical model's ability to forecast local weather conditions by relating model outputs to local observations. In this study, by processing the wind information from the available weather stations across Sweden, a new formulation of the gust parametrization is defined; subsequently this formulation is evaluated using stations across Norway to further test the validity of the new parametrization.

Statistics for comparison
To evaluate the match/mismatch between the observed wind series and the estimated ones from different reanalysis products or from different gust parametrizations, the following statistical tests are used: (1) Pearson's correlation coefficient, (2) Root Mean Square Error (RMSE), and (3) bias. These statistics are calculated for wind series resolved at different time scales: (1) hourly or 3-hourly series; (2) daily series, calculated as daily wind means; (3) daily maxima (hereafter daily max) series, i.e. the sequence of maximum wind values recorded each day; and (4) monthly series, calculated as monthly wind averages.
To measure the degree of association (or linear relationship) between two wind series X and Y (i.e., observed and estimated wind series), the Pearson's correlation coefficient is calculated as in Eq. (1) (Gibbons and Chakrborti 2003): where cov is the covariance between X and Y, and var stands for the variance. Before testing for the existence of a correlation between X and Y, it is necessary to remove traces of periodicity which may lead to unreasonably large ρ values (Marzban and Schaefer 2001). For this reason, mean diurnal and annual cycle are removed from hourly series before ρ is calculated; in monthly series, mean seasonal cycle is cleared away. For daily and daymax series, the mean annual cycle is removed prior to the ρ computation.
A good wind estimation must produce simulated wind in the proximity of the true (observed) wind realization (Von Storch and Zwiers 1999). To mathematically express the vicinity of such realizations, this study uses the RMSE for defining such "distance". The RMSE of the estimated wind ŵ t for times t, with observed wind w t over T times, is calculated as in Eq. (2): The RMSE allows different wind estimations to be compared. The estimator (reanalysis product or parametrization) that has RMSE less than all other estimators is considered the best. However, other properties are desirable, such as the lack of bias. The bias of the estimated wind ŵ t for times t, with observed wind w t over T times, is computed using Eq. (3): Bias tells us if there is a tendency to deviate from the true wind realization.

Classification of measuring stations
As wind climatology is strongly dependent on terrain characteristics of the region (Jimenéz and Dudhia 2012), measuring stations across Sweden are classified according to the terrain characteristics around the stations. By plotting mean WG or WS conditions for each station during 2013-2017 versus the station distance to the sea and its elevation (Fig. 2), measuring stations across Sweden have then been classified into three groups: (1) coast stations, if located no further than 20 km from the sea; (2) inland stations, if located further than 20 km from the sea and with an elevation lower than 750 m above sea level (a.s.l.); and (3) mountain stations, if located further than 20 km from the sea and with an elevation greater than 750 m a.s.l.. This classification was suggested to be able to group weather stations influenced by similar physical processes (Achberger et al. 2006): (1) coast stations are dominated by regional and local circulations due to the difference between sea and land (e.g., sea breezes, turbulence due to step-like change in surface roughness; Gustavsson et al. 1995;Borne et al. 1998); (2) inland stations are strongly influenced by land surface processes and land-atmosphere interactions; and (3)

ERAINT vs ERA5
ERAINT and ERA5 wind outputs are compared with observed WS and WG series to evaluate the skills of these reanalysis products in simulating wind across Sweden. The reliability of the comparison of station data with gridded datasets of 80 (ERAINT) and 31 (ERA5) km of horizontal spatial resolution is justified by the spatial scale analysis (see Sect. 3.1 for the spatial scale definition). As shown by Fig. 3, both the fitted decay functions ζ and ψ have Pearson's correlation ~ 0.8 and ~ 0.9 for the ERAINT (80 km) and the ERA5 (31 km) spatial resolution, respectively. These high correlation values support the use of ERAINT and ERA5 to explore the spatial variability of wind across Sweden because, with their horizontal resolution, they are able to reasonably capture the spatial coherence of the observed datasets.
To evaluate the performance of ERA5 in simulating wind conditions across Sweden, the seasonal (Fig. 4) and the diurnal (Fig. 5) cycles (defined as the monthly and hourly mean, respectively) of both WS and WG are plotted for observations, ERAINT and ERA5 for 2013-2017. Figure 4 shows the mean WS and WG seasonal cycle for the various datasets across the three regions. In the coast stations, both ERA products reproduce well the observed seasonal cycles of WS and WG well, with the strongest winds in winter (e.g., December) and the weakest winds in summer (e.g., July). However, both reanalyses show a weak overestimation of up to 1.5 m s −1 , particularly ERAINT. Across inland stations, the mean seasonal cycle has a less season-dependent signal and the two reanalyses also overestimate wind conditions by up to 2 m s −1 for WS and 4 m s −1 for WG, with ERA5 always being closer to the observed seasonal cycle. For WG, differences between ERA5 and ERA5_inst are small (~ 0.2 m s −1 ): ERA5_inst always shows slightly lower values than ERA5, with fewer deviations when compared to the observations. This is because, by definition, ERA5_inst only selects one gust parameter, excluding all the other possibly greater calculated gusts during the previous hour, which, if higher, would lead to a greater WG output in ERA5. For the mountain stations, observed WS and WG also show a seasonal cycle as displayed for the coast ones, with the strongest and weakest winds during the cold (i.e., December) and warm (i.e., July) months, respectively. Over mountain areas, all reanalysis products show the largest underestimation among the regions, particularly for WS from ERA5.
The mean diurnal cycle for each region is plotted in Fig. 5. Observed and reanalyzed mean daily cycles of WS and WG show higher WS and WG during the day compared to the night-time hours, with a maximum peak around 15 CEST (Central European Summer Time; local time for Sweden) and minimum peak around 4 CEST. For the coastal stations this cycle is flattest, with less pronounced differences between day and night. For the costal stations, less pronounced differences between day and night are found, particularly for WS. In coast stations, both ERAINT and ERA5 (and ERA5_inst) slightly overestimate the observed WS and WG diurnal cycles. For the inland stations, this overestimation becomes even higher throughout the day. ERA5 resembles the observed daily cycle better than ERAINT, with ERA5_inst always performing the best for the same reason mentioned above for the seasonal cycle. In contrast to the coast and inland regions, WS and WG diurnal cycles are strongly underestimated (~ 3 m s −1 ) for the mountain stations for all reanalyses and most of the day, with the ERAINT diurnal cycle less negatively biased than ERA5. Observed diurnal variations in mean WG across mountain areas are also underestimated by the two ERA products, particularly for ERA5 and ERA5 inst, with the largest biases occurring during night-time hours. In general, simulated winds in all regions show the lowest biases occurring during the day compared to night, as daytime heating (especially over land surfaces) increases the turbulent state of the boundary layer and the downward momentum mixing (Borne et al. 1998;Achberger et al. 2006). It is noteworthy to mention that the WS diurnal cycle in ERA5 displays a drop between 11 and 12 CEST, which is not shown for WG. This behavior has been observed for other locations (Hersbach et al. 2018), but the cause is still unknown. An irregular behavior between Scatter-points are clustered into three groups: (1) coast stations (bluegreen), (2) inland stations (yellow), and (3) mountain stations (violet) analysis cycles is probably the reason behind this undesirable feature because the discontinuity occurs exactly between two 12-h analysis windows and no jumps are observed for the forecasts connecting the analyses (Hersbach et al. 2018).
Statistics for both WS (Fig. 6) and WG (Fig. 7) show the best agreement with observations for ERA5. In particular, the mean Pearson's correlation coefficient is higher for ERA5 across all regions and the mean RMSE and mean bias is smaller compared to ERAINT. Similar results are found for ERA5 and ERA5_inst in comparison with observed WG.
Furthermore, it is noteworthy to notice that biases are strongly dependent on WS and WG intensity. As Fig. 8 shows, the high observed WS and WG intervals are the ones which display the strongest negative biases in both ERAINT and ERA5; while larger positive biases occur under weak WS or WG intervals. Therefore, stations with the highest mean WS and WG (Fig. 9), i.e. coast and mountain stations (see Fig. 2), are negative biased. Stations with weak wind conditions on average, as mostly inland stations, show large positive biases.
To summarize, in all statistics ERA5 shows better performance than ERAINT in representing both WS and WG across Sweden, most likely due to the improvements made in the representation of the model processes and data assimilated.
For WG there are no significant differences between ERA5 and ERA5_inst, and the identified discrepancies are due to how WG is calculated for ERA5_inst.

Decomposition of gust estimates by ERA5
To test if the parametrization used by ERA reanalyses can realistically represent the physical processes behind the origin of gusts across Sweden, a close look at the gust parametrization is needed. As mentioned in Sect. 1, the parameter gust (hereafter, GUST) is calculated at each model integration step and its maximum during the latest hour (or during 3-h for ERAINT) is archived as WG (ECMWF 2016). In particular, the gust parameter is the sum of the instantaneous WS (WS instantaneous ) and the total gustiness (ΔWG gustiness ): The turbulent gustiness (ΔWG turbulent ) and the convective (ΔWG convective ) terms contribute to the total gustiness ΔWG gustiness : with (4) GUST = WS instantaneous + ΔWG gustiness .
(5) GUST = WS instantaneous + ΔWG turbulent + ΔWG convective Fig. 3 Definition of spatial scale through the scatterplot of Pearson's correlation coefficient between the all 90 WS (left) and WG (right) series across Sweden for 2013-2017 plotted against the station distance. The red line is the fitted exponential function which defines the typical decays of the spatial correlation of WS and WG series with increasing distance. Circles marked with a black cross are the points that were used for the fitting of the exponential function ΔWG turbulent expresses the gustiness driven by the shearstress (function of the friction velocity) and the stability of the boundary layer (function of the Monin-Obukhov length) and is derived by the field experiments of Panofsky et al. (1977) over uniform surfaces and based on the similarity theory (Arya 2001). It is calculated as: with where C turbulent = 7.71, u * is the surface friction velocity, L is the Monin-Obukhov length, and z i is the height (set up as z i = 1000 m).
Following Bechtold and Bidlot (2009), the convective gustiness contribution (ΔWG convective ) represents gusts generated in deep convective situations (which would have not (6) ΔWG gustiness = ΔWG turbulent + ΔWG convective . (7) , been captured by ΔWG turbulent ) and is calculated as proportional to the non-negative low-level wind shear ΔWS + shear : with where C convective = 0.6 is a tunable convective "mixing" parameter, and WS 850 WS 950 is the difference between the 850 hPa and 950 hPa wind speeds, representing the lowlevel wind shear. In ERA5, instead of pressure level data, wind speeds from model levels 114 and 124 are respectively used for 850 hPa and 950 hPa wind speeds (Hewson 2019, personal communication). Following this formulation, for each station and at each hour during 2013-2017, the parametrized WG is calculated using the observed WS instead of WS instantaneous and using ERA5 outputs (e.g., friction velocity) of the closest grid point to calculate ΔWG turbulent and ΔWG convective . In this way, the parametrized WG does not suffer from biases or errors due to the modelled WS and the contribution of each gustiness term (ΔWG turbulent and ΔWG convective ) can be addressed.  (black), ERAINT (violet), and ERA5 (orange) series; for WG, the mean monthly cycle is also calculated using ERA5_inst (yellow) Figure 10 shows the contribution of the different WG parametrization terms to the total mean seasonal cycle of WG for the different regions during 2013-2017. For the coast stations, the major contribution to the parametrized WG comes from the WS. By adding next the ΔWG turbulent and ΔWG convective inputs, the mean seasonal cycle of the parametrized WG matches the observed cycle well, although it is underestimated slightly. Across inland stations, similar contributions from WS and ΔWG turbulent are shown. Although there is a close match between the observed and parametrized WG during winter months, differences are found during the rest of the year. The observed WG shows seasonal variability, being strongest especially during March-June, whereas the parametrized WG does not change much during the year on average. For the mountain stations, WS contributes mostly to the parametrized WG. With the added ΔWG turbulent and ΔWG convective terms, the parametrized seasonal cycle, marked by stronger WG in summer than winter (with a peak observed during June), resembles the observed one, even thought is constantly underestimated up to 2 m s −1 in each month. Notice how the ΔWG convective contribution in all the three regions does not follow the seasonality of occurrence of deep-convection situations and mesoscale convective systems at the origin of convective gusts in a country like Sweden. For example, their frequency of occurrence is in fact much higher during the warm months (April-September) than the cold months (October-March) across Finland (Jeong et al. 2011;Punkka and Bister 2015) and this does not agree with the seasonal cycle of ΔWG convective shown in Fig. 10.

Towards a new gust parametrization
The identified discrepancy in the mean seasonal cycle between the observed and simulated WG (see Sect. 4.2) suggests that an improved parametrization of WG across Sweden is required. For this reason, the WG parametrization was adjusted by tuning for each station for both the C turbulent of the ΔWG turbulent term and the C convective of ΔWG convective . Using multi-regression, for each station and using all time steps available, the function in Eq. (11) is fitted following the parametrized WG formulation of Eq. (5): with (11) Y = C turbulent ⋅ ΔWG friction+stability + C convective ⋅ ΔWS + shear + C unexplained Fig. 4, but for the mean diurnal cycle where C unexplained is a coefficient also calculated by multiregression which represents what the turbulent and the convective contributions cannot explain in the difference between observed WG and WS. As discussed in Sect. 4.2, for each station and at each hour, ERA5 outputs (e.g., friction velocity) of the closest grid point are used to calculate ΔWG friction + stability and ΔWS + shear . The different coefficients (C turbulent , C convective and C unexplained ) of Eq. (11) are calculated by different tuning procedures in order to evaluate the importance of each term in the parametrization:
• Tuning procedure 1-By using the fixed C convective = 0.6 (standard value suggested by Bechtold and Bidlot 2009), the coefficients C turbulent and C unexplained are calculated for each station. • Tuning procedure 2-The coefficients C convective and C unexplained are calculated with fixed C turbulent = 7.71 (as suggested by Panofsky et al. 1977). • Tuning procedure 3-C turbulent , C convective and C unexplained are all calculated by multi-regression.
In Fig. 11 the different tuned coefficients C turbulent , C convective and C unexplained are plotted, for each tuning procedure, against each respective station elevation. When plotting the calculated C turbulent against the station elevation, the where h is the station elevation. A and B are the coefficients of the best-fit curve. The coefficient values, obtained by interpolating the curve, are shown in Table 1. C convective does not display any elevation-dependency in its distribution, but its average value (both in tuning procedure 2 and tuning procedure 3) is much lower than the 0.6 suggested by Bechtold and Bidlot (2009). Overall, when an elevation-dependency is displayed in C turbulent , the average unexplained contribution of C unexplained is small (~ 0.1 m s −1 for tuning procedure 1 and ~ 0.3 m s −1 for tuning procedure 3) and its distribution does not vary with station elevation. In contrast to tuning procedure 2 where C turbulent is fixed, the unexplained contribution increases slightly with increase in elevation: this linear increase could be included in the turbulent term when C turbulent is allowed to vary with the elevation (as for tuning procedure 1 and tuning procedure 3).
Following the results discussed so far, we aim here to further improve the parametrization of WG by: (1) implementing an elevation-dependency in the ΔWG turbulent contribution through a C turbulent function which is able to take elevation differences among different regions (coast, inland and mountain) into account; and (2) tuning the ΔWG convective term through the adjustment of the C convective coefficients according to observed WG climate across Sweden.  (2009): (1) a parametrization (hereafter, turbulence function) where, following tuning procedure 1, an elevation-dependent C turbulent is implemented but no change to the convective contribution is applied; (2) a parametrization (hereafter, convection tuned) where, following tuning procedure 2, C turbulent is fixed to 7.71 and the convective contribution is adjusted for Sweden by tuning C convective ; (3) a parametrization (hereafter, turbulence function + convection tuned) where, following tuning procedure 3, C turbulent is expressed as a function of elevation and C convective is tuned. The main features of these four parametrizations (standard, tuned 1, tuned 2 and tuned 3) are summarized in Table 1. When an elevation-dependency is applied, Eq. (13) is used to express the C turbulent function. To tune the convective contribution, the mean value of the C convective coefficients calculated for each station is adopted as the C convective parameter in the parametrization. C unexplained is not implemented in any of the tested parametrizations for a better comparison with the standard parametrization (where no C unexplained is included), primarly because the physical meaning of this additional term cannot be explained.
The four different WG parametrizations are tested against the observed WG. Figure 12 shows the mean seasonal cycles of WG from observations and various parametrizations for 2013-2017. No considerable discrepancies are shown between the observed WG and WG predicted from the different parametrizations for the coast stations, although an average underestimation of ~ 0.5 m s −1 is shown. When a turbulence function is applied, the parametrized WG simulate the observed seasonal cycle better. Across inland and especially for the mountain regions, if no elevation-dependency is included (as in WG standard and convection tuned), the observed seasonal variability is strongly underestimated (up to ~ 2 m s −1 across mountain). On the other hand, when the convective contribution is tuned, the shape of the seasonal pattern of the parametrized WG better resembles the observed shape: the convective tuning can improve the intermonthly variability but it cannot get rid of the underestimation unless also an elevation-dependency is included (as in WG turbulence function + convection tuned).
Statistics to enable a detailed comparison are shown in Fig. 13. Mean Pearson's correlations are higher for both convection tuned and turbulence function + convection tuned parametrizations compared to the standard parametrization: when the convective gust contribution is tuned, the observed variability is improved, followed by the parametrized WG at all temporal scales. But for the convection tuned WG, there are no improvements in both mean RMSE and bias compared to the standard parametrization: the C convective tuning cannot improve the evident negative bias. This negative bias can be reduced, especially across the mountain region, when an elevation-dependency is implemented in the turbulent contribution. This is seen in the mean bias statistics of the turbulent tuned and turbulence function + convection tuned parametrizations. To summarize, Fig. 13 shows that: (1) by implementing an elevation-dependency in the turbulent contribution, the negative bias displayed by the standard (2) by tuning the convective term, higher correlation with the observed WG can be reached. Better performance in the simulation of WG can be achieved when both a turbulence function is implemented and the convective gust term is tuned. In comparison to the standard parametrization, the most appreciable improvements are achieved in the turbulence function + convection tuned parametrization: both higher mean correlation and lower average bias are shown for this new parametrization.

Tests for Norwegian stations
To test if the developed turbulence function + convection tuned parametrization could also be applied successfully to other regions, it was tested in Norway, situated alongside Sweden and characterized by a more complex topography (Fig. 1). In fact, Norway is dominated by the Scandes mountain ranges with an average elevation of 500 m a.s.l.: only one-fifth of Norway's total areas is situated below 150 m (https ://www.natio nsenc yclop edia.com/Europ e/Norwa y-TOPOG RAPHY .html; last accessed 8 May 2020). Its high plateaus are traversed by rugged mountains broken by numerous valleys, and the coastlines are characterized by the fjords, narrow inlets with steep sides created by glaciers. On the contrary, the Scandes only dominate central and northern Sweden, especially in the western areas close to the Norwegian border (Achberger et al. 2006). Coastal regions and southern Sweden are mostly flat and characterized by lowlands. Due to such topography differences, Norway was chosen to test the formulated WG parametrizations in more complex terrain.
The created elevation-dependent C turbulent function developed for Sweden and the mean C convective coefficient calculated for the Swedish stations are used to compute parametrized WG across Norway (hereafter, tuned − Sweden design). Note that the stations across Norway have been selected so that their elevation is not higher than around 900 m, i.e. the maximum elevation registered for the Swedish stations used in this study. In this way the C turbulent function, created using the Swedish stations in the elevation-range 0-900 m, should also be valid for the same elevation-range across Norway. Furthermore, the 15 stations across Norway are also used to design an elevation-dependent C turbulent and a tuned C convective in the same way that those coefficients have been developed for Swedish stations (Sect. 4.3): these different C turbulent and tuned C convective functions are used to parametrize WG (hereafter, tuned − Norway design). The detailed differences in the designed turbulent function and in the calculated mean C convective coefficient are shown in Table 2. Figure 14 shows the mean seasonal cycle of WG calculated using observations and the different parametrizations (standard, tuned − Sweden design and tuned − Norway design) across Norway. The observed seasonal cycle across Norway shows higher WG during winter than summer, with a peak in June. The standard parametrization underestimates measured WG, with varying differences for different months (from ~ 0.5 m s −1 of January to ~ 1.5 m s −1 of April). The mean seasonal cycle of the tuned WG follows the variation of the observed WG better than that of WG standard does, with an almost constant negative bias (~ 0.75 m s −1 for tuned − Sweden design and ~ 0.5 m s −1 for tuned − Norway design) throughout the year. These features are confirmed by Fig. 15: the tuned − Sweden design and the tuned − Norway design parametrizations perform better than the standard one. Among the three parametrizations, the tuned ones show Fig. 10 Contribution of the different WG parametrization terms to the mean seasonal cycle across different region for 2015-2017. In particular, in dark red the contribution of the observed WS is rep-resented, in light blue the one of the calculated ΔWG turbulent , and in yellow-green the added input of ΔWG convective . In black the mean seasonal cycle from the observed WG is shown the highest Pearson's correlation coefficients, with RMSE errors largely explained by a negative bias. Such negative bias is greater for WG tuned − Sweden design than WG tuned − Norway design and could be caused by the different durations of the measured gust. In fact, WG in Norwegian stations is calculated as the maximum 3 s wind speed in the previous hour, while across Sweden WG refers to the maximum 2 s wind speed in the previous hour. The elevation-dependent C turbulent function, designed by using Swedish stations with 2 s maximum gusts, could cause an additional underestimation when adopted to simulate WG across Norway, where WG is recorded as the 3 s maximum wind speed. But this hypothesis is not supported by the work of Safaei Pirooz et al. (2020, under revision), where the effect of gust duration on recorded gust wind speed is quantified. Safaei Pirooz et al. (2020, under revision) reports that, by increasing the gust duration from 2 to 3 s, both gust and peak factors decreases, meaning that lower gust wind speeds are recorded when employing a 3 s gust duration compared to shorter gust durations, such as 2 s. According to the analysis of Safaei Pirooz et al. (2020, under revision), a positive difference between WG tuned − Sweden design and WG tuned − Norway design (i.e., WG tuned − Sweden design greater than WG tuned − Norway design) should have been Fig. 11 Distributions of the C turbulent (left column), C convective (middle column) and C unexplained (right column) coefficients versus the measuring station's elevation in the tuning procedure 1 (top row), tuning procedure 2 (middle row), and tuning procedure 3 (bottom row). In red the linear regression which fit best to the distribution is shown. The dashed blue line is the mean value of the distribution. If a coefficient was fixed in the tuning procedure, its value is displayed as blue line detected across Norway if the difference in gust duration was the reason for the discrepancy. Therefore, the greater negative bias in the tuned − Sweden design compared to the tuned − Norway design WG may be explained by other factors such as the more complex topography across Norway.

Generalization of the new WG parametrization
The results have shown evident improvements in the parametrization of WG across Sweden by adding an elevationdependency to the turbulent gustiness and through calibration of the convective contribution. The validity of such improved formulation is confirmed as, even with calibrations using only the Swedish data, better WG simulations are also seen in Norway. For this reason, future studies should investigate how WG can be parametrized better in other regions when an elevation-dependent turbulent contribution and a convective tuning are implemented. Only in this way can the new gust parametrization designed in this study be reliably generalized and possibly be used to simulate WG for different regions or even globally. In order to evaluate the universality of such elevationdependent turbulent gusts, the C turbulent function should be created and tested for a greater number of stations with a larger range of elevations: here only three mountain (i.e., elevation higher than 750 m) stations were available across Sweden, with a maximum elevation of only ~ 900 m a.s.l..
It is worth mentioning that this study has improved the WG parametrization currently used by ERA5 through a statistical post-processing technique, as explained in Sect. 3.2. Such technique has evidenced the need of (1) including an elevation-dependency to the turbulent

Possible explanation of the elevation-dependency in turbulent gustiness
This study has shown the need for implementing an elevation-dependent C turbulent function to more realistically simulate WG across both Sweden and Norway. But what effect can explain such an elevation-dependency? As mentioned in Sect. 4.2, the formulation of Panofsky et al. (1977) to express turbulent gusts requires a horizontally homogeneous surface:  all the data analyzed by Panofsky et al. (1977) were in fact obtained over such flat surfaces where the Monin-Obukhov similarity hypothesis is valid (Arya 2001). In line with this limitation, Panofsky et al. (1977) clearly stated that "it is quite possible that mesoscale terrain features can influence the scale of horizontal velocity components". Over the land grid of a climate model, such as that used for ERA5, the condition of homogeneity cannot be guaranteed. Any atmospheric model, with its inherent spatial discretization, has to smooth the surface physical properties such as orography (Jiménez and Dudhia 2012). The approximation of smoother topography is less accurate especially over mountain areas, where the heterogeneity driven by valleys and mountains becomes relevant. Unresolved topographic features can introduce biases in wind simulations across complex terrain regions, as shown by Jiménez and Dudhia (2012) across the north-east of the Iberian Peninsula. By including unresolved terrain characteristics through the use of the concept of a momentum sink term, and the use of the standard deviation of the subgrid-scale orography as well as the Laplacian of the topographic field, Jiménez and Dudhia (2012) were able to improve WS simulations over a complex mountain terrain region. Following such results, the elevation dependency of the turbulent gusts found in the current study across both Sweden and Norway points to the need of including unresolved topography in the WG parametrization. This can be better achieved in future studies only if the possible mechanisms by which gustiness is affected by adjacent topography are explained and the terrain influence can be quantitatively addressed, similarly to what has been done by Zhao et al. (2020) in relation to the roles of moisture budget under the influence of topography on summer extreme precipitation over Northern China.

Proxies of convective gustiness
Convective wind gusts are generated by rapidly descending air masses in thunderstorms or downdrafts (Holleman 2001). When impacting the ground, the descending air is deflected by the surface and causes wind gusts as outflow. Convective gusts can be generated in mesoscale convective systems, which are organized assemblages of thunderstorms where the cumulative effects of convective updrafts and downdrafts interacting on scales of 100 km generate distinct mesoscale circulations persisting for several hours (Haberlie and Ashley 2019). A sheared environment is needed to increase the potential for a mesoscale convective system to produce severe surface winds through a stronger and more organized convective structure (Cohen et al. 2007). In general, the vertical wind shear can organize deep convective systems (e.g., thunderstorms, squall lines, and other mesoscale features) and extend their lifetime (Anber et al. 2014). For this reason, forecasters have adopted the gradient wind, as lowlevel wind shear, to estimate the maximum gusts likely to occur near the surface (Bradbury et al. 1994;Nakamura et al. 1996). Based on this understanding, Bechtold and Bidlot (2009) parametrized gusts from deep-convection using lowlevel wind shear (expressed as WS 850 WS 950 ) as it enhances downdraft and is thus related to the maximum gust speed. This study has improved the gust forecasting ability also by calibrating of the convective contribution with observational data. Statistically, such a tuning procedure works, but it cannot be considered adequate: it should be explored a better formulation of ΔWG convective which expresses more realistically the relevant physical processes at the origin of convective gustiness across Sweden. This is supported by the fact that ΔWG convective seasonality across Sweden is not well represented by the current parametrization (where maximum gustiness is related to low-level wind shear) as it is not in line with the seasonal occurrence of deep-convection situations and mesoscale convective systems at the origin of convective gustiness (see Sect. 4.2). Low-level wind shear expressed as the difference between 850 and 950 hPa may not be suitable for a region like Scandinavia: convective gusts may be calculated in relation to a different definition of low-level wind shear or by different proxies of maximum gustiness. Indeed, significant wind shear at mid-levels of the troposphere can also influence lower-level flow strength through supercell development by removing precipitation from updrafts and inducing vertical perturbation pressure gradients (Rotunno and Klemp 1982). As a result, a deeperlayer shear (as 0-6 km) is often needed to assess severe wind events (Johns and Doswell 1992;Pilortz et al. 2016). This is the case in work by Cohen et al. (2007), where wind shear in a shallow layer (0-2 km) is not found to be as good as the 0-6 km and 0-10 km shears in discriminating between different mesoscale convective systems. Future studies should investigate the possibilities of forecasting surface gusts by using a deeper wind shear or even by relating convective gusts to other relevant meteorological variables, such as convective available potential energy (CAPE), midlevel lapse-rate or vertical differences in equivalent potential air temperature (Cohen et al. 2007).

Summary and conclusions
In this study, the new ERA5 reanalysis product has been compared with the hourly WS and WG observations across Sweden for the period 2013-2017. It was found that ERA5 reproduces both the WS and WG measurements more realistically than its predecessor ERAINT, most likely due to higher resolution, better model physics, more data assimilated, and its more advanced assimilation method. However, evident discrepancies are still found across the inland and mountain regions. For this reason, the WG parametrization was further investigated. In ERA5 wind gusts are parametrized as the sum of three terms: (1) the WS; (2) the turbulent gustiness ΔWG turbulent , which expresses the effects of the surface friction and the BL stability; and (3) the convective contribution ΔWG convective to include the gustiness in deep convection. Improved formulations of ΔWG convective were explored by: (1) implementing an elevation-dependent coefficient in the turbulent contribution to include elevation differences among the stations; (2) adjusting the convective gust contribution through a tuned coefficient C convective . Just by including this elevation-dependency, noticeable improvements were seen, especially over the mountain regions: the elevation-dependent function can explain most of the negative bias in the simulated WG. Even closer agreement between the observed and parametrized WG were shown when both the elevation-dependent function and the tuned convection term were implemented.
To test if these implementations could forecast WG better in other regions other than Sweden, the new formulated WG parametrization was tested using 15 stations across Norway during 2015-2017. Here the topography is more complex, dominated by the Scandes mountain range that runs through the whole country. For Norway, when the elevation-dependent turbulent coefficient and the tuned convective contribution (both calibrated with Sweden data) are added, the observed WG variability is more realistically simulated, although it is still negatively biased. This bias is reduced, but not eliminated, when the elevation-dependency and the tuning of C convective are designed using Norwegian stations. The difference in bias between Sweden-tuned and Norway-tuned parametrized WG cannot be explained by the differences in gust duration on the recorded WG in Sweden and Norway. Such errors must be due to other reasons such as the more complex topography over Norway.