Recent trends in wind speed across Saudi Arabia, 1978–2013: a break in the stilling

We analyse recent trends and variability of observed near‐surface wind speed from 19 stations across Saudi Arabia (SA) for 1978–2013. The raw wind speed data set was subject to a robust homogenization protocol, and the stations were then classified under three categories: (1) coast, (2) inland and (3) mountain stations. The results reveal a statistically significant (p < 0.05) reduction of wind speed of −0.058 m s−1 dec−1 at annual scale across SA, with decreases in winter (−0.100 m s−1 dec−1) and spring (−0.066 m s−1 dec−1) also detected, being non‐significant in summer and autumn. The coast, inland and mountain series showed similar magnitude and significance of the declining trends across all SA series, except for summer when a decoupled variability and opposite trends of wind speed between the coast and inland series (significant declines: −0.101 m s−1 dec−1 and −0.065 m s−1 dec−1, respectively) and the high‐elevation mountain series (significant increase: +0.041 m s−1 dec−1) were observed. Even though wind speed declines dominated across much of the country throughout the year, only a small number of stations showed statistically significant negative trends in summer and autumn. Most interestingly, a break in the stilling was observed in the last 12‐year (2002–2013) period (+0.057 m s−1 dec−1; not significant) compared to the significant slowdown detected in the previous 24‐year (1978–2001) period (−0.089 m s−1 dec−1). This break in the slowdown of winds, even followed by a non‐significant recovery trend, occurred in all seasons (and months) except for some winter months. Atmospheric circulation plays a key role in explaining the variability of winds, with the North Atlantic Oscillation positively affecting the annual wind speed, the Southern Oscillation displaying a significant negative relationship with winds in winter, spring and autumn, and the Eastern Atlantic negatively modulating winds in summer.


