The contribution of large-scale atmospheric circulation to variations of observed near-surface wind speed across Sweden since 1926

This study investigates the centennial-scale (i.e., since 1926) variability of observed near-surface wind speed across Sweden. Results show that wind speed underwent various phases of change during 1926–2019, i.e., (a) a clear slowdown during 1926–1960; (b) a stabilization from 1960 to 1990; (c) another clear slowdown during 1990–2003; (d) a slight recovery/stabilization period for 2003–2014, which may continue with a possible new slowdown. Furthermore, the performance of three reanalysis products in representing past wind variations is evaluated. The observed low-frequency variability is properly simulated by the selected reanalyses and is linked to the variations of different large-scale atmospheric circulation patterns (e.g., the North Atlantic Oscillation). However, the evident periods of decreasing trend during 1926–1960 and 1990–2003, which drive most of the stilling in the last century, are missing in the reanalyses and cannot be realistically modeled through multiple linear regression by only using indexes of atmospheric circulation. Therefore, this study reveals that changes in large-scale atmospheric circulation mainly drive the low-frequency variability of observed near-surface wind speed, while other factors (e.g., changes in surface roughness) are crucial for explaining the periods of strong terrestrial stilling across Sweden.


Introduction
Near-surface wind speed (hereafter, NSWS) is a key factor in transferring heat, moisture, energy, and momentum between the Earth's surface and the atmosphere (Abhishek et al. 2010). By controlling the evaporation demand, surface winds can alter the hydrological cycle and partly affect agriculture productivity (Rayner 2007;McVicar et al. 2012 regulates as well the accumulation and dispersion of air pollutants near emission sources, such as traffic in urban areas (Grundström et al. 2015); it can also contribute to the cooling in urban heat islands (Bing et al. 2021), and it strongly impacts soil erosion over many dry land regions ). In addition, as the development of renewable energy resources is central to energy scenarios that help keeping warming below 2 °C (IPCC 2014), society relies on electricity production from wind farms when it comes to the reduction of greenhouse gas emissions (Zeng et al. 2019). For example, in Sweden, wind power stood for about 16% of the overall energy production in 2020 (Berard 2021).
For all the mentioned reasons, it is crucial to understand how NSWS has been and will be affected by a warmer climate, so that society can adapt to the expected new wind scenarios (IPCC 2013). Therefore, over the last few decades, various studies have investigated multidecadal changes in observed NSWS across various areas of the globe (e.g., Azorin-Molina et al. 2014;Laapas and Venäläinen 2017). It has been reported a general slowdown in NSWS over land in most northern midlatitude regions during the last 30-50 years (McVicar et al. 2012). Such a general slowdown in terrestrial winds, which has been named "stilling" by Roderick et al. (2007), differs from the opposite increase in NSWS revealed over large parts of the oceans (Zheng et al. 2016). However, recent studies have also evidenced a recovery in the observed terrestrial NSWS decline during the last decades (Kim and Paik 2015;Zhang and Wang 2020). Such a break in the terrestrial stilling became prominent since around 2010 across the Northern Hemisphere, especially in Europe, North America, and East Asia (Zeng et al. 2019).
Even if the exact reasons behind the terrestrial stilling and its reversal are unclear, different possible causes have been proposed. Vautard et al. (2010) showed that the increase in surface roughness (e.g., urbanization, land-use changes, forest growth) may be the main reason behind the terrestrial decrease in NSWS. But this differs from what was found by Zeng et al. (2018), who argued that the increase in vegetation cover has had only a limited influence on the wind stilling. Large-scale atmospheric circulation has been identified by many studies as a key driver for NSWS variations and changes (Wu et al. 2018). In fact, in a changing climate, regional warming differences control variations of regional surface pressure gradients, which modulate circulation patterns (Lin et al. 2013). For instance, phase changes of the North Atlantic Oscillation (Azorin-Molina et al. 2018a;Minola et al. 2021) and of the Pacific Decadal Oscillation (Zeng et al. 2019) have been linked to the reversal in the terrestrial stilling. The aging of measuring instruments has also been suggested as a potential cause of the NSWS decrease (Azorin-Molina et al. 2018b).
Previous studies have already investigated NSWS changes and variations across Sweden. In particular, Minola et al. (2016) looked at the multidecadal variability in observed NSWS during 1956-2013. An overall statistically significant slowdown was found, but strong differences in the seasonal trends were also found (i.e., weak increase in winter and strong stilling in spring, summer, and autumn). Such seasonal differences could be explained by the impact of large-scale atmospheric circulation, in particular by the influence of the North Atlantic Oscillation, especially when it comes to modulate winter NSWS variability. In addition, Minola et al. (2021) investigated NSWS changes across Sweden for the recent past two decades (i.e., 1997-2019), which were only partly covered by Minola et al. (2016). Here it is shown that the stilling ceased in 2003, and no clear trend can be detected afterwards. Such stilling-reversal pattern was linked to changes in the North Atlantic Oscillation. Furthermore, this study revealed that changes in surface roughness (e.g., changes in forest cover) contributed to the general slowdown detected for 1997-2019.
As NSWS varies on multiple temporal scales, sufficiently long data series are necessary for a meaningful assessment of climate variability (IPCC 2021). Even though climate journals may store valuable old not-digitized climate data , studies on NSWS changes have only focused over the last ~60 years as reliable, continuous, and easy-to-access observations lack over a longer time window. Because NSWS observations are almost absent at the centennial scale, alternative datasets, such as twentieth century reanalyses (i.e., developed to cover the entire twentieth century), have been largely used to investigate the centennial-scale variability in NSWS (Bett et al. 2017;Shen et al. 2021a). In fact, by using a forecast model in which information from observations are assimilated (Dee et al. 2011), reanalyses provide spatially complete and physically coherent simulated data. However, any climate-variability result obtained by using reanalysis products must be evaluated against in situ measurements data as the reanalysis performance is strongly dependent on the considered time period and the selected study area (Ramon et al. 2019;Yu et al. 2019).
For all the mentioned reasons, this study aims at investigating centennial-scale (since 1926) variations and changes in observed NSWS across Sweden. No other dataset of observed NSWS has been investigated before for such a long time period. By comparing observed series against the ones from reanalysis products, the performance of current climate reanalyses in reproducing past changes in NSWS across Sweden on a centennial time scale is also evaluated. In addition, the impact of large-scale atmospheric circulation on the detected variations is explored by quantifying its contribution to the NSWS variability.

NSWS observations
NSWS observations across Sweden are available for downloading at: https:// www. smhi. se/ data/ utfor skaren-oppna-data/ (last accessed 9 March 2023). At this link, it is possible to access only wind data which has been measured, digitized, and later stored at the archive servers of the Swedish Meteorological and Hydrological Institute (SMHI). Reliable and continuous NSWS measurements that can be downloaded here start only in the 1950s at the earliest (Minola et al. 2016). But there are still valuable wind observations that have not been digitized yet, which are stored on paper climate journals at the SMHI archive in Norrköping. Therefore, to create the longest-available century-long NSWS dataset, through the work of Engström et al. (2022), wind measurements were rescued and digitized by scanning old weather journals from 1920 to 1940. The resulting dataset is constituted by rescued NSWS series from 13 measuring stations, which contain observations since the early 20s. To remove any non-climatic change point, a four-step homogenization procedure was later applied to the monthly mean observed NSWS series. More details about how the observed series were homogenized can be found in Zhou et al. (2022). However, the rescued and homogenized 13 stations do not all cover the same time period (e.g., only two stations are still active today, most stations have observations until the 90s) and contain large periods of missing values. Therefore, we select only 7 stations that cover NSWS observations for the 71-year 1926-1996 period with a maximum of 5% of missing data (i.e., 42 months). The 1926-1996 time window was selected as it was identified as the time period that most of the stations can cover with not too many missing values. Figure 1 shows the location of the 7 stations selected in the NSWS dataset which will be analyzed here. Unfortunately, all the stations from the west coast of Sweden were discarded (as they contain long periods of missing data).
Since the rescued NSWS dataset covers only 1926-1996, in order to investigate NSWS changes also after 1997, we complement this dataset with NSWS measurements from two previous studies: (i) homogenized NSWS observations for 1956-2013 from Minola et al. (2016); and (ii) homogenized NSWS observations for 1997-2019 from Minola et al. (2021). By comparing those three datasets, we can thus investigate changes and variations in NSWS across Sweden during 1926-2019, i.e., over almost a century. To make an accurate comparison, from the datasets of Minola et al. (2016) and Minola et al. (2021) we only select the same stations from the 1926-1996 dataset. If the same station does not appear in those two datasets, we select the nearest station available (see Fig. 1). For Holmögadd station, the nearest station in the 1956-2013 dataset is Bjuröklubb (at 106 km distance), which is already included as the nearby station. Therefore, for the dataset 1956-2013, only 6 stations are used for representing the regional mean.

NSWS from reanalyses
This study compares observed NSWS with wind outputs of three different reanalyses: (a) the European Center for Medium-Range Weather Forecasts (ECMWF) twentieth century reanalysis (hereafter, ERA20C); (b) the fifth generation ECMWF reanalysis (hereafter, ERA5); and (c) the National Oceanic and Atmospheric Administration-Cooperative Institute for Research in Environmental Sciences (NOAA-CIRES) twentieth century reanalysis, version 2c (hereafter, 20CR). ERA20C is the ECMWF's first atmospheric reanalysis of the twentieth century, covering the 1900-2010 time period (Poli et al. 2016). By assimilating only observations of surface pressure and surface marine winds, it models various climatic variables at a horizontal resolution of ~125 km, with a 3-hourly temporal resolution. Instead, ERA5 is the latest modern reanalysis product of ECMWF (Hersbach et al. 2020). It delivers hourly outputs at a horizontal grid spacing of ~31 km from 1950 to the present. Contrary to ERA20C, which assimilates only surface pressure and surface marine wind observations, ERA5 combines a large amount of different historical observations using advanced modeling and data assimilation systems. It also assimilates terrestrial vertical wind profiles from satellites, and radio-and aircraft-sondes, but it does not include NSWS observations over land as they cannot be fully interpreted by the data assimilation system (Dee et al. 2011). 20CR is a global gridded reanalysis produced by the collaboration between NOAA's Physical Sciences Laboratory and CIRES at the University of Colorado (Compo et al. 2011). It assimilates only surface observations of synoptic pressure into the NOAA's Global Forecast System and prescribes sea surface temperature and sea ice distribution for estimating various climate variables since 1836 on a 2.0° × 2.0° latitude-longitude grid.

Indexes of atmospheric circulation patterns
To explore the possible influence of large-scale atmospheric circulation on the NSWS variability across Sweden, the relationship between wind series and different atmospheric circulation modes is analyzed. In particular, here we look at these 4 climate indexes: (1) the North Atlantic Oscillation (hereafter, NAO) index; (2) the Arctic Oscillation (hereafter, AO) index; (3) the Scandinavian Pattern (hereafter, SCA) index; and (4) the East Atlantic Pattern (hereafter, EA) index. NAO drives the shift of the Atlantic storm tracks across Europe, thus affecting how stormy the weather conditions are in Northern or Southern Europe (Hurrel et al. 2003). The NAO series for this study is obtained from https:// cruda ta. uea. ac. uk/ cru/ data/ nao/ (last accessed 9 March 2023). AO affects the north-to-south location of the storm-steering mid-latitude jet stream (Scott 2021). The centennial-scale AO series (i.e., until 2002) is downloaded from https:// www. atmos. colos tate. edu/ ~davet/ ao/ Data/ ao_ index. html (last accessed 9 March 2023). To be able to cover the more recent years (i.e., from 2002 to the present), AO index for 1950-2019 was obtained from https:// www. cpc. ncep. noaa. gov/ produ cts/ precip/ CWlink/ daily_ ao_ index/ ao. shtml (last accessed 9 March 2023). The SCA pattern influences the exit region of the Atlantic jet stream from its climatological mean position, thus determining the preferred region of cyclone growth (Blackburn and Hoskins 2001). EA may play a role in positioning the North Atlantic storm track by modulating the location and strength of the NAO dipole (Moore et al. 2011). Centennial-scale SCA and EA series are retrieved from Comas-Bru and Hernández (2018) (https:// doi. panga ea. de/ 10. 1594/ PANGA EA. 892768? format= html# downl oad, last accessed 9 March 2023).

Trend analyses
Trends in NSWS are represented by the slope of the applied regression analysis (i.e., linear regression) and are expressed as changes per decade (i.e., dec −1 ). The magnitude of the linear trends is calculated using the non-parametric Sen's method (Gilbert 1987). Their significance is defined using the modified Mann-Kendall test (for considering the effect of statistically significant autocorrelation; Hamed and Ramachandra Rao 1998) and is reported at three different levels of significance: (1) significant at p < 0.05; (2) significant at p < 0.10; and (3) non-significant at p > 0.10 (Azorin-Molina et al. 2014). On top of the linear trends, the Gaussian low-pass filter with a 15-year window is used to show the low-frequency variability of the time series. When possible, NSWS series are expressed as anomaly series: in this way, when calculating regional means, more windy series do not dominate the regional NSWS series (Azorin-Molina et al. 2014). For example, for the 1926-1996 dataset, anomalies are calculated as the deviation from 1926 to 1996. But when it comes to the comparison between the datasets 1926-1996, 1956-2013, and 1997-2019, NSWS series cannot be expressed as anomalies as there is no common time period for all the series. Analyses are carried out for annual and seasonal series.

Statistics of comparison
To mathematically evaluate the agreement between different series, the following statistical tests are used: (1) Pearson's correlation coefficient (hereafter, r), to measure the degree of association (i.e., linear relationship; Gibbons and Chakraborti 2003); (2) root mean squared error (hereafter, RMSE), to mathematically express the vicinity between two series (Von Storch and Zwiers 1999); and (3) bias, to identify the tendency to constant deviate of a realization compared to another.

Multiple linear regression
To quantify the contribution of large-scale atmospheric circulation to NSWS variability, variance based on multiple linear regression (hereafter, MLR) is used (Pedroni 1999;Shi et al. 2019). In particular, by using the indexes of teleconnection patterns (one or more of them) as independent variables and the NSWS as the true realization, regression coefficients are calculated based on multiple regression theory. Those regression coefficients, in combination with the circulation indexes, are then used to calculate the simulated NSWS, which expresses NSWS variability only driven by large-scale atmospheric circulation variations. The match/mismatch between raw and reconstructed NSWS series is quantified through the coefficient of determination R 2 (Von Storch and Zwiers 1999). Under the hypothesis that the two series are equal, R 2 quantifies the ability of the simulated NSWS series to explain the variation of the raw NSWS series. R 2 values closer to 1 indicate a higher level of association. Figure 2 shows the variability of the mean (i.e., average over all stations) annual and seasonal NSWS anomaly series during 1926-1996, while Table 1 summarizes their trends. Annually, an overall negative trend of −0.11 m s −1 dec −1 (statistically significant at p < 0.05) has been found during 1926-1996, even though periods of different signs of change occurred. In particular, three phases of NSWS changes can be identified during 1926-1996 (see Table 1): (a) a clear slowdown of −0.20 m s −1 dec −1 (p < 0.05) during 1926-1960; (b) a stabilization from 1960 to 1990 (+0.01 m s −1 dec −1 , non-significant at p < 0.05); and (c) the possible start of a new slowdown since 1990. Seasonally, winter is the season that shows the highest interannual variability, while summer is the lowest one (Fig. 2). All seasons display an overall significant (p < 0.05) slowdown of ~ −0.11 m s −1 dec −1 during 1926-1996 (Table 1), which is caused by a distinct decrease for 1960-1970, followed by a period of stabilization/recovery until ~ 1990. Such a decreasing trend is more evident for summer and autumn, while for winter and spring, even if a general slowdown appears, the NSWS series are also dominated by interannual variabilities and by cycles of changes. To notice that the great negative NSWS anomaly observed in winter of 1996 strongly affects the winter (and therefore also the annual) trend for 1990-1996, which appears to be strongly negative, but not significant at p < 0.05.

NSWS variability since 1926
To better understand the NSWS variability over Sweden during the last century, the 1926-1996 series of this study are compared with the ones from Minola et al. (2016) and Minola et al. (2021), able to cover the more recent decades. Therefore, by adding the 1956-2013 and 1997-2019 datasets to the rescued dataset, we can investigate NSWS variability from 1926 until 2019. As shown in Fig. 3, the 1926-1996 and the 1956-2013 mean (i.e., averaged over all the stations in the dataset) series display a similar variability during the common period 1956-1996: r have high values both annually and seasonally (e.g., for annual series, r is 0.85, statistically significant at p < 0.05). The 1956-2013 and 1997-2019 mean NSWS series also show a high correlation for the common 1997-2013 time period (e.g., annual r of 0.72, statistically significant at p < 0.05).
Even though the three datasets correlate well with each other, they display differences in the magnitude of NSWS. Such biases may be related to the different homogenization techniques adopted and the different reference series used to detect breakpoints and adjust them. However, the high correlation values confirm that the series complement each other and can be used all together for investigating NSWS variabilities and changes during the extended period of 1926-2019. Annually, the stilling detected from 1990 continues until around 2003; it is then followed by a period of stabilization and weak recovery from such a strong decreasing trend. Those phases of change can be identified for all the seasons, even in winter and autumn under a stronger interannual variability   during 1926-1996, 1956-2013, and 1997-2019. Annually and seasonally, both the reanalyses do not capture the strong decrease observed during 1926-1960. Instead, for 1900-1990, they simulate a positive increase in NSWS (e.g., for annual series in ERA20C, 0.03 m s −1 dec −1 during 1900-1990, significant at p < 0.05). Similarly, also the evident stilling in observed NSWS for 1990-2013 does not have the same magnitude in the simulated trend of ERA20C and 20CR: in fact, from 1990 to 2010, the reanalyses also show a general decline, but such slowdown is not as evident and strong as the one in the observed dataset. Table 2 further explores the performance of the two reanalyses in representing observed NSWS by showing various statistics for comparison. While ERA20C tends to underestimate observed wind speed and 20CR overestimates NSWS, reanalyzed NSWS series well correlate with the observed series during winter, autumn, and summer (lower r values annually and during summer), even if they do not display the observed strong stilling during 1926-1960 and 1990-2003. Overall, although ERA20C and 20CR cannot reproduce the two stilling periods, they can represent well the overall interdecadal variability with its fluctuations in the change. For instance, the low-frequency variations of observed NSWS in spring during  are properly simulated by all the reanalyses. Figure 4 also shows the variability of NSWS from ERA5 during 1950-2019 in comparison with the observed one for 1926-1996, 1956-2013, and 1997-2019, while Table 2 investigates the performance of this reanalysis in simulating observed NSWS variabilities. Like ERA20C and 20CR, this reanalysis cannot simulate the strong negative trend for 1926-1960 and 1990-2013. Instead, ERA5 well captures the weak variations occurring Table 1 Trends (m s −1 dec −1 ) of the mean (i.e., average over all stations across Sweden) annual and seasonal NSWS anomalies for 1926-1996 and the sub-periods 1926-1960, 1960-2990, and 1990-1996. Statistically significant trends are shown in boldface for p < 0.05 and in italic for p < 0.10 1926-1996 1926-1960 1960-1990 1990-1996  Overall, even if ERA20C, 20CR, and ERA5 differ in the observation data that they assimilate (see Sect. 2.2), all the reanalyses well agree in how they simulate the variations in NSWS over Sweden (Fig. 4 and Fig. S1 in the supplementary material). Annually and seasonally, correlation during the common time period 1950-2010 is always higher than 0.75 and significant at p < 0.05 (Table S1 in the supplementary material). As shown in Fig. 4, the reanalyses do not show the significant and strong decrease in NSWS as displayed by the observation dataset. Instead, they show periods of variations: those phases of change are similar in ERA20C, 20CR, and ERA5. Even if the correlation is high, ERA20C and ERA5 differ in a constant bias of ~0.6 m s −1 , ERA5 and 20CR in ~0.6 m s −1 , and ERA20C  (Table S1 in the supplementary material). Overall, ERA20C simulates less strong wind compared to ERA5, while the modeled wind of 20CR is the strongest among the dataset here. From the comparison of Table 2, it is evident that ERA5 performs better than ERA20C and 20CR in simulating observed NSWS variability, most likely due to its higher resolution (which helps in the comparison between one grid point and single in-situ data), better model physics, more data assimilated, and its more advanced assimilation method (see Sect. 2.2).  Table 2 Annual and seasonal statistics for comparison between mean (i.e., average over all stations across Sweden) observed NSWS and reanalyzed NSWS from ERA20C for 1926-1996(and 1956), 20CR for 1926-1996(and 1956), and ERA5 for 1950-1996(1956. In particular, r, RMSE, and bias are shown. r coefficients statistically significant at p < 0.05 are shown in boldface  Table 3 Annual and seasonal Pearson's correlation coefficients between NSWS anomalies and the NAO index, the AO index, the SCA index, and the EA index for 1926-1996, 1956-2013, and 1997-2015 The relationship between observed NSWS and the 4 different indexes of atmospheric circulation is investigated by calculating the Pearson's correlation for 1926-1996 annual and seasonal series (Table 3). NAO displays a significant (p < 0.05) positive correlation annually and for all the seasons, except for summer. When plotting the NAO index versus the observed NSWS (Fig. S2 in the supplementary material), it is evident that NAO well follows the low-frequency variability of observed wind (especially for winter, when r is The main reason for the mismatch between AO and wind speed series is that the periods of strong stilling in observed NSWS do not appear in the AO variations. Like NAO and AO, EA correlates well with the observed NSWS series for winter, spring, and autumn. SCA is the only circulation pattern that displays a significant (p < 0.05) negative correlation (−0.27) for summer, while for the other seasons and annually, r values are low and not significant. Table 3 further explores the correlation between the 4 indexes of atmospheric circulation modes and the observed NSWS anomaly series of Minola et al. (2016Minola et al. ( , 2021Minola et al. ( ) for 1956Minola et al. ( -2013Minola et al. ( and 1997Minola et al. ( -2015  large-scale atmospheric circulation may play a crucial role in explaining the changes and variations of observed wind, especially when it comes to its low-frequency variability.

Reconstructing NSWS variability
Given the strong influence of large-scale atmospheric circulation on NSWS variability (see Sect. 4.3), we use all 4 indexes of atmospheric circulation modes (i.e., NAO, AO, SCA, and EA) to reconstruct observed NSWS using the MLR. The reconstructed observed NSWS (hereafter, RWS-obs) is calculated using all 4 indexes as the combination of both NAO, AO, SCA, and EA can explain a higher variance (R 2 ) of the observed NSWS compared to the use of just one index or a combination of two or three of them (e.g., R 2 up to 0.44 for spring). Figure 5 compares the mean annual and seasonal series of observed NSWS and RWS-obs for 1926-1996. Annually, the RWS-obs captures the interannual variability and the small interdecadal fluctuations (e.g., variability with a period of change of around 5 years, as during 1980-1990) of the observed NSWS series. But RWS-obs does not display the clear slowdown for , which therefore cannot be explained by the large-scale atmospheric circulation changes. Seasonally, a similar pattern can be noticed: RWS-obs series follow the cyclical changes of observed NSWS, but do not simulate periods of evident decreasing trend, such as the 1926-1960 slowdown. For example, for autumn, the second period with a marked slowdown starts earlier than for the other seasons (1985rather than 1990see Fig. 3); when looking at the period 1926-1996, it is already possible to identify 11 years (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)  Recalling what is seen in Sect. 4.2, ERA20C, 20CR, and ERA5 cannot reproduce as well the evident observed decrease during 1926-1960 and 1990-2003. But all these reanalyses go along with the observed low-frequency variability. This resembles the behavior of the RWS-obs calculated using the teleconnection-pattern indexes, which also does not simulate the periods of significant decrease. Therefore, we also reconstruct the NSWS of ERA20C (hereafter, RWS-ERA20C), ERA5 (hereafter, RWS-ERA5), and 20CR (RWS-20CR). Those RWS series well reproduce the variability of NSWS in ERA20C (R 2 = 0.64), 20CR (R 2 = 0.62), and ERA5 (R 2 = 0.58) (Fig. 6). Instead, RWS-obs can only explain 0.12 of the variance of observed NSWS (mostly because RWS-obs does not capture the 1926-1960 slowdown, not-driven by large-scale atmospheric circulation changes). This confirms that all the selected reanalyses can simulate NSWS variations driven by large-scale atmospheric circulation, but cannot include the 1926-1960 and 1990-2003 periods of marked stilling, which may not be caused by changes in large-scale atmospheric circulation patterns.
To further confirm that the strong 1926-1960 decreasing trend is not driven by largescale atmospheric circulation changes, the mean observed NSWS series has been detrended (hereafter, obs_detrend), which means that the 1926-1960 negative trend has been removed from the mean series. The reconstructed NSWS series for obs_detrend (hereafter, RWS-obs_detrend) well follows the variability of the NSWS obs_detrend series, with a R 2 of 0.44, much greater than the R 2 of 0.12 for RWS-obs (Fig. 7). RWS-obs_detrend agrees as well with RWS-ERA20, RWS-ERA5, and RWS-20CR. This further confirms that the lowfrequency fluctuations in NSWS driven by the large-scale atmospheric circulation variability are captured by all the reanalyses. But the strong decreasing trends (e.g., −0.24 m s −1 dec −1 for 1926-1960, significant at p < 0.05), which were observed both during 1926-1960 and 1990-2003 and not simulated by the reanalysis products, should not be driven by changes in the large-scale atmospheric circulation. (d) a slight recovery/stabilization period for 2003-2014, which may continue with a possible new slowdown. When comparing the wind variability with the one of different indexes of atmospheric circulation patterns (i.e., NAO, AO, SCA, and EA), only the low-frequency variability of observed NSWS correlates well with the large-scale circulation variations: large-scale atmospheric circulation changes cannot explain the strong stilling during 1926-1960 and 1990-2003. Instead, they can explain small fluctuations, as different phases of changes in NSWS observed after 2003. In a similar way, also the reanalysis products (i.e., ERA20C, 20CR, and ERA5) do not simulate the robust decreasing trends during 1926-1960 and 1990-2003, but, instead, can only reproduce low-frequency variations. In fact, the 1926-1996 reconstructed NSWS, calculated by multiple regression using the different indexes of atmospheric circulation patterns, can explain only 0.12 of the observed wind variability. But when the strong decreasing trend during 1926-1960 is removed, the explained variability of the reconstructed wind increases to 0.44: large-scale atmospheric circulation changes are the most likely main drivers for the relative-small fluctuations of NSWS occurring over the past 100 years, but they cannot explain the significant decreasing trends of wind speed.

Summary and discussion
This study highlights that the 1926-1960 and 1990-2003 periods of evident slowdown drive most of the stilling in the last century. Those negative trends appear in every season, showing that the reason behind such stilling is not something occurring seasonally. In addition, the robust trends of those periods are not found in the reanalysis products: as reanalyses have proven to be able to properly simulate large-scale atmospheric circulation variations (Poli et al. 2016), such evident decreasing trends cannot be fully attributed to changes in atmospheric circulation patterns (Torralba et al. 2017). Instead, large-scale atmospheric changes can only explain the low-frequency variations of NSWS. In particular, no single index of atmospheric circulation pattern can fully reconstruct the relative-small variations in NSWS. Each index may be able to better explain the NSWS variability to some extent at a specific season (e.g., strong NAO influence during winter; Minola et al. 2016). But the indexes all together generate a higher amount of explanatory power in the reconstructed wind series (Zeng et al. 2019). This is because different teleconnections may interact with each other (e.g., EA modulates the location and the strength of the NAO dipole; Moore et al. 2011), and are complementary indicators in explaining climate changes and variations (Vicente-Serrano et al. 2016). Therefore, the low-frequency variability of NSWS cannot be simply linked to just a single atmospheric circulation pattern, but it should rather be resolved by the combined effects of variations in various teleconnection patterns (Shen et al. 2021a). Shen et al. (2021a) found that centennial-scale changes of NSWS over China are mainly driven by the combined effects of large-scale ocean-atmosphere circulations. Instead, across Sweden, large-scale atmospheric circulation changes can only explain the lowfrequency variations of NSWS, with the internal variability not being the reason behind the clear slowdown during 1926-1960 and 1990-2003. At the centennial scale, variations in large-scale atmospheric circulation patterns seem to drive cyclical decadal changes in NSWS, while the contribution of other factors (i.e., changes in surface friction) to the stilling across Sweden seems to be considerable (as it may explain the periods of evident decreasing trends). This is in line with what was found by model simulations at a global scale, where internal variability acted mainly to increase the NSWS from 1968 to 2014 through cyclical phases of change, while the contributions of other factors to the global terrestrial stilling after the 1960s were considerable (Shen et al. 2021b). Among the possible factors which could have affected past wind changes across Sweden, Minola et al. (2021) have already investigated the contribution of changes in forest cover to the observed wind variability during 1997-2019. Unfortunately, this study cannot attribute the observed stilling to such factors. For example, the lack of orthophotos at a centennial scale does not allow quantifying the contribution of forest cover changes to the observed trend for the time period considered in this study. Other factors which can be included among the possible factors at the origin of the evident periods of stilling are land-use change driven by urbanization , anthropogenic aerosol emissions (Li et al. 2016), and aging of measuring instruments (Azorin-Molina et al. 2018b).
This study has also shown that current reanalysis products are not able to properly reproduce the actual changes and variations in observed NSWS. In fact, wind series obtained from reanalysis datasets only display decadal low-frequency variations, but do not reproduce the evident slowdown during 1926-1960 and 1990-2003. This is in line with what is shown by Torralba et al. (2017) and Minola et al. (2021), which evidenced the inability of reanalysis datasets in reproducing the observed decline of surface winds as they do not include some drivers of wind speed variability such as changes in surface roughness. To overcome the lack of global coverage and long historical records of wind observations, researchers have largely relied upon reanalysis data for wind resource assessment to estimate revenues of electricity production from wind farms (Pryor et al. 2009;Holt and Wang 2012). Therefore, it is necessary to careful weight their conclusions as those analyses based on reanalysis products may have overestimated wind statistics as the wind series they have used may not reflect actual changes in surface winds (i.e., they do not include the past observed stilling).
To conclude, as already revealed by Minola et al. (2021), the variability of observed NSWS across Sweden is driven by both changes in large-scale atmospheric circulation and other factors such as modification in surface roughness, which interplay in causing the observed wind variations. This study adds that low-frequency variability can be attributed mainly to variations in atmospheric patterns, which cannot explain the periods of evident stilling: these strong slowdown trends that drive most of the stilling in the last century should be driven by other factors including land-use and cover change and/or anthropogenic aerosol emissions. Future studies should focus on the quantification of the contribution of these factors.