Estimating centennial-scale changes in global terrestrial near-surface wind speed based on CMIP6 GCMs

A global terrestrial stilling in recent decades has been reported, but the centennial-scale changes in global terrestrial near-surface wind speed (NSWS) and the potential contributing factors are yet to be revealed. Consequently, in this study, centennial-scale changes in global terrestrial NSWS are investigated based on Coupled Model Intercomparison Project phase 6 datasets, and that the potential factors causing those changes are detected. The results show that the global annual mean NSWS increased from 1850 to 1967 (+0.0045 m s−1 decade−1, p< 0.01), with significant increases in North America, Europe, Africa, and South Asia. However, the NSWS decreased from 1968 to 2014 (−0.0044 m s−1 decade−1, p < 0.01), significantly so in the mid-to-high latitudes of the Northern Hemisphere. The seasonal mean NSWS also increased before the 1960s and decreased thereafter. However, the NSWS over South America and most of Southern Africa increased during the study period. The changes in NSWS were caused mainly by changes in the number of strong windy days. The increase in NSWS from 1850 to 1967 could be attributed to internal variability, and the decrease in NSWS from 1968 to 2014 could be attributed to natural, aerosol, and greenhouse-gas forcings. However, internal variability acted mainly to increase the NSWS from 1968 to 2014, and so it is suggested that the contributions of external forcings to the global terrestrial stilling after the 1960s were considerable.