Introduction
Multidecadal studies of observed near-surface (i.e. ∼10-m height) wind speed variability and trends have become more frequent in the scientific literature during the last decade; for a comprehensive global review of 148 studies reporting wind speed trends see McVicar et al. (2012). This recent interest in changes of near-surface wind speed is due to the average reduction of wind speed of −0.140 meter per second per decade (m s −1 dec −1 ) in the last 30-50 years detected over mid-latitude terrestrial surfaces (e.g. Australia, China or the United States), particularly in the Northern Hemisphere (McVicar et al., 2012). Roderick et al. (2007) evoked the term 'stilling' to refer to this declining trend. Two major drivers of stilling have been proposed: (1) an upward trend in land-surface roughness (Vautard et al., 2010;Bichet et al., 2012;Wever, 2012); and (2) decadal variability of atmospheric circulation (Lu et al., 2007;Azorin-Molina et al., 2014. This phenomenon is also found to be in line with the widespread reduction in pan evaporation, which has been attributed to stilling (e.g. McVicar et al., 2012). Contrary to this terrestrial weakening of winds, some recent studies (Kim and Paik, 2015;Dunn et al., 2016;Azorin-Molina et al., 2017a) have reported a slight recovery of wind speed after the 2010s, particularly in central and East Asia. Moreover, the complexity of the multidecadal variability of wind speed is revealed by the increasing trends observed over global ocean surfaces (Wentz et al., 2007;Tokinaga and Xie, 2011;Young et al., 2011), or the decoupled trends of winds below and above the boundary layer (Azorin-Molina et al., 2017b), which demonstrates the existence of uncertainties behind the stilling phenomenon. Therefore, there is the need of further assessing and attributing the variability e967 of wind speed over underexplored regions in the global map of McVicar et al. (2012), such as the Arabian Peninsula, because trends and associated causes likely vary spatio-temporarily across Earth.
Previous studies on wind speed changes over the Middle East are scarce. Table 2 in McVicar et al. (2012) summarizes the near-surface wind speed trends reported across neighbouring countries to our target region, with six and seven of the 14 studies conducted in the region exhibiting increases and decreases, respectively, and one of them reporting a non-trending behaviour. This review demonstrated large discrepancies in wind speed changes in space and time, particularly this latter factor as the observed sign and magnitude of trends are especially sensitive to the time-period analysed (Troccoli et al., 2012). For Saudi Arabia (hereafter SA), Rehman (2004) and Rehman et al. (2007) conducted two preliminary studies on the variability of winds for site-specific stations (i.e. Yanbo 24 ∘ N, 38 ∘ E; and Rafha 30 ∘ N, 44 ∘ E) and over different time-spans (i.e. 1970-1983; and 1970-2004), finding a decline of −2.222 m s −1 dec −1 and an increase of +0.040 m s −1 dec −1 , respectively. Recently, Rehman (2013) assessed annual wind speed trends from 1970 to 2006 using 28 stations and found that 18 show declines and the remaining 10 stations increases, with an overall trend of −0.185 m s −1 dec −1 . No other studies have examined wind speed changes across SA, with the majority of the published works performing wind data analysis for wind power potential assessment (e.g. Rehman, 2004;Al-Abbadi, 2005). Over long-term time scales (i.e. >30-years), a better understanding of near-surface wind speed changes across the Middle East is important for assessing: (1) wind-power energy (Al-Abbadi, 2005), (2) hydrological processes including the frequency and severity of droughts (since wind speed modulates atmospheric evaporative demand; El-Nesr et al., 2010) and (3) air pollution dispersal (Lin et al., 2015) or dust emissions (Klingmüller et al., 2016), among many other socio-economic and environmental issues.
Our main objective is to comprehensively assess the spatial and temporal variability and trends of observed near-surface winds over SA, by specifically: (1) applying robust quality control and homogenization protocol to the observations, (2) covering a 36-year period (1978-2013) and (3) detecting trends at annual, seasonal and monthly scales. Our study extends and improves Rehman's (2013) initial study in several ways: the key step (1) above on quality control and homogenization of wind speed data was previously omitted, (2) including more recent 2006-2013 data thus evaluating more recent changes in wind speed in the last decade than was previously performed and (3) not only assessing annual scales as done previously. Our secondary objective is to better understand the role played by atmospheric circulation analysing the influence of teleconnection indices on the observed wind speed changes. Figure 1 shows the geographical setting and topographic features of SA, occupying an area between 15 ∘ and 33 ∘ N, and 34 ∘ and 56 ∘ E. This Middle East country represents the largest territory of the Arabian Peninsula, covering 2.15 million km 2 , and is located between the Red Sea in the west and the Persian Gulf to the east. It is dominated by desert, semi-desert and shrub lands. The territory can be divided into three regions: (1) the central plateau rising abruptly from the Red Sea, with a narrow coastal plain, and highlands, which descend gradually towards the Persian Gulf, (2) the Jabal al-Hejaz Mountains in the southwest, with the highest peak reaching 3133 m a.s.l. in Jabal Sawda, and (3) the southeast sand desert known as the 'Empty Quarter' (647 500 km 2 ), an inhabited land where no wind speed observations are available. Because it is tropical and sub-tropical latitudinal location and the role of the Saudi Arabian high pressure, the climate is 'desert climate' according to Köppen and Geiger, except for the mountain areas. This means that high day time air temperatures, particularly in summer with extremely high values (average of 45 ∘ C), cold nights and large diurnal range air temperatures prevail. Annual precipitation is very low (typically <100 mm year −1 ) associated with few convective events, except for the southwestern region where Indian monsoons brings around 300 mm of rainfall between October and March.

Anemometer observations
Observed wind speed data were supplied by the Presidency of Meteorology and Environment (PME; https://www.pme .gov.sa/; accessed 1 September 2017) of the Kingdom of SA covering an original (i.e. raw) data set, which comprises 29 stations. Wind speed data were measured by different types of anemographs and anemometers, with the RM Young Model 03002 wind sentry anemometer and vane (http://www.youngusa.com/products/11/77.html; accessed 1 September 2017) currently being used to record wind speed and wind direction. The reliability of wind speed data sets used in this study is ensured by the fact that all stations are fully manned (24 h per 7 days per week), which means that periodic inspections were conducted in case of for example anemometer rotor damage or bearing malfunctions (Pindado et al., 2014). Moreover, personal communication with staff from the PME confirmed that calibrations of anemometer devices were carried out once a year, that wind speed measurements were made at 10 m above ground level and all of stations are located at the airports (i.e. well-exposed sites with less proximal environment changes; Azorin-Molina et al., 2014).
The original daily-mean near-surface wind speed series (in m s −1 ) were obtained from 24 hourly mean wind speed values from 0000 till 2400 UTC. Azorin-Molina et al. (2017c) concluded that 24-h wind measurements are more reliable than only four-synoptic times per day for the estimation of long-term wind speed trends. All 29 raw series were used for the robust quality control and   Table 1) of the 29 stations with wind speed data and the main geographic features. Stations are classified into three categories: (1) coast, (2) inland and (3) mountain observation sites.
homogenization protocol described in Section 2.3; however, for the trend analyses we only used those series that corresponded to the longest series with few missing data (i.e. no more than 42 months or 3.5-years) in the 36-year 1978-2013 period. This meant that 10 stations were excluded due to large data gaps (see Figure 1 and Table 1). Moreover, even though nine stations have data availability since 1970, 1976 was unobserved for all stations due to a shutdown ordered by administrative issues. Therefore, to comprehensively assess observed wind speed variability we decided to not infill this year of missing data with reanalysis output to avoid non-real variability because of their shortcomings in simulating near-surface layer processes (McVicar et al., 2008;Pryor et al., 2009;Vautard et al., 2010). This decision set the study period to be the 36-year 1978-2013 span.

Quality control and homogenization
We applied the R package Climatol (Guijarro, 2017; http:// www.climatol.eu/; accessed 1 September 2017), which performs quality control, relative homogenization and missing data infilling on data sets of any climate variable. Climatol works with 'normalized' values, and here the normal ratio was used (i.e. every wind speed series was divided by their mean value). All wind speed data in each series (either present or missing) are estimated as a mean of their nearby data (optionally weighted by an inverse distance function), and differences between observed and estimated values are used for quality control and homogeneity analysis by applying the standard normal homogeneity test (SNHT; Alexandersson, 1986). Outliers greater than a prescribed threshold are deleted and the most inhomogeneous series are split at their more statistically significant break-points (p < 0.05). Thresholds for both outlier rejection and break correction can be adaptively adjusted by assessing histograms of standardized anomalies and maximum SNHT values in all series. This process is repeated iteratively until no statistically significant outliers or break-points greater than thresholds set by the user remain in the series, and a final missing data filling generate the homogenized and complete series. Data completeness is crucial for assessing long-term wind speed trends . The first step to apply this methodology was to aggregate the daily wind speed data into monthly values, allowing a maximum of five missing days in each month as proposed by Azorin-Molina et al. (2014), by means of a function of the same Climatol package; if not, the monthly value was set as missing (441 months were missed and filled for the 19 selected stations). Averaging daily data into monthly values is almost mandatory to improve the detection of inhomogeneities, as the much higher variability of the daily series can compromise this task (Szentimrey, 2013). The homogenization of the monthly wind speed series detected four zero values. After deleting them, a new run resulted in the detection of 76 break-points with SNHT > 25 and the correction of 14 outliers with a departure from their estimated values greater than 4 standard deviations. Figure 2 number of stations that have been found homogeneous (seven stations) or have suffered a different number of splits (up to eight in one station) because of the detection of break-points on them. Figure 2(b) displays an example of break-point detection on the series of anomalies, namely in the Turaif series (id. 27 in Table 1), where the highest SNHT reaches a value of 85 in 1992. Other statistically significant break-points were detected in successive iterations within the program run. Finally, Figure 2(c), top, presents the reconstruction of the series of this station after having split it at the four break-points with SNHT > 25. Original data are plotted in black, and the five homogeneous sub-periods generate complete adjusted series, which are plotted in different colours. The lower part of the figure shows the correction factors applied to the data in each of the adjusted series. After completing this robust quality control and homogenization protocol using all available 29 raw series, 19 homogeneous wind speed series are used in this study; a description of them is given in Table 1 and its spatial distribution in Figure 1.

Atmospheric circulation modes
Following findings from Naizghi and Ouarda (2017), who examined the wind speed variability in the neighbouring United Arab Emirates (UAE) and its relationships with teleconnection indices, we selected five atmospheric circulation patterns to analyse the interplay between large-scale atmospheric circulation and the wind speed variability. These are: (1) the North Atlantic Oscillation (NAO; Jones et al., 1997), as the station-derived normalized pressure between the Gibraltar and Reykjavik as obtained from the Climate Research Unit (https:// crudata.uea.ac.uk/cru/data/nao/; accessed 1 September 2017), (2) the East Atlantic (EA) teleconnection as the north-south dipole of anomaly centre spanning the North Atlantic from east to west, as shown in Barnston and Livezey (1987) and retrieved from the National Oceanic and Atmospheric Administration-National Centers for Environmental Prediction (NOAA-NCEP) (http://www .cpc.ncep.noaa.gov/data/teledoc/ea.shtml; accessed 1 September 2017), (3) the Southern Oscillation (SO) defined as the normalized pressure difference between Tahiti and Darwin as calculated by Ropelewski and Jones (1987)   In all, 29 stations had raw time series yet 10, shaded in grey, did not pass the robust quality control and homogenization protocols and are not used herein.
(Spain) and Padova (Italy) as defined by Martin-Vide and Lopez-Bustins (2006)  We also computed annual SLP composite anomaly maps (from the 1981-2010 mean) from the National Centers for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (NCAR) reanalysis (http:// www.esrl.noaa.gov/psd/; accessed 1 September 2017) for the five strongest positive and negative years of the teleconnection indices having a largest influence on the wind speed changes during 1978-2013. This helps to characterize atmospheric pressure configurations that strengthen or weaken surface winds across SA.

Statistical approach
Monthly homogenized data from each of the 19 stations for 1978-2013 were first expressed as deviations (in m s −1 ) from the 1981-2010 mean. These wind speed anomalies avoid that some windy series dominate the regional series (see below) when averaging different stations. The trend analysis is based in the application of the Sen's slope method (Gilbert, 1987) to retrieve the magnitude (in m s −1 dec −1 ) of the wind speed trend, and a 10-year Gaussian low-pass filter was also computed to illustrate multidecadal variability of wind speed changes. Trends shown in this study were computed for the recent 36-year 1978-2013 study period, and we also reported trends for a 1979-2008 sub-period for comparison purposes against the reference global assessments of Vautard et al. (2010). This also helps to assess the sensitivity of the sign and magnitude of trends to the study period considered (Troccoli et al., 2012;Azorin-Molina et al., 2014;Minola et al., 2016).
The non-parametric correlation coefficient of Mann-Kendall's tau-b (Kendall and Gibbons, 1990) was used for assessing the statistical significance of linear trends at different time-scales: annual, seasonal and monthly. We follow the conventional definition of the four seasons: winter [December-February (DJF)], spring [March-May (MAM)], summer [June-August (JJA)] and autumn [September-November (SON)]. A pre-whitening procedure for removing any statistically significant autocorrelation on the wind speed series was applied to account for the autocorrelation function (von Storch, 1995). Furthermore, we evaluated the statistical significance of wind speed trends at three p-level thresholds, following (1) significant at p < 0.05, (2) significant at p < 0.10 and (3) not significant at p < 0.10. Looking at these three p level grades, instead of at a subjective binary threshold (e.g. at p < 0.05), helps readership to assess wind speed trends from a 'process and importance' perspective, and not only from a 'statistically significant' point of view (Weatherhead et al., 1998;Nicholls, 2001). Furthermore, we also calculated the field significance of the detected statistically significant trends at the 95% confidence level by applying the methods of Livezey and Chen (1983) and Wilks (2006), to evaluate whether the number of stations with statistically significant trends have occurred by chance (Dadaser-Celik and Cengiz, 2014). For determining if and when a change (i.e. break-point year) in wind speed series occurred, we applied the change-point analysis, which is an effective and powerful statistical tool (Gavit et al., 2009). Lastly, the Pearson's correlation coefficient (r) at the three above-mentioned p-level thresholds was computed to quantify the links between the five atmospheric circulation indices and the observed wind speed anomalies.
We not only analysed wind speed trends for each of the 19 stations based on the homogenized data, but also for four averaged anomaly series being: (1) country-average using all 19 stations, (2) mean for the coast region (five stations), (3) mean for the inland region (eight stations) and (4) mean for the mountain region (six stations). These latter three categories refer to a classification of stations as a function of its location along the coastline, in the inland plateaus below 1000 m a.s.l., and above this elevation over mountainous areas. Table 1 describes each of the homogeneous wind speed stations used in this study. Assessing wind speed trends for the three different geographical and topographical environments is crucial since wind speed dynamics differ between them. As examples, the stations near the Red Sea and Persian Gulf are affected by local sea-or land-breeze circulations, the inland stations by local turbulent diurnal-nocturnal winds, and the stations in the Jabal al-Hejaz Mountains mostly by large-scale synoptic winds. Previous works have identified different dynamics e972 C. AZORIN-MOLINA et al.
Units are m s −1 dec −1 . Statistically significant trends were defined as those p < 0.10 (in bold) and p < 0.05 (in bold and in parenthesis). and trends in wind speed depending on geographical settings and local environment. For instance, McVicar et al. (2010) highlighted that wind speed is declining more rapidly at higher elevations than lower elevations, while Azorin-Molina et al. (2017b) found a decoupling of trends below and above the trade-wind inversion layer in the Canary Islands. Previous studies also detected a more statistically significant atmospheric stilling over coastal stations (e.g. Griffin et al., 2010;Minola et al., 2016).  Table 1) located far to the north registered the highest annual mean wind speed with 4.4 m s −1 . Second, the geographical dependence of wind speed, since wind speed is stronger for coastal areas (3.8 m s −1 ; persistent sea breeze circulations) than the inland ones (3.6 m s −1 ; local winds) or than mountainous stations (3.0 m s −1 ; stable layers of the Saudi Arabian subtropical high pressure) located near the Jabal al-Hejaz Mountain range. Seasonally, this latitudinal and geographical dependence of wind speed is also observed for each map: spring is the windiest (3.8 m s −1 ) season, followed closely by summer (3.7 m s −1 ), whereas both winter (3.3 m s −1 ) and autumn (3.1 m s −1 ) wind speeds tend to be much weaker. This same seasonal pattern is also detected for each of the coast, inland and mountain series.

Climatology of wind speed
To summarize the rest of maps shown in Figure 3, the highest variability of wind speed occur in winter and over the inland and mountain series as shown by the standard deviation; mean maximum wind speed values are greater in spring and summer and in the coast and inland stations; whereas, on the contrary, mean minimum wind speed statistics are recorded in winter and autumn, and particularly reflected in the mountain series.

Annual, seasonal and monthly trends
Near-surface wind speeds have significantly declined across SA over 1978-2013 (Table 2). Annual wind speeds reduced by −0.058 m s −1 dec −1 (p < 0.05). The slowdown is more marked for coastal (−0.089 m s −1 dec −1 ; p < 0.05) and inland (−0.074 m s −1 dec −1 ; p < 0.05) stations, and non-significant in the mountain series. Seasonal statistics reveal a temporally dynamic behaviour in the magnitude and statistical significance of the trends throughout the year. Winter experienced the strongest decline, with an all-SA slowdown of −0.100 m s −1 dec −1 (p < 0.05) and negative trends of the same magnitude and significance exhibited for the coast, inland and mountain series. The all-SA spring stilling is of lesser magnitude (−0.066 m s −1 dec −1 ) and statistically significant (p < 0.05), with similar results for the same coast and inland series, whereas non-significant declining trends occurred in mountain series. This dominance of statistically significant all-SA negative wind speed trends was interrupted for summer (−0.042 m s −1 dec −1 , p > 0.10) and autumn (−0.030 m s −1 dec −1 , p > 0.10) seasons. However, when analysing wind speed trends between the three regional series, we found e973 Figure 4. Mean annual and mean seasonal wind speed anomalies (in m s −1 ) for Saudi Arabia, coast, inland and mountain series from 1978 to 2013. A 10-year Gaussian low-pass filter is also shown with a black solid line only for Saudi Arabia to highlight multi-decadal variability. All series are expressed as anomalies from the 1981-2010 mean. Dark grey dashed lines represent the break in the stilling, with numbers showing positive (blue) and negative (red) trends (in m s −1 dec −1 ) before (1978-2001) and after (2002-2013) the break-point for Saudi Arabia; statistically significant trends were defined as those where p < 0.05 (in bold and in parenthesis) and p < 0.10 (in bold).
a characteristic decoupled variability and opposite trends of wind speed between statistically significant declines in the coast (−0.101 m s −1 dec −1 , p < 0.10) and inland series (−0.065 m s −1 dec −1 , p < 0.05) and increase in the mountain series (+0.041 m s −1 dec −1 , p < 0.10). Wind speed changes are negligible in autumn; being negative for all series but not statistically significant (p > 0.10). Also, we detected a clear intra-annual cycle, with declines for all months, but being more negative and statistically significant (p < 0.05) between November and March, and a weak and non-significant decline between April and September; with the exception of a statistically significant slowdown in August. This is the same trend pattern for the coast and inland series, whereas the mountain series displayed an opposite positive sign between April and September. These increasing trends are statistically significant (p < 0.05) in June and August. Figures 4 and 5 represent the temporal variability of annual/seasonal and monthly wind speed series, respectively, for the four series. Annually, the 10-year Gaussian low-pass filter uncovers three phases in the multidecadal variability of wind speed across SA: (1) a non-trending phase from 1978 until around 1990, (2) a marked reduction phase from 1991 until 2001 and (3) a recovery of wind speed phase since 2002 until present. This recent tendency towards positive anomalies after two decades of reductions is a major finding that is also seen in the seasonal and monthly plots. For instance, we calculated trends in wind speed before (  break-point year (at p < 0.05) detected from the annual regional series in 2001, finding an opposite sign in the magnitude of trends between the two periods. This particularly occurs annually, in spring, summer and autumn, and from February to November. On the contrary, winter and December-January showed a consistent declining trend for both the 24-year and 12-year periods. This break in the stilling detected in the last decade is summarized in Figure 6. It is therefore noticeable from this graph the dominance of positive trends in wind for all series during the last 12-year (2002-2013) period compared to the negative trends detected in the previous 24-year  period. However, this recent recovery is only statistically significant in late summer (August) and autumn (September-October) for the all-SA, coast and inland series, and for late spring (May) and summer (July and August) for the mountain series.

Spatiotemporal distribution of trends
Figures 7 and 8 display the spatial distribution of the sign, magnitude and statistical significance of wind speed trends on the annual, seasonal and monthly scales for the 19 stations. Table 3 summarizes the percentage of these stations showing statistically significant positive or negative trends. Annually (Figure 7(a)), wind stilling occurred at 94.7% of stations, with 22.2% (p < 0.05) and 50.0% (p < 0.10) of these stations being statistically significant. It is noticeably that significant negative trends only occurred from 24 ∘ N northward. Seasonally, we found a well-defined cycle of the sign of changes of surface winds throughout the year, with a dominance of declining trends for all seasons, but exhibiting a weak to moderate reduction in the percentage of stations showing negative trends in summer and autumn. In winter (Figure 7(b)), except for the most eastern station, the slowdown of wind dominated over all the rest of locations (94.7%, with 33.3 and 44.4% being statistically significant at p < 0.05 and p < 0.10, respectively), particularly in western SA. In spring (Figure 7(c)), we found similar percentage of stations showing wind stilling (89.5%) but higher proportions of them being statistically significant (35.3 and 52.9% being significant at p < 0.05 and p < 0.10, respectively) and there is a high concentrating in the western and northwestern parts. For summer map (Figure 7(d)), there is a clear change of trends compared to winter and spring, with a smaller percentage of stations exhibiting declines (63.5%; with 25.0% being statistically significant for both at p < 0.05 and p < 0.10, respectively), with significant (at p < 0.05) increases over the Jabal al-Hejaz Mountains in southwestern SA. Finally, in autumn (Figure 7(e)) the percentage of negative trends increase (73.5%), but few of them are statistically significant (21.4 and 28.6% being significant at p < 0.05 and p < 0.10, respectively).
This seasonal spatial trend pattern is also very clear when analysing the monthly maps in Figure 8 and relative frequencies shown in Table 3. For instance, the percentage of stations showing declines greater than 70% from November till May, reaching 100% in January, whereas this wind stilling dropped below the 70% threshold from June till October, with a minimum percentage of stations showing negative trends in August (52.6%) and September (47.4%). Moreover, it is worth noting that while negative trends are statistically significant (p < 0.10) for a moderate proportion of stations in winter, spring and late autumn, positive trends are observed during the warmer months at several stations (particularly, the two mountainous stations in the Jabal al-Hejaz Mountains in June and August). Figure 9 shows the spatial distribution of Pearson's correlation links between the observed wind speed anomalies at the 19 stations (and for the all-SA, coast, inland and mountain series) and the five atmospheric circulation modes at annual and seasonal time scales. Tables 4 and 5 summarize the most correlated teleconnection patterns for each of the four series and the trends of each of the teleconnection indices analysed for 1978-2013, respectively. Overall, we found that atmospheric circulation changes show a statistically significant (at both p < 0.05 and p < 0.10 levels) relationship with the recent variability of wind speeds across SA. As shown in Figure 9 and summarized in Table 4, the North Atlantic Oscillation Index (NAOI) exerted its major influence on winds at annual scale, the Southern Oscillation Index (SOI) modulating surface wind changes in winter, spring and autumn, with the Eastern Atlantic Index (EAI) in summer. The rest of atmospheric circulation modes also partly influenced winds (e.g. the WEMOI annually) and for specific land-stations and regions at different temporal scales, but their relationships with wind speed trends are weaker compared to the NAOI, EAI and SOI. Therefore, we mainly focus our description below on these three most correlated atmospheric circulation indices. Annually the NAOI exerted the strongest and statistically significant link (r, 0.45, p < 0.05) compared to the rest of teleconnection indices. The positive relationship between the NAOI and wind speed occurred across the whole country as shown by significant r values for the  four series. Therefore, the annual stilling of winds could be linked to the negative trend of the NAOI (−0.208, p < 0.05) over the last 3-4 decades (see Table 5); with three of the strongest negative years of the NAOI occurring in the last decade (not shown). It is worth highlighting that the EAI (r − 0.35, p < 0.05) and SOI (r − 0.32, p < 0.10) also exerted a statistically significant but negative relationship on winds annually. In contrast to the NAOI, the EAI (+0.151, p < 0.05) and SOI (+0.263, p < 0.05) showed a tendency towards more positive phases for 1978-2013, with the strongest positive years of the EAI and SOI taking place in the last decade (not shown). This could be partly causing the annual weakening of wind speed.

Relationship with large-scale atmospheric circulation modes
For winter, the SOI exerted a negative relationship on the variability of wind speeds across SA (r −0.53, p < 0.05). The impact of the SOI on the multidecadal variability of winds is particularly strong over the mountain stations (r −0.60, p < 0.05). When looking at the long-term tendency of the winter SOI, we detected a statistically significant (p < 0.05) increase of +0.373 in the last 36-years. Therefore, the strongest wind stilling observed in winter might be associated with this tendency towards a positive phase of the SOI, since the rest of atmospheric circulation modes do not show any significant (p > 0.10) trend. The negative relationship of the SOI on wind speed changes is also statistically significant in spring (r −0.29, p < 0.10) and autumn (r −0.36, p < 0.05). Positive phases of the SOI e977 Figure 8. The same as Figure 7, but at monthly basis. are also becoming more frequent in spring (+0.260, yet not significant at p < 0.10) and autumn (+0.330, p < 0.10), partly contributing to the weakening of winds in both seasons.
Lastly, for summer, the EAI (and to a lesser extent the NAOI) seems to affect the variability of wind speeds, showing a statistically significant negative relationship (r −0.32, p < 0.05), which is particularly strong for the coast (r −0.55, p < 0.05) and the mountain (r 0.44, p < 0.05) regions. Note that this latter positive relationship found between the mountain wind speed series and the EAI in summer has an opposite sign when compared to the e978 C. AZORIN-MOLINA et al. other three series. This decoupled sign of the correlation (also detected for NAOI) might be able to explain the different wind speed trends found in summer between low-elevation coast and inland stations (declines) and the high-elevation mountain stations (increases). Over 1978-2013, the summer EAI (and the NAOI) has become significantly (p < 0.05) positive in +0.327 (negative in −0.406). This means that more positive (negative) phases of the EAI (the NAOI) could have contributed to a weakening of surface winds at low-elevation (coast and e979 Values are expressed as standardized sea level pressure difference. Significant trends were defined as those where p < 0.05 (in bold and in parenthesis) and p < 0.10 (in bold).
inland series), and a strengthening of high-elevation winds (mountain series) as shown in Table 2. SLP anomaly composites for the five strongest positive and negative years of the three climate modes are shown in Figure 10. To summarize, positive anomalies of SLP over the Middle East and negative anomalies over the Western Siberia dominate under positive NAOI and negative SOI phases. This might strength subtropical anticyclonic circulations westward of the region and cyclonic ones northward, thus reinforcing northerly and northwesterly winds across SA, while for the EAI, nearly zero or positive SLP anomalies under negative phases also tend to enhance winds over the region.

Discussion
Recent 1978-2013 trends in annual near-surface wind speed revealed a statistically significant (p < 0.05) decline of −0.058 m s −1 dec −1 across SA, with a clear dominance (94.7%) of stations showing negative trends. This finding is in accordance with the stilling phenomenon, but of lesser magnitude to that quantified by McVicar et al. (2012) in −0.140 m s −1 dec −1 on average over mid-latitude terrestrial surfaces in the last 30-50 years. Moreover, our results also revealed new evidence regarding the wind stilling in an understudied region, for example filling the map shown by Vautard et al. (2010) over the common 1979-2008 period (see linear trends in Table 2), we calculate that wind speed significantly (p < 0.05) declined −0.089 m s −1 dec −1 over all SA then. Our trend results for this sub-period are consistent with Vautard et al. (2010) who found a declining trend of −0.080 m s −1 dec −1 for the neighbouring region of South Asia. However, in general, changes in near-surface wind speed found in here are of lesser magnitude with respect to other studies, particularly, to the statistically significant decreasing trend of −0.185 m s −1 dec −1 found by Rehman (2013)  to our results can be associated with the sensitivity to the different time span (Troccoli et al., 2012) and set of stations used by each study and, most importantly, that most of these studies (e.g. Rehman, 2013) did not apply any quality control and homogenization method to the raw wind speed data set. Using data with high uncertainty and inhomogeneity probably caused overestimating the stilling across SA as Zahradníček and Štèpánek (2017) concluded.
Scientists have been traditionally sceptical about wind speed trends retrieved from anemometer observations (Wan et al., 2010) and therefore studies of long-term variability of surface winds have been scarce until last decade. These data uncertainty are due to climate monitoring issues such as instrumental malfunctions, anemometer changes and technological improvements, shifts in measurement height and sites and so on (Pryor et al., 2009), and the absence of metadata to accurately detect break-points (Li et al., 2011). Recently, some advances on the homogenization of such non-climatic changes in wind speed series have been achieved by the scientific community (Wan et al., 2010;Azorin-Molina et al., 2014, allowing one to make accurate wind speed trend assessments from quality-controlled and homogenized data sets as for SA.
Another interesting feature found in SA, a country located within the dominance of the Saudi Arabian subtropical high pressure is the seasonal-monthly cycle in the magnitude and statistical significance of trends, with significant (p < 0.05) declines in surface wind speed in winter (−0.100 m s −1 dec −1 ) and spring (−0.066 m s −1 dec −1 ), and non-significant in summer and autumn. Therefore, this slowdown in near-surface wind speed was stronger in winter, as also found, e.g. Dadaser-Celik and Cengiz (2014) for Turkey. Our results are in agreement with the recent study of Parvathi et al. (2017) who projected a consistent reduction of the winter monsoon winds over the Arabian Sea at the end of the 21st century using all Coupled Model Intercomparison Project phase-5 (CMIP5) models. This intra-annual cycle is clearly discernible in the percentage of stations showing declines, being widespread in winter (100% in January) yet affecting less stations in summer-autumn (August 52.6% and September 47.4%). This distinct seasonal-monthly trend pattern, with a dominance of negative trends in winter-spring and late autumn months, and non-significant declines or even positive in summer and early autumn have only been recently reported in other mid-latitude regions of the North Hemisphere such as, for example, Spain and Portugal (i.e. between 35 ∘ and e980 C. AZORIN-MOLINA et al.
Here, we also detected a clear decadal variability of surface winds displaying three well-defined phases: (1) a non-trending tendency from 1978 until 1990, (2) a marked wind decline dominating from 1990 until 2001, followed by (3) a recent recovery of wind speed since 2002. This change in the tendency of wind speed observed in SA, particularly in the summer months from May to August, is in agreement with the recent recovery reported by Kim and Paik (2015) for South Korea since 2003. The same was detected in the global assessment conducted by Dunn et al. (2016) and Azorin-Molina et al. (2017a) who observed a global strengthening in terrestrial wind speed, but since 2013 after a continuous slowdown from the beginning of the global HadISD2 anemometer records in 1973 (Dunn et al., 2016). It is noticeable that this recent wind speed recovery has, to date, been poorly analysed as discussed Azorin-Molina et al. (2017a), and therefore specific future attribution studies on this particular break in the stilling or even recovery of winds are strongly needed. For SA, Al Senafi and Anis (2015) found an increase of Shamal days (strong N and NW winds) and dust storms in summer for 1998 till 2012, which agrees with our finding of rebound of wind speed in summer and autumn since 2002. However, in a region such as the Arabian Gulf where local sea-breeze circulations occur in more than 70% of the days in a year (Eager et al., 2008), attribution analysis will be complex due to interplaying local factors and therefore is beyond the scope of the current spatio-temporal assessment of wind speed variability. Future work should fill this knowledge gap in the recent rebound of winds.
Because different atmospheric dynamics dominate coast, inland and mountain geographical regions, we also analysed wind speed trends averaging anomalies series for stations classified into these three categories. Contrary to McVicar et al. (2010), who concluded that near-surface wind speeds are declining more rapidly at higher elevations than lower elevations, we reported a similar magnitude of negative trends for the three regions, being even less negative for the mountain series (e.g. in spring and autumn). More interestingly, we found a decoupled variability and opposite trends of wind speed in summer between the low-elevation coast and inland series, which display statistically significant declines (at p < 0.10 for coast and p < 0.05 for inland), and the high-elevation mountain series showing significant increases (p < 0.10). Recently, this finding has been noticed by Azorin-Molina et al. (2017b), who observed a decoupling wind speed tendency when comparing stations below and above the trade-wind inversion layer in the Canary Islands (Spain; between 27.6 ∘ and 29.4 ∘ N, and 13.3 ∘ and 18.2 ∘ W), with a dominant positive tendency of wind speed observed in the mountainous station. Azorin-Molina et al. (2017b) attributed these elevation-based differences to different altitude-dependent atmospheric dynamics. Therefore, our results contrast McVicar et al. (2010), who linked the strong declining trend in high-elevation stations to the observed increasing altitude of the tropopause (Santer et al., 2003). These few studies, and limited knowledge, suggest the need for comprehensive assessments and attribution of wind speed variability at different elevations (i.e. mountains) and at different altitudes (i.e. levels of the troposphere), as recently suggested by Azorin-Molina et al. (2017b). Another interesting feature was the observation of a latitudinal dependence of the sign, magnitude and significance of wind speed trends, with a clear dominance of significant (negative) trends in the northern half part of SA, with non-significant changes in the southern half except for the positive trends detected over the Jabal al-Hejaz Mountains. Moreover, a characteristic feature when analysing the spatial distribution of trends is that changes can be different between nearby stations (Figures 7 and 8), which is a consequence of complex mechanisms driving trends or, even the homogenization applied, derived from instrumental and observational issues in measuring wind speed.
Our analyses of the relationship of the observed wind speed anomalies with five teleconnection indices showed close agreement with the findings of Naizghi and Ouarda (2017), who analysed links between the long-term wind speed variability in the neighbouring United Arab Emirates and teleconnection indices. For instance, we detected that atmospheric circulation has important relationship with the variability of terrestrial near-surface wind speed trends across SA, with the NAOI strongly affecting winds annually, the SOI in winter, spring and autumn, and the EAI (also NAOI) in summer. Naizghi and Ouarda (2017) found similar associations with the NAOI and EAI modulating wind speed in summer, and the ENSO and the Indian Ocean Dipole in winter and autumn. Dasari et al. (2018) also found an impact of the ENSO on the interannual variability of the Red Sea convergence zone and associated rainfall, denoting an increase of southeasterly winds (20-30%) over the southern Red Sea and a decrease of northwesterly winds (10-15%) over the northern Red Sea under El Niño phases. Trends shown by atmospheric circulation modes might be associated with the observed trends in wind speed. For instance, the statistically significant 1978-2013 annual stilling could be partly explained by the negative trend of the NAOI. The same association for the SOI in winter, spring and autumn, which showed a tendency towards more positive phases in the last 36-years and might be exerting to the observed e982 C. AZORIN-MOLINA et al. weakening of winds in both seasons. Lastly, the summer EAI has tended to become more positive, contributing to the slowdown of winds observed at low-elevations and a recovery of winds measured at the high-elevation mountain stations. It is worth highlighting that this positive relationship between the high-elevation mountain wind speed series and the EAI is of opposite sign compared to the rest of regional (coast and inland) series (see Figure 9). This decoupled sign of the correlation (also detected for NAOI) might be explaining the different wind speed trends found in summer between stations located in low-elevation stations (declines) and the high-elevation stations (increases) (see Table 2). This is also shown in the different spatio-temporal correlations (Figure 9) found between the teleconnection indices and wind speed series for coast, inland and mountain stations, denoting the different atmospheric dynamics dominating these geographical regions. Moreover, we found that positive (i.e. anticyclonic) SLP anomalies tend to strengthen winds, which is in agreement with Yu et al. (2015) who concluded that anomalously anticyclonic circulation over the Arabian Peninsula enhances Shamal winds (northwesterly winds). Other previous studies also highlighted the major role played by changes in atmospheric circulation on the variability of wind speeds: e.g. Dadaser-Celik and Cengiz (2014) in Turkey, Azorin-Molina et al. (2014 in Portugal and Spain; to name but a few. Therefore, it is obvious that atmospheric circulation plays a key major role in explaining the multidecadal variability of surface winds, estimated in 10-50% by Vautard et al. (2010) and associated with a weakening in the inter-hemispheric SLP gradient force (Guo et al., 2011;Parvathi et al., 2017), but other factors partly contribute to the observed trends, including changes in surface roughness (Wu et al., 2016), instrumentation (Wan et al., 2010), in particular the degradation of anemometer bearings in a dusty environment such as SA (Azorin-Molina et al., 2018), different time intervals (Azorin-Molina et al., 2017c), among many other factors contributing simultaneously and varying spatio-temporally. Therefore, a comprehensive attribution analysis of the complex relationships of for example atmospheric circulation (teleconnection indices) on wind speed variability and the recent recovery is strongly needed.
Lastly, SA is planning to diversify its national energy mix by supplementing the energy production through renewable sources of energy like wind, solar photovoltaic, geothermal and waste-to-energy (Al-Maamary et al., 2017). Hence, the present study will help in potential site selection for wind power deployment projects and accurate wind power assessment process (Rehman et al., 2007).

Conclusion
The major findings of the 1978-2013 surface winds in Saudi Arabia (SA) are: 1. Wind speed over the country declined significantly −0.058 m s −1 dec −1 (p < 0.05) at the annual scale, which was contributed by significant declines in winter (−0.100 m s −1 dec −1 ) and spring (−0.066 m s −1 dec −1 ). The strong trend in winter is particularly influenced by declines in February and March. The trends in the northern SA are particularly statistically significant. 2. For summer, a decoupled variability and opposite trends of wind speed were detected between the low-elevation coast and inland series (significant declines: −0.101 m s −1 dec −1 and −0.065 m s −1 dec −1 , respectively), and the high-elevation mountain series (significant increase: +0.041 m s −1 dec −1 ). 3. A clear intra-annual cycle in the number (percentage) of stations showing declines was found, with the highest % in winter (100% in January) and the lowest in summer-autumn (August 52.6% and September 47.4%). 4. The wind stilling dominated from 1978 until 2001, followed with a recent break in the stilling or even rebound associated with the detected strengthening of winds in summer and autumn for 2002-2013. 5. Atmospheric circulation variability plays a key role in explaining the interannual variability of the surface winds, reflected by the statistically significant correlations between the annual wind speed and NAOI (positive relationship), winter-spring-autumn SOI (negative relationship) and summer EAI (negative relationship).