Introduction
Investigating changes in near-surface wind speed (NSWS) facilitates the understanding of atmospheric circulation, improves climate analysis and prediction. As more countries commit to emissions reductions by mid-century to curb anthropogenic climate change, decarbonization of the electricity sector becomes a first-order task in reaching this goal. Renewable, particularly wind power, will be predominant component of this transition (Sherman et al 2017(Sherman et al , 2021. A rapid increase in wind power capacity has been predicted in many countries and regions (Pryor et al 2020). The NSWS changes have a considerable influence on the wind energy resource (He et al 2010). NSWS also shows effects on global and regional evapotranspiration (Roderick et al 2007, Liu et al 2014, visibility (Sun et al 2018, Zhang et al 2020a, agriculture, and ecosystems (Deng et al 2021). The intensification of NSWS may exacerbate soil erosion, thereby generating more-severe dust storms (Alizadeh-Choobari et al 2014, Guan et al 2017, Segovia et al 2017, Zhang et al 2019a, 2020b. Consequently, investigating the long-term changes in NSWS is beneficial to the assessment of wind energy and environmental and ecological issues related to NSWS.
Previous studies have reported long-term decrease in global NSWS (Vautard et al 2010, Zhang et al 2019b. The reduction in NSWS has also been observed at regional scale, including in China (Guo et al 2011, Lin et al 2013, 2021a, Australia (McVicar et al 2012), Czech Republic (Brazdil et al 2009), Switzerland (McVicar et al 2010, United Kingdom (Earl et al 2013), Turkey (Dadaser-Celik and Cengiz 2014), South Korea (Kim and Paik 2015), Spain andPortugal (Azorin-Molina et al 2014, 2016), and Brazil (Gilliland and Keim 2018). Roderick et al (2007) referred to these declining trends in NSWS as 'stilling' . The long-term reduction in NSWS is also revealed by reanalyses, but which is weaker in reanalyses than that in observations (You et al 2010, Zhang and. Recent studies mainly analyzed the climatology, interannual variability, and trends of NSWS at global and regional scales (Torralba et al 2017, Ramon et al 2019, Wohland et al 2019, and the spatiotemporal characteristics of global NSWS at centennial scale are yet to be revealed. Some studies also discovered a weak increase in NSWS over the past decade. Zeng et al (2019) referred to this phenomenon as 'reversal' . The reversal of stilling has also been discovered at regional scale. For instance, the annual mean NSWS increased over South Korea after 2003 (Kim and Paik 2015), southwestern China after 2000 (Yang et al 2012), northwestern China after 1993 , and eastern China after 2010 (Zha et al 2021b). Zeng et al (2019) highlighted the major role of ocean-atmosphere oscillations in explaining stilling versus reversal. We proposed that stilling and reversal could be a cyclical, decadal pattern of NSWS at the centennial-scale (Shen et al 2021).
Compared to the investigation of historical NSWS, the projections of NSWS are investigated rarely at global scale. Based on the Coupled Model Intercomparison Project Phase 5 (CMIP5), the wind speed at the end of this century are projected to reduce in most regions of globe (Chen et al 2012, Kumar et al 2015. Karnauskas et al (2018) pointed out that the wind power reduction in future was located mainly in the mid-latitudes across the Northern Hemisphere (NH) and increase over the tropics and Southern Hemisphere (SH). Based on a multi-model ensemble of CMIP5, Sherman et al (2021) suggested that available wind resource over China and India will decline slightly. However, what are the main features of NSWS at centennial-scale during the historical period? The investigation of historical wind speed changes at centennial-scale could provide a scientific basis for the projections of wind speed.
The NSWS changes are influenced by internal variability and external forcings (Wu et al 2018a).
Internal variability affecting NSWS are attributed to the changes of large-scale ocean-atmosphere circulations. For instance, the NSWS changes over the NH could be influenced by the Tropical Northern Atlantic (TNA), the Pacific Decadal Oscillation (PDO), and the North Atlantic Oscillation (NAO) (Zeng et al 2019). El Niño/Southern Oscillation (ENSO) is a result of sea-air interaction on multiple space-time scales, which has a significant effect on the Asian monsoons and PDO (Ni et al 1995, Fu et al 2011, Lu et al 2016. We have revealed that the ENSO has a considerable effect on NSWS in China (Shen et al 2021). External forcings impact on NSWS is attributed to anthropogenic activities, including land use and cover change , Zhang et al 2019b, anthropogenic aerosol emissions , and greenhouse gas (GHG) emissions (Jiang et al 2010, Zha et al 2020. Deng et al (2021) analyzed the global NSWS based on CMIP6 and suggested that the NSWS in NH decreased from 1980 to 2010, and that the attribution analysis was focused mainly on revealing the differences between land and ocean NSWSs. Unfortunately, the centennial-scale changes in terrestrial NSWS and the potential factors contributing to those changes were not investigated.
Consequently, the novelty of this study compared to previous ones is twofold: (a) spatiotemporal characteristics of the centennial-scale changes in global terrestrial NSWS are investigated based on CMIP6 models; (b) potential factors causing decadal changes in NSWS at centennial scale are detected. The results presented herein help to improve the understanding of global terrestrial NSWS changes.

Datasets
Multi-model simulations of NSWS from CMIP6 are used. The simulations are distinguished by their 'ripf ' index, which denotes the initial states (r), initialization methods (i), physics versions (p), and the forcing datasets (f) for CMIP6 (Grose et al 2020, Jian et al 2020, Parding et al 2020. This study uses the first ensemble member ('r1i1p1f1') each from a group of 33 models that were available at the time of writing (table S1). CMIP6 offers improved spatial resolution and physical parameterization compared to CMIP5 (Eyring et al 2016(Eyring et al , 2019. To reveal the factors inducing centennial-scale variability in NSWS, the Detection and Attribution Model Intercomparison Project (DAMIP) is used. Only ten models of DAMIP were available at the time of writing (table S1). Based on DAMIP, four main forcings affecting NSWS changes are discussed: (a) aerosol forcing (anthropogenicaerosol-only historical simulations (BC, OC, SO 2 , SO 4 , NO X , NH 3 , CO, NMVOC): hist-aer experiment in DAMIP), (b) GHG forcing (well-mixed GHGonly historical simulations: hist-GHG experiment in DAMIP), (c) natural forcing (natural-only historical simulations (solar irradiance, stratospheric aerosol): hist-nat experiment in DAMIP), and (d) internal variability (Historical simulation experiment minus main external forcing experiment (hist-GHG, histaer, hist-CO 2 )). It is worth noting that the internal variability is assuming linearity of the GHG, aerosol and natural variability responses. Details for DAMIP are given in Gillett et al (2016).
Four centennial-scale reanalyses with daily mean wind speed are used to compare with the CMIP6: (a) the European Centre for Medium-Range Weather Forecasts (ECMWF) 20th-century reanalysis assimilating surface observations only (ERA-20C) that is produced with ECMWF Integrated Forecasting System version Cy38r1 (Poli et al 2013), (b) the ECMWF ten-member ensemble of coupled climate reanalyses of the 20th century (CERA-20C) that is produced with Integrated Forecasting System version Cy41r2 (Hersbach et al 2015), (c) the National Oceanic and Atmospheric Administration (NOAA)-Cooperative Institute for Research in Environmental Sciences (CIRES) 20th-century reanalysis version 2 (NOAA-CIRES-20CR-V2C) that is produced by the Earth System Research Laboratory Physical Sciences Division from NOAA and the University of Colorado CIRES (Compo et al 2011), and (d) NOAA-CIRES-DOE 20th-century reanalysis version 3 (NOAA-CIRES-20CRV3) that is the new version of the 20th-century reanalysis (Slivinski et al 2019). Details for these reanalyses are presented in table S2.

Methods
A piecewise linear function (PWLF) is used to fit the trends of NSWS during different periods. PWLF automatically detect the optimal turning point (TP) and allowing multiple linear models to be fitted to each distinct section of the time series (Fyllas et al 2009, Jekel andVenter 2019). A PWLF can be described as follows where b 1 is the location of the first breakpoint, b 2 is the location of the second breakpoint, and so forth until the last breakpoint b n b . α 1 , α 2 , · · · , α n b −1 and β 1 , β 2 , · · · , β n b −1 are the regression coefficients. In cases where the TPs are unknown, global optimization is used to find the best set of TPs that minimizes the overall sum-of-square of the residuals (Storn and Price 1997). Multiple TPs can be set; however, in order to avoid overfitting, we assume that there is one and only one TP during the entire study period. The necessity of introducing the TP is examined with the significance t-test under the null hypothesis that 'β is not different from zero' . If the maximum of the absolute value for β is occurred in a given year, that year is determined to be the TP. Statistical significance for the regression include the goodness-of-fit, the P value for the entire model, and the trends before and after the TP (Zha et al 2021b).
The variance of wind speed is calculated based on annual mean value of NSWS (equation (2)) where V x denote the variance of NSWS, n denote the total time period, x i denote the NSWS, andx denote the mean value of NSWS. We calculate the multimodel ensemble (MME) mean of NSWS firstly and then calculate the variance. A Gaussian low-pass filter is used to determine the decadal signals of NSWS with a nine-year time window and an empirical orthogonal function (EOF) is used to partition the wind fields into a set of orthogonal modes consisting of spatial structures and corresponding time series (Wu et al , 2018b(Wu et al , 2018c. All datasets are interpolated onto the 1.0 • × 1.0 • resolution grid using bilinear interpolation. To evaluate the frequency trends for different categories of windy days, four windy days are defined based on their wind speed percentiles: (a) light windy days (<25th percentile), (b) gentle windy days (between 25th and 50th percentiles), (c) moderate windy days (between 50th and 75th percentiles), and (d) strong windy days (>75th percentile) (Zha et al 2021b). The trends of annual/seasonal NSWS and frequency trends for different categories of windy days are calculated based on the least-square method. We calculate the ensemble mean of multi-models for NSWS firstly and then calculate the trends and variances . Four seasons are defined as boreal winter (December, January, and February (DJF)), spring (March, April, and May (MAM)), summer (June, July, and August (JJA)), and autumn (September, October, and November (SON)).

Spatiotemporal characteristics of global NSWS at centennial scale
Temporal evolutions of NSWS in the MME mean at centennial scale are shown in figure 1. Annual mean NSWS increased from 1850 to 1967 (+0.0045 m s −1 decade −1 , p < 0.01) and decreased thereafter (−0.0044 m s −1 decade −1 , p < 0.01) ( figure 1(a)). Seasonally, the NSWS also increased before 1967 and decreased after, with the strongest increase in winter (+0.0060 m s −1 decade −1 , p < 0.01) ( figure 1(b)) and the strongest decrease in spring (−0.0063 m s −1 decade −1 , p < 0.01) ( figure 1(c)). The TPs were inconsistent among the seasons but mainly located in and around the 1960s. Previous studies have reported a decrease in global NSWS since 1960s (Vautard et al 2010, Zhang et al 2019b. The spatial patterns of trend changes among models were also compared (figures S1 and S2 (available online at stacks.iop.org/ERL/16/084039/ mmedia)), and most models also mainly exhibited an increase in NSWS from 1850 to 1967 (figure S1) and a decrease in NSWS from 1968 to 2014 (figure S2). However, the differences of trends among models are considerable.
The EOF was used to assess whether the aforementioned centennial-scale changes were the predominate signals in NSWS (figure S3). The spatial pattern of the first mode of the EOF mainly exhibited positive values across globe from 1850 to 1967, and the time series of the first mode exhibited an increase with the explained variance of 32.86%, which passed the significance North test (figures S3(a) and S3(b)). From 1968 to 2014, the spatial pattern of the first mode of the EOF mainly exhibited negative values over the mid-to-high latitudes of NH and positive values over SH; furthermore, the time series also exhibited an increase with the explained variance of 18.57%, which passed the significance North test (figures S3(c) and S3(d)). Consequently, the NSWS changes from increasing to decreasing should be the predominant signals in NSWS changes at centennial scale. To demonstrate whether the aforementioned characteristics of NSWS based on the CMIP6 models can be produced in global reanalysis products, NSWS changes in four analyses are also analyzed ( figure S4). The NSWS showed a characteristic from increasing to decreasing in the selected four centennial-scale reanalyses, although the interannual variability, trends and TPs were inconsistent with those from CMIP6.
The mean values and variances of NSWS from 1850 to 1967 and from 1968 to 2014 displayed similar spatial pattern. There was no distinct difference in NSWS climatology during the two periods. The large values of NSWS mainly located in North America, North Africa, Central Asia, and Australia, and the small values of NSWS mainly located in South America, Southern Africa, and East Asia (figures 2(a) and (b)). The large values of variance mainly located in eastern part of North America, Europe, and Central Asia (figures 2(c) and (d)), implied that the NSWS showed strong interannual and interdecadal changes in these regions. A significant difference was found in the NSWS trends during the two periods. The NSWS mainly increased from 1850 to 1967, with significant increases over North America, Europe, Africa, and South Asia (figure 2(e)). From 1968 to 2014, the NSWS decreased over the NH (e.g. in North America, Europe, Central Asia, and South Asia) and increased over the SH (e.g. in South America and Southern Africa) (figure 2(f)). The significant differences of seasonal NSWS changes were also found in the trends (figure 3); however, the mean values and variances of seasonal NSWS showed similar spatial patterns during the two periods (figures S5 and S6). From 1850 to 1967, the global NSWS tended to increase in each season (figures 3(a), (c), (e) and (g)); however, it mainly decreased over the NH and increased over the SH from 1968 to 2014 (figures 3(b), (d), (f) and (h)). These characteristics were more significant in summer (figure 3(f)) and autumn (figure 3(h)) than other seasons.
Overall, the changes in global NSWS generated a transition since 1960s that occurred mainly at midto-high latitudes in the NH. In the SH, particularly in South America and Southern Africa, the NSWS has been increasing in all seasons. The global mean NSWS increased before 1960s and decreased after, which was the primary features of NSWS changes. However, the NSWS did not change uniformly throughout world. The global mean NSWS changes from increasing to decreasing at centennial scale did not mean that this characteristic presented everywhere. The NSWS changes showed regional difference. , and autumn (SON), respectively. Different windy days are defined as follows: light windy days (<25th percentile), gentle windy days (between 25th and 50th percentiles), moderate windy days (between 50th and 75th percentiles), and strong windy days (>75th percentile). * , * * , and * * * indicate that the trends pass the significance t-test at the 0.10, 0.01, and 0.001 levels, respectively. The black bars denote the uncertainties of the results.

Trends for different categories of windy days during the two periods at centennial scale
The decadal changes in NSWS at centennial scale could cause the changes in frequencies of different windy days, so the frequency trends for different windy days are investigated (figure 4). From 1850 to 1967, the number of strong windy days increased (+3.54 days decade −1 ; p < 0.001), while the numbers of light and moderate windy days decreased, at the rates of −3.00 (p < 0.001) and −1.98 days decade −1 (p < 0.001), respectively. The number of strong windy days also increased in all seasons from 1850 to 1967, with the strongest increase in winter (+2.06 days decade −1 ; p < 0.001) and the weakest increase in autumn (+1.08 days decade −1 ; p < 0.001). Compared to figure 1, the significant increases in annual and seasonal mean NSWSs from 1850 to 1967 were due mainly to the increase in the number of strong windy days and decrease in the number of light windy days. From 1968 to 2014, the numbers of all categories of windy days decreased except for light windy days (+3.97 days decade −1 ; p < 0.001) ( figure 4(a)). Seasonally, the number of light windy days increased, with the most significant increase in summer (+2.96 days decade −1 ; p < 0.01) ( figure 4(d)). The numbers of moderate and strong windy days in all seasons decreased except for strong windy days in winter, with the most significant decreases in spring (−1.20 days decade −1 ; p < 0.1) for moderate windy days and summer (−1.94 days decade −1 ; p < 0.1) for strong windy days. Accordingly, the pronounced decrease in NSWSs from 1968 to 2014 were due mainly to the decreases in the numbers of strong and moderate windy days.

Potential factors causing centennial-scale changes in NSWS
The potential factors causing the centennial-scale changes in NSWS are detected based on the ten DAMIP models of CMIP6. Actually, these models can also capture the characteristics of NSWS changes during the study period, i.e. increasing before 1960s and decreasing after (figure S7).
Temporal evolutions of NSWS in different forcings are shown in figure 5. In the experiment with internal variability, the NSWS increased significantly from1850 to 1967, consistent with the historical NSWS changes. With aerosol and natural forcings, the NSWS neither increased nor decreased significantly from 1850 to 1967. With GHG forcing, the NSWS decreased during the period ( figure 5(a)). However, with only aerosol forcing experiment, the NSWS decreased in the MME of DAMIP after the 1960s. The NSWS trends among different models from 1850 to 1967 and from 1968 to 2014 are also  (c)). The increase in historical NSWS from 1850 to 1967 was due mainly to internal variability, which was discovered in nine of the ten DAMIP models, and no other forcing gave rise to pronounced NSWS strengthening from 1850 to 1967 ( figure 5(b)). However, from 1968 to 2014, the decrease in historical NSWS was due mainly to natural, aerosol, and GHG forcings; for aerosol forcing in particular, most of the DAMIP models displayed decreased NSWS ( figure 5(c)).
How the different forcings affect the probability density function (PDF) of NSWS is investigated ( figure 6). The PDF of NSWS anomalies with internal variability was closer to that of the historical NSWS during the two periods compared to the other forcings, but the PDF moved toward larger positive NSWS anomalies after the 1960s. Consequently, it was mainly internal variability that increased NSWS after the 1960s. During the two periods, the NSWS PDFs with aerosol, natural, and GHG forcings were higher and narrower than the historical NSWS PDF. From 1850 to 1967, the mean values of NSWS anomalies with aerosol, GHG, and natural forcings were close to zero, but not so with internal variability. By contrast, after the 1960s, the mean values of NSWS anomalies with aerosol and GHG forcings were negative, whereas the contributions of internal variability to NSWS changes were positive. Therefore, external forcings (aerosol and GHG forcings) decreased the NSWS after the 1960s.
Spatial patterns of NSWS trends under different forcings are also analyzed ( figure 7). From 1850 to 1967, accompanied by GHG, natural, and aerosol forcings, the NSWS mainly decreased, strongly so in North America and Eurasia; however, the increase in NSWS from 1850 to 1967 was due mainly to internal variability ( figure 7(g)). Compared to figure 2(e), the increase in global NSWS from 1850 to 1967 was due mainly to internal variability. From 1968 to 2014, the NSWS also mainly decreased accompanied by GHG, natural, and aerosol forcings, but the decrease in NSWS accompanied by these factors from 1968 to 2014 was more pronounced than that from 1850 to 1967. Note that the increase in NSWS from 1968 to 2014 was due mainly to internal variability, especially in East Asia, South America, and Africa. Therefore, the NSWS would decrease more significantly from 1968 to 2014 if the effects of internal variability were excluded. Compared figure 7 with figure 2(f), the reduction of NSWS in figure 7 was weaker than that in figure 2(f). Consequently, we proposed that the GHG, natural, and aerosol forcings could induce the NSWS reduction from 1968 to 2014; however, there was not any single factor that significantly dominate the reduction in NSWS. It is worth noting that aerosol and GHG affect the NSWS change at the centennialscale could show uncertainty due to the internal variability cannot be well excluded based on the MME of CMIP6. Furthermore, the uncertainty of aerosol and GHG influence on NSWS may also be generated from the model bias (Bichet et al 2012, Deng et al 2021.

Conclusions and discussions
The centennial-scale variability in global terrestrial NSWS is discussed, and the potential factors causing centennial-scale changes in NSWS were revealed based on CMIP6 global climate models in this study. The main results are as follows.
(a) From 1850 to 1967, the global annual mean terrestrial NSWS increased, and the significant increases were found in North America, Europe, Africa, and South Asia. From 1968 to 2014, the global annual mean terrestrial NSWS decreased, and the significant decreases were found in North America, Europe, Central Asia, and South Asia. The mean values and variances of NSWS during the two periods showed similar spatial patterns, and therefore there was no distinct difference in NSWS climatology during the two periods. The significant difference of global terrestrial NSWS changes at the centennial scale was shown in the NSWS trends. Global annual mean and seasonal mean NSWSs increased were caused by the increase in the number of strong windy days from 1850 to 1967, and which decreased were caused by the decreases in the numbers of moderate and strong windy days from 1968 to 2014. (b) The increase in historical NSWS from 1850 to 1967 was due mainly to internal variability, which was found in nine of the ten DAMIP models, and no other forcing gave rise to pronounced NSWS strengthening from 1850 to 1967. From 1968 to 2014, the decrease in historical NSWS was due mainly to natural, aerosol, and GHG forcings. The NSWS would decrease more significantly from 1968 to 2014 if the effects of internal variability were excluded. The NSWS PDFs with aerosol, natural, and GHG forcings were higher and narrower than the historical NSWS PDF.
In this study, the results show that the aerosols could decrease NSWS. Some studies have proposed that aerosols can induce the NSWS changes. On a continental scale, aerosols reduce surface insolation and weaken the land-ocean thermal contrast, these inhibiting the development of monsoons. Locally, aerosol radiative effects alter the thermodynamic stability and convective potential of the lower atmosphere leading to reduced temperatures and increased atmospheric stability . An increase in air stability due to interactions between aerosols and radiation can reduce vertical mixing, which in turn declines the vertical flux of horizontal momentum. Since winds are generally higher aloft than at the surface, weakened vertical mixing reduces the transfer of fast winds aloft to the surface, slowing surface winds compared with those aloft and atmospheric circulations (Jacobson and Kaufman 2006). Therefore, the effects of aerosols on NSWS could be attributed to the impacts of aerosols on monsoon circulations and thermodynamic stability. The Asian monsoon region is a primary source of emissions of diverse species of aerosols; therefore, the NSWS reduction under the aerosol forcing is stronger in East Asia and South Asia than that in other regions. However, the physical mechanisms of aerosols affect NSWS and monsoons are complex, and which should be investigated systematically in future. The GHG forcing could affect the global NSWS changes through modulating the meridional atmospheric circulation (Deng et al 2021).
Several theories have tried to identify potential mechanisms describing how internal variability influence wind speed (Wang 2002, Timmermann et al 2018. For instance, the positive phase of TNA is linked with a weakened Hadley circulation and leads to a southward component of surface wind and a reduction of wind speed in the mid-latitudes (Zeng et al 2019). The temperature gradients during the negative and positive phases of PDO generate the easterly and westerly components of NSWS, which decrease and increase the prevailing westerly winds in the mid-latitudes, respectively . The negative and positive phases of NAO have different jet stream configurations and wind systems over Europe. Strong and weak East Asia summer monsoon are associated with increases and decreases in the surface air temperature difference between land and sea, respectively, and then which could cause the changes in NSWS (Huang et al 2018). Therefore, the NSWS changes over the monsoon regions could be influenced by the monsoon circulations (Xu et al 2006). However, the interaction and modulation among different internal variabilities are pronounced; therefore, the NSWS changes could not be simply linked with just one internal variability but should be determined by the combined effects of variations in multiple internal variabilities (Zeng et al 2019). In this study, we just detect that the internal variability could influence the NSWS; however, we cannot guarantee that all internal variability signals satisfy a linear relationship due to the different internal variability signals could not be independent. The mechanism and contributions of internal variability affects NSWS should be performed in the future.
The decreases of NSWS since the 1960s in CMIP6 and the four global reanalysis products are weaker than those actually observed. Previous studies have discussed possible causes of this difference. Including, the decrease in observed NSWS is a manifestation of changes in surface roughness that are not included in the surface boundary conditions used in the climate models (Chen et al 2012), the current models have only a relatively weak capacity for representing some aspects of atmospheric flow (Zeng et al 2019), the observed trends are partly the product of non-climate-related factors, such as inhomogeneities in station settings or instrumentation (Zha et al 2021a), the data assimilation systems of reanalyses could contain inappropriate model topography and inaccuracies in the atmospheric boundary layer processes (Fan et al 2021), and the spatial resolution of the model is too coarse (Shen et al 2021).
The results of the present study offer new insights or at least a perspective on the changes of global terrestrial NSWS, and which could provide a scientific basis for the projections of wind speed and wind energy. However, some limitations must be mentioned. We determined potential factors inducing centennial-scale changes in NSWS, but we did not quantify their contributions and reveal the physical mechanisms of different forcings affect the centennial-scale variability in NSWS. These aspects could be carried out based on the large ensembles with different Earth system models (ESMs) due to the comparisons with large ensembles of different ESMs can reduce the uncertainty of results. Additionally, the different forcings influencing the NSWS show regional differences, and the predominant forcing factors causing regional NSWS should be quantified based on dynamical downscaling in the future.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// esgf-node.llnl.gov/projects/cmip6/. Assistant Project of Chinese Academy of Sciences, the Program for Key laboratory in University of Yunnan Province, and the Chinese Jiangsu Collaborative Innovation Center for Climate Change